常微分方程

Ordinary differential equation，簡稱ODE，自變量只有一個的微分方程。
例子1： $\frac{dy}{dx}=f(x,y)$ ,$f(x,y)$是已知函數

偏微分方程

Partial differential equation，簡稱PDE，自變量有多個的微分方程。
例子2：$u_t?a^2{u}_{xx}=0,a>0$為常數（熱傳導方程，拋物型方程的典型代表）

顯式（Explicit）

第n步結果可以從n-1, n-2, …1步的結果直接推導出來，迭代時每步的計算量很小，但迭代增量也有限制，不能太大，否則會出現發散。

隱式（Implicit）

第n步的計算結果不能直接從前面的結果推導出來，必須做進一步的求解，這樣，迭代時每步的計算量很大。

泰勒公式

$$f(x)=\frac{f(a)}{0!}+\frac{f^{\prime }(a)}{1!}(x?a)+\frac{f^{\prime \prime }(a)}{2!}(x?a{)}^2+...+\frac{f^{(n)}(a)}{n!}(x?a{)}^n+R_n(x)$$
稱為 f(x) 在 x = a 點關于 x 的冪函數展開式，又稱為 Taylor 公式，式中Rn（x）叫做 Lagrange 余項。

歐拉方法

考慮一階常微分方程的初值問題
$$?\begin{array}{rl}& \frac{dy}{dx}=f(x,y)\\ & y(x_0)=y_0\end{array}$$

前向歐拉法

$$y_{n+1}=y_n+hf(x_n,y_n),n=0,1,...$$

后向歐拉算法

$$y_{n+1}=y_n+hf(x_{n+1},y_{n+1}),n=0,1,...$$
前后向歐拉法的推理、舉例及穩定性對比的超棒分析！地址入口

梯形公式

$$y_{n+1}=y_n+\frac{h}{2}[f(x_n+y_n)+f(x_{n+1}+y_{n+1})],n=0,1,...$$

改進的歐拉算法

$$?\begin{array}{rl}& {\hat{y}}_{n+1}=y_n+hf(x_n,y_n)\\ & y_{n+1}=y_n+\frac{h}{2}[f(x_n,y_n)+f(x_{n+1},{\hat{y}}_{n+1})]\end{array}$$
也可以表示成
$$?\begin{array}{rl}& y_f=y_n+hf(x_n,y_n)\\ & y_b=y_n+hf(x_{n+1},y_f)\\ & y_{n+1}=\frac{1}{2}(x_f+x_b)\end{array}$$

其中$y_f$表示利用向前(顯式)歐拉公式的近似值，$y_b$表示利用向后(隱式)歐拉公式的近似值(利用了$y_f$)，最后取平均值。
上面利用$f(x_{n+1},{\hat{y}}_{n+1})]$是有誤差的，也可通過多次迭代減少誤差，具體公式如下
$$?\begin{array}{rl}& {\hat{y}}_{n+1}^{(0)}=y_n+hf(x_n,y_n)\\ & y_{n+1}^{(k+1)}=y_n+\frac{h}{2}[f(x_n,y_n)+f(x_{n+1},{\hat{y}}_{n+1}^{(k)})],k=0,1,2,...\end{array}$$

其中，后向歐拉算法和梯形公式是隱式算法，前向歐拉算法和改進的歐拉算法是顯式算法。歐拉算法計算容易，但是精度低，梯形公式精度高，但是是隱式形式，不易求解。將兩式結合，則可以得到改進的歐拉公式。先用歐拉公式求出yn+1的一個粗糙的估計值，再用梯形方法進行精確化，稱為校正值。

四階龍格庫塔算法（Runge-kutta method）

$$?\begin{array}{rl}& k_1=f(x_n,y_n)\\ & k_2=f(x_n+\frac{h}{2},y_n+\frac{h}{2}\times k_1)\\ & k_3=f(x_n+\frac{h}{2},y_n+\frac{h}{2}\times k_2)\\ & k_4=f(x_n+h,y_n+h\times k_3)\\ & y_{n+1}=y_n+\frac{(k_1+2k_2+2k_3+k_4)}{6}\times h\end{array}$$
龍格-庫塔算法是一種在工程上應用廣泛的高精度單步算法，算法精度較高，如果預先取四個點就是四階龍格-庫塔算法。
真的！龍格庫塔算法的超強解析！地址入口

四階龍格庫塔函數代碼

runge_kuttx0_o4.m
代碼參考的博客地址入口
總結一下代碼的優點~
1、使用函數句柄，方便復用~
2、模仿了ode45的函數輸入變量，和ode45用起來差不多~

function [t,y,n]=runge_kuttx0_o4(ufunc,tspan,y0,h)%參數表順序依次是微分方程組的函數名稱，時間起點終點,初始值，步長（參數形式參考了ode45函數） if nargin<4h=0.01; end if size(tspan)==[1,2]t0=tspan(1);tn=tspan(2); elseerror(message('MATLAB:runge_kuttx0_o4:WrongDimensionOfTspan')); end n=floor((tn-t0)/h);%求步數 t(1)=t0;%時間起點 y(:,1)=y0;%第一列賦初值，可以是向量，代表不同的初始值 for i=1:n t(i+1)=t(i)+h; k1=ufunc(t(i),y(:,i)); k2=ufunc(t(i)+h/2,y(:,i)+h*k1/2); k3=ufunc(t(i)+h/2,y(:,i)+h*k2/2); k4=ufunc(t(i)+h,y(:,i)+h*k3); y(:,i+1)=y(:,i)+h*(k1+2*k2+2*k3+k4)/6; %按照龍格庫塔方法進行數值求解 end

test1.m

clear clc test_fun=@(t,y)(y+3*t)/t^2; tspan=[1 4]; y0=-2; h=1; [t1,y1]=ode45(test_fun,tspan,y0); [t2,y2]=runge_kuttx0_o4(test_fun,tspan,y0,h); plot(t1,y1,'r',t2,y2,'g') legend('ode45函數效果','自編四階龍格庫塔函數效果') xlabel('t'); ylabel('y'); title('效果對比圖')

test.m運行結果

步長選擇

之前討論的所有的龍格-庫塔方法都是以$\Delta t$定步長來展開的，但從$x_i?x_{i+1}$單步遞推過程來說，步長$\Delta t$越小，局部截斷誤差越小(方法確定情況下)，但是隨著步長的縮小，不但會引起計算量的增加，而且也有可能引起舍入誤差的嚴重積累；但步長$\Delta t$太大又不能達到預期的精度要求，所以選擇合適的步長$\Delta t$，在實際計算中也是比較重要的。其實有時候在實際使用中步長并不需要算法確定，而是需要根據數據幀率來確定的，比如imu數據。??
下面給出求解步長的步驟：

1、以步長$\Delta t$開始，利用龍格-庫塔公式計算$x_i?x_{i+1}$得到一個近似值$x_{i+1}^{\Delta t}$；
2、然后步長減半為$\Delta t/2$，利用龍格-庫塔公式分兩步計算$x_i?x_{i+\frac{1}{2}}?x_{i+1}$得到一個近似值$x_{i+1}^{\Delta t/2}$；
3、計算$\mid x_{i+1}^{\Delta t/2}?x_{i+1}^{\Delta t}<?\mid$是否成立，如果成立直接步長選擇$\Delta t$，否則繼續步長減半，重復上訴步驟直到滿足精度要求。
test2.m

clear clc test_fun=@(t,y)(y+3*t)/t^2; tspan=[1 15]; y0=-2; h=1; acc=0.0001; [t1,y1,n1]=runge_kuttx0_o4(test_fun,tspan,y0,h); tf=t1;yf=y1;hf=h; h=h/2; [t2,y2,n2]=runge_kuttx0_o4(test_fun,tspan,y0,h); while(sum(abs(y1(2:n1+1)-y2(3:2:n2+1)))/n1>=acc) t1=t2;y1=y2;n1=n2; h=h/2; [t2,y2,n2]=runge_kuttx0_o4(test_fun,tspan,y0,h); end plot(tf,yf,'r',t2,y2,'g') legend(['四階龍格庫塔初始步長下的效果,h=',num2str(hf)],['四階龍格庫塔精度約為0.01的效果,h=',num2str(h*2)]) xlabel('t'); ylabel('y'); title('效果對比圖')

總結

以上是生活随笔為你收集整理的四阶龙格库塔算法及matlab代码的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。