求解初值问题
简介前期准备欧拉法改进的欧拉法龙格-库塔法标准四阶显式Kutta公式三级三阶显式公式四级四阶显式Kutta公式四级四阶显式Gill公式 示例MATLAB代码结果简介
通过求解简单的初值问题：
{ d u d x = f ( x , u ) ( 1 ) u ( x 0 ) = u 0 ( 2 ) \begin{cases} \dfrac{du}{dx}=f(x,u)&&&&&&(1)\\ u(x_0)=u_0&&&&&&(2) \end{cases} ⎩⎨⎧dxdu=f(x,u)u(x0)=u0(1)(2)
引入欧拉法、改进的欧拉法、龙格-库塔法等。
前期准备
数值解法的基本思想就是先对 x x x和 u ( x ) u(x) u(x)在区间 [ x 0 , ∞ ) [x_0,\infty) [x0,∞)上进行离散化，然后构造递推公式，再进一步得到 u ( x ) u(x) u(x)在这些位置的近似取值。
取定步长 h h h，令 x n = x 0 + n h ( n = ± 1 , ± 2 , ⋯ ) x_n=x_0+nh(n=\pm1,\pm2,\cdots) xn=x0+nh(n=±1,±2,⋯)得到离散的位置： x 1 , x 2 , ⋯ , x n , ⋯ x_1,x_2,\cdots,x_n,\cdots x1,x2,⋯,xn,⋯ u ( x ) u(x) u(x)在这些点精确取值为： u ( x 1 ) , u ( x 2 ) , ⋯ , u ( x n ) , ⋯ u(x_1),u(x_2),\cdots,u(x_n),\cdots u(x1),u(x2),⋯,u(xn),⋯利用数值解法得到的这些点的近似取值，记为 u 1 , u 2 , ⋯ , u n , ⋯ u_1,u_2,\cdots,u_n,\cdots u1,u2,⋯,un,⋯
欧拉法
欧拉法的核心就是将导数近似为差商。
将导数近似为向前差商，则有：
d u d x ∣ x = x n ≈ u ( x n + 1 ) − u ( x n ) h \left.\dfrac{du}{dx}\right|_{x=x_n} \approx \frac{u(x_{n+1})-u(x_n)}{h} dxdu∣∣∣∣x=xn≈hu(xn+1)−u(xn)
代入(1)式，有：
u ( x n + 1 ) = y ( x n ) + h f ( x n , u ( x n ) ) u(x_{n+1})=y(x_{n})+hf(x_n,u(x_n)) u(xn+1)=y(xn)+hf(xn,u(xn))
用 u n + 1 u_{n+1} un+1和 u n u_n un代替 u ( x n + 1 ) u(x_{n+1}) u(xn+1)和 u ( x n ) u(x_n) u(xn)，得：
u n + 1 = u n + h f ( x n , u n ) u_{n+1}=u_{n}+hf(x_n,u_n) un+1=un+hf(xn,un)
因此，若知道 u 0 u_0 u0我们就可以递归出 u 1 , u 2 , ⋯ u_1,u_2,\cdots u1,u2,⋯
如果将导数近似为向后差商:
d u d x ∣ x = x n ≈ u ( x n ) − u ( x n − 1 ) h \left.\dfrac{du}{dx}\right|_{x=x_n} \approx \frac{u(x_{n})-u(x_{n-1})}{h} dxdu∣∣∣∣x=xn≈hu(xn)−u(xn−1)
类似的，就可以得到：
u n − 1 = u n − h f ( x n , u n ) u_{n-1}=u_{n}-hf(x_n,u_n) un−1=un−hf(xn,un)
这样，若知道 u 0 u_0 u0我们就可以递归出 u − 1 , u − 2 , ⋯ u_{-1},u_{-2},\cdots u−1,u−2,⋯
改进的欧拉法
对 ( 1 ) (1) (1)式在 [ x n , x n + 1 ] [x_n,x_{n+1}] [xn,xn+1]上积分，可得：
u ( x n + 1 ) = u ( x n ) + ∫ x n x n + 1 f ( x , u ) d x u(x_{n+1})=u(x_{n})+\int^{x_{n+1}}_{x_n}f(x,u)dx u(xn+1)=u(xn)+∫xnxn+1f(x,u)dx
其中， n = 0 , 1 , ⋯ n=0,1,\cdots n=0,1,⋯。用不同方式来近似上式的积分运算，就会得到不同的递推公式。若使用左端点计算矩形面积并取近似：
∫ x n x n + 1 f ( x , u ) d x ≈ h f ( x n + 1 , u ( x n + 1 ) ) \int^{x_{n+1}}_{x_n}f(x,u)dx\approx hf(x_{n+1},u(x_{n+1})) ∫xnxn+1f(x,u)dx≈hf(xn+1,u(xn+1))
代入上式得：
u n + 1 = u n + h f ( x n , u n ) u_{n+1}=u_{n}+hf(x_{n},u_{n}) un+1=un+hf(xn,un)
若使用梯形的面积做近似：
∫ x n x n + 1 f ( x , y ) d x ≈ h 2 [ f ( x n , u ( x n ) ) + f ( x n + 1 , u ( x n + 1 ) ) ] \int^{x_{n+1}}_{x_n}f(x,y)dx\approx\frac{h}{2}[f(x_{n},u(x_{n}))+f(x_{n+1},u(x_{n+1}))] ∫xnxn+1f(x,y)dx≈2h[f(xn,u(xn))+f(xn+1,u(xn+1))]
得到：
u n + 1 = u n + h 2 [ f ( x n , u n ) + f ( x n + 1 , u n + 1 ) ] u_{n+1}=u_{n}+\frac{h}{2}[f(x_{n},u_{n})+f(x_{n+1},u_{n+1})] un+1=un+2h[f(xn,un)+f(xn+1,un+1)]
欧拉法虽然精度偏低，但它是显式的，可直接得到结果。而梯形公式是隐式的，虽然精度较高，却无法通过一步计算得到结果，若用迭代法计算，运算量较大。综合这两种方法，可以相得益彰：先用显式格式却低精度的欧拉法计算得到一个粗略的预测值 u ˉ n + 1 \bar{u}_{n+1} uˉn+1，再将这个预测值代入梯形公式进行修正，得到较高精度的结果 u n + 1 u_{n+1} un+1。
{ u ˉ n + 1 = u n + h f ( x n , u n ) u n + 1 = u n + h 2 [ f ( x n , u n ) + f ( x n + 1 , u ˉ n + 1 ) ] \begin{cases} \bar{u}_{n+1}=u_{n}+hf(x_{n},u_{n})\\ u_{n+1}=u_n+\dfrac{h}{2}[f(x_{n},u_{n})+f(x_{n+1},\bar{u}_{n+1})] \end{cases} ⎩⎨⎧uˉn+1=un+hf(xn,un)un+1=un+2h[f(xn,un)+f(xn+1,uˉn+1)]
龙格-库塔法
将以上两种方法分别写成如下形式：
{ u n + 1 = u n + h K 1 K 1 = f ( x n , u n ) \begin{cases} u_{n+1}=u_{n}+hK_1\\ K_1=f(x_{n},u_{n}) \end{cases} {un+1=un+hK1K1=f(xn,un)
{ u n + 1 = u n + h 2 ( K 1 + K 2 ) K 1 = f ( x n , u n ) K 2 = f ( x n + h , u n + K 1 ) \begin{cases} u_{n+1}=u_{n}+\dfrac{h}{2}(K_1+K_2)\\ K_1=f(x_{n},u_{n})\\ K_2=f(x_{n}+h,u_{n}+K_1) \end{cases} ⎩⎪⎪⎨⎪⎪⎧un+1=un+2h(K1+K2)K1=f(xn,un)K2=f(xn+h,un+K1)
上述方法都是通过 f ( x , u ) f(x,u) f(x,u)在不同位置的线性组合来计算 u n + 1 u_{n+1} un+1的值，所考虑的位置越多，精度也越高。类似的，就得到龙格-库塔法的思想：如果用 f ( x , u ) f(x,u) f(x,u)在更多位置的线性组合来构造递推公式，将会得到更高的精度。
这样，递推公式将有如下形式：
{ u n + 1 = u n + h ∑ i = 1 r R i K i K 1 = f ( x n , u n ) K i = f ( x n + a i h , u n + ∑ j = 1 i − 1 b i j K j ) , i = 2 , 3 , ⋯ , r \begin{cases} u_{n+1}=u_{n}+h\sum\limits_{i=1}^{r} R_i K_i\\ K_1=f(x_{n},u_{n})\\ K_i=f(x_{n}+a_i h,u_{n}+\sum\limits_{j=1}^{i-1} b_{ij} K_j), i=2,3,\cdots,r \end{cases} ⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧un+1=un+hi=1∑rRiKiK1=f(xn,un)Ki=f(xn+aih,un+j=1∑i−1bijKj),i=2,3,⋯,r
其中， R i , a i , b i j R_{i},a_i,b_{ij} Ri,ai,bij为待定常数。(利用 T a y l o r Taylor Taylor展开就可以确定待定系数)
标准四阶显式Kutta公式
{ y n + 1 = y n + h 6 ( K 1 + 4 K 2 + K 3 ) , K 1 = f ( x n , y n ) , K 2 = f ( x n + 1 2 h , y n + 1 2 h K 1 ) , K 3 = f ( x n + h , y n − h K 1 + 2 h K 2 ) ; \begin{cases} y_{n+1}=y_n+\dfrac{h}{6}(K_1+4K_2+K_3),\\ K_1=f(x_n,y_n),\\ K_2=f(x_n+\frac{1}{2}h,y_n+\frac{1}{2}h K_1),\\ K_3=f(x_n+h,y_n-h K_1+2h K_2); \end{cases} ⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧yn+1=yn+6h(K1+4K2+K3),K1=f(xn,yn),K2=f(xn+21h,yn+21hK1),K3=f(xn+h,yn−hK1+2hK2);
三级三阶显式公式
{ y n + 1 = y n + h 4 ( K 1 + 3 K 3 ) , K 1 = f ( x n , y n ) , K 2 = f ( x n + 1 3 h , y n + 1 3 h K 1 ) , K 3 = f ( x n + 2 3 h , y n + 2 3 h K 2 ) ; \begin{cases} y_{n+1}=y_n+\dfrac{h}{4}(K_1+3K_3),\\ K_1=f(x_n,y_n),\\ K_2=f(x_n+\frac{1}{3}h,y_n+\frac{1}{3}h K_1),\\ K_3=f(x_n+\frac{2}{3}h,y_n+\frac{2}{3}h K_2); \end{cases} ⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧yn+1=yn+4h(K1+3K3),K1=f(xn,yn),K2=f(xn+31h,yn+31hK1),K3=f(xn+32h,yn+32hK2);
四级四阶显式Kutta公式
{ y n + 1 = y n + h 8 ( K 1 + 3 K 2 + 3 K 3 + K 4 ) , K 1 = f ( x n , y n ) , K 2 = f ( x n + 1 3 h , y n + 1 3 h K 1 ) , K 3 = f ( x n + 2 3 h , y n − 2 3 h K 1 + h K 2 ) , K 4 = f ( x n + h , y n + h K 1 − h K 2 + h K 3 ) ; \begin{cases} y_{n+1}=y_n+\frac{h}{8}(K_1+3K_2+3K_3+K_4),\\ K_1=f(x_n,y_n),\\ K_2=f(x_n+\frac{1}{3}h,y_n+\frac{1}{3}h K_1),\\ K_3=f(x_n+\frac{2}{3}h,y_n-\frac{2}{3}h K_1+h K_2),\\ K_4=f(x_n+h,y_n+h K_1-h K_2+h K_3); \end{cases} ⎩⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎧yn+1=yn+8h(K1+3K2+3K3+K4),K1=f(xn,yn),K2=f(xn+31h,yn+31hK1),K3=f(xn+32h,yn−32hK1+hK2),K4=f(xn+h,yn+hK1−hK2+hK3);
四级四阶显式Gill公式
{ y n + 1 = y n + h 6 ( K 1 + ( 2 − 2 ) K 2 + ( 2 + 2 ) K 3 + K 4 ) , K 1 = f ( x n , y n ) , K 2 = f ( x n + 1 2 h , y n + 1 2 h K 1 ) , K 3 = f ( x n + 1 2 h , y n + 2 − 1 2 h K 1 + ( 1 − 2 2 ) h K 2 ) , K 4 = f ( x n + h , y n − 2 2 h K 2 + ( 1 + 2 2 ) h K 3 ) ; \begin{cases} y_{n+1}=y_n+\frac{h}{6}(K_1+(2-\sqrt{2})K_2+(2+\sqrt{2})K_3+K_4),\\ K_1=f(x_n,y_n),\\ K_2=f(x_n+\frac{1}{2}h,y_n+\frac{1}{2}h K_1),\\ K_3=f(x_n+\frac{1}{2}h,y_n+\frac{\sqrt{2}-1}{2}h K_1+(1-\frac{\sqrt{2}}{2})h K_2),\\ K_4=f(x_n+h,y_n-\frac{\sqrt{2}}{2}h K_2+(1+\frac{\sqrt{2}}{2})h K_3); \end{cases} ⎩⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎧yn+1=yn+6h(K1+(2−2 )K2+(2+2 )K3+K4),K1=f(xn,yn),K2=f(xn+21h,yn+21hK1),K3=f(xn+21h,yn+22 −1hK1+(1−22 )hK2),K4=f(xn+h,yn−22 hK2+(1+22 )hK3);
示例
求解常微分方程:
{ d y d x = x 3 − y x , y ( 1 ) = 2 5 . \begin{cases} \dfrac{dy}{dx}=x^3-\dfrac{y}{x},\\ y(1)=\dfrac{2}{5}. \end{cases} ⎩⎪⎨⎪⎧dxdy=x3−xy,y(1)=52.
要求步长为 h = 0.1 h=0.1 h=0.1，其中，精确解为
y = 1 5 x 4 + 1 5 x . y=\frac{1}{5}x^4+\frac{1}{5x}. y=51x4+5x1.
MATLAB代码
clear all,close all,clcf=@(x,y)x^3-y/x;h=0.1;%% Euler methodx=[1:h:2];N=size(x,2)-1;y1=[2/5,zeros(1,N)];for n=1:Ny1(n+1)=y1(n)+h*f(x(n),y1(n));end%% Improved Euler methody2=[2/5,zeros(1,N)];for n=1:Ny2(n+1)=y2(n)+h*f(x(n),y2(n));y2(n+1)=y2(n)+h/2*(f(x(n),y2(n))+f(x(n+1),y2(n+1)));end%% Standard fourth-order explicit Kutta formulay3=[2/5,zeros(1,N)];for n=1:NK1=f(x(n),y3(n));K2=f(x(n)+1/2*h,y3(n)+1/2*h*K1);K3=f(x(n)+h,y3(n)-h*K1+2*h*K2);y3(n+1)=y3(n)+h/6*(K1+4*K2+K3);end%% Three-level three-order explicit formulay4=[2/5,zeros(1,N)];for n=1:NK1=f(x(n),y4(n));K2=f(x(n)+1/3*h,y4(n)+1/3*h*K1);K3=f(x(n)+2/3*h,y4(n)+2/3*h*K2);y4(n+1)=y4(n)+h/4*(K1+3*K3);end%% Fourth-level fourth-order explicit Kutta formulay5=[2/5,zeros(1,N)];for n=1:NK1=f(x(n),y5(n));K2=f(x(n)+1/3*h,y5(n)+1/3*h*K1);K3=f(x(n)+2/3*h,y5(n)-1/3*h*K1+h*K2);K4=f(x(n)+h,y5(n)+h*K1-h*K2+h*K3);y5(n+1)=y5(n)+h/8*(K1+3*K2+3*K3+K4);end%% Fourth-level fourth-order explicit Gill formulay6=[2/5,zeros(1,N)];for n=1:NK1=f(x(n),y6(n));K2=f(x(n)+1/2*h,y6(n)+1/2*h*K1);K3=f(x(n)+1/2*h,y6(n)+(sqrt(2)/2-0.5)*h*K1+(1-sqrt(2)/2)*h*K2);K4=f(x(n)+h,y6(n)-sqrt(2)/2*h*K2+(1+sqrt(2)/2)*h*K3);y6(n+1)=y6(n)+h/6*(K1+(2-sqrt(2))*K2+(2+sqrt(2))*K3+K4);endy=1/5*x.^4+1./(5*x); % Exact solutionplot(x,y,'k',x,y1,'xr',x,y2,'ob','Markersize',10,'LineWidth',1.5)axis([1 2.1 0 4.5]),xlabel x,ylabel ulegend('Exact','Euler','Improved Euler')%plot(x,y,'k',x,y3,'*r',x,y4,'+b',x,y5,'-y',x,y6,'or','Markersize',10,'LineWidth',1.5)%axis([1 2.1 0 6]),xlabel x,ylabel u,set(gca,'Fontsize',18)%legend('Exact','34Kutta','33kutta','44Kutta','44Gill')
