欧拉法、改进的欧拉法、龙格-库塔法求解初值问题

欧拉法、改进的欧拉法、龙格-库塔法求解初值问题

简介

通过求解简单的初值问题：

[egin{cases}dfrac{du}{dx}=f(x,u)&&&&&&(1)\u(x_0)=u_0&&&&&&(2)end{cases}​ ]

引入欧拉法、改进的欧拉法、龙格-库塔法等。

前期准备

数值解法的基本思想就是先对x和u(x)在区间[x0,∞)上进行离散化，然后构造递推公式，再进一步得到u(x)u(x) u(x)u(x)u(x)u(x)在这些位置的近似取值。

• 取定步长h，令(x_n=x_0+nh(n=±1,±2,⋯))
• 得到离散的位置：(x_1,x_2,⋯,x_n,)
• u(x)在这些点精确取值为：(u(x_1),u(x_2),⋯,u(x_n))
• 利用数值解法得到的这些点的近似取值，(u_1,u_2,cdots,u_n)

欧拉法

欧拉法的核心就是将导数近似为差商。

将导数近似为向前差商，则有：

[left.frac{d u}{d x} ight|_{x=x_{n}} approx frac{uleft(x_{n+1} ight)-uleft(x_{n} ight)} {h} ]

代入(1)式，有：

[uleft(x_{n+1} ight)=yleft(x_{n} ight)+h fleft(x_{n} | uleft(x_{n} ight) ight) ]

用(u_{n+1})​和 (u_n)代替(u(x_{n+1}))和(u(x_n))，得：

[u_{n+1}=u_{n}+h fleft(x_{n}, u_{n} ight) ]

因此，若知道(u_0)我们就可以递归出(u_1,u_2,cdots)
如果将导数近似为向后差商:

[left.frac{d u}{d x} ight|_{x=x_{n}} approx frac{uleft(x_{n} ight)-uleft(x_{n-1} ight)}{h} ]

类似的，就可以得到：

[u_{n-1}=u_{n}-h fleft(x_{n}, u_{n} ight) ]

这样，若知道(u_0)我们就可以递归出(u_{-1}, u_{-2} cdots)

改进的欧拉法

对((1))式在([x_n,x_{n+1}])上积分，可得：

[uleft(x_{n+1} ight)=uleft(x_{n} ight)+int_{x_{n}}^{x_{n+1}} f(x, u) d x ]

其中，(n=0,1,cdots)用不同方式来近似上式的积分运算，就会得到不同的递推公式。若使用左端点计算矩形面积并取近似：

[int_{x_{n}}^{x_{n+1}} f(x, u) d x approx h fleft(x_{n+1}, uleft(x_{n+1} ight) ight) ]

代入上式得：

[u_{n+1}=u_{n}+hf(x_{n},u_{n}) ]

若使用梯形的面积做近似：

[int_{x_{n}}^{x_{n+1}} f(x, y) d x approx frac{h}{2}left[fleft(x_{n}, uleft(x_{n} ight) ight)+fleft(x_{n+1}, uleft(x_{n+1} ight) ight) ight] ]

得到：

[u_{n+1}=u_{n}+frac{h}{2}left[fleft(x_{n}, u_{n} ight)+fleft(x_{n+1}, u_{n+1} ight) ight] ]

欧拉法虽然精度偏低，但它是显式的，可直接得到结果。而梯形公式是隐式的，虽然精度较高，却无法通过一步计算得到结果，若用迭代法计算，运算量较大。综合这两种方法，可以相得益彰：先用显式格式却低精度的欧拉法计算得到一个粗略的预测值(ar{u}_{n+1})，再将这个预测值代入梯形公式进行修正，得到较高精度的结果(u_{n+1})。

[left{egin{array}{l} ar{u}_{n+1}=u_{n}+h fleft(x_{n}, u_{n} ight) \ u_{n+1}=u_{n}+frac{h}{2}left[fleft(x_{n}, u_{n} ight)+fleft(x_{n+1}, ar{u}_{n+1} ight) ight] end{array} ight.]

龙格-库塔法
将以上两种方法分别写成如下形式：

[left{egin{array}{l} u_{n+1}=u_{n}+h K_{1} \ K_{1}=fleft(x_{n}, u_{n} ight) end{array} ight.]

[left{egin{array}{l} u_{n+1}=u_{n}+frac{h}{2}left(K_{1}+K_{2} ight) \ K_{1}=fleft(x_{n}, u_{n} ight) \ K_{2}=fleft(x_{n}+h, u_{n}+K_{1} ight) end{array} ight.]

上述方法都是通过(f(x,u))在不同位置的线性组合来计算(u_{n+1})​的值，所考虑的位置越多，精度也越高。类似的，就得到龙格-库塔法的思想：如果用(f(x,u))在更多位置的线性组合来构造递推公式，将会得到更高的精度。
这样，递推公式将有如下形式：

[left{egin{array}{l} u_{n+1}=u_{n}+h sum_{i=1}^{r} R_{i} K_{i} \ K_{1}=fleft(x_{n}, u_{n} ight) \ K_{i}=fleft(x_{n}+a_{i} h, u_{n}+sum_{j=1}^{i-1} b_{i j} K_{j} ight), i=2,3, cdots, r end{array} ight.]

其中，(R_{i},a_i,b_{ij})为待定常数。(利用(Taylor)展开就可以确定待定系数)

标准四阶显式Kutta公式

[left{egin{array}{l} y_{n+1}=y_{n}+frac{h}{6}left(K_{1}+4 K_{2}+K_{3} ight) \ K_{1}=fleft(x_{n}, y_{n} ight) \ K_{2}=fleft(x_{n}+frac{1}{2} h, y_{n}+frac{1}{2} h K_{1} ight) \ K_{3}=fleft(x_{n}+h, y_{n}-h K_{1}+2 h K_{2} ight) end{array} ight.]

三级三阶显式公式

[left{egin{array}{l} y_{n+1}=y_{n}+frac{h}{4}left(K_{1}+3 K_{3} ight) \ K_{1}=fleft(x_{n}, y_{n} ight) \ K_{2}=fleft(x_{n}+frac{1}{3} h, y_{n}+frac{1}{3} h K_{1} ight) \ K_{3}=fleft(x_{n}+frac{2}{3} h, y_{n}+frac{2}{3} h K_{2} ight) end{array} ight.]

四级四阶显式Kutta公式

[left{egin{array}{l} y_{n+1}=y_{n}+frac{h}{8}left(K_{1}+3 K_{2}+3 K_{3}+K_{4} ight) \ K_{1}=fleft(x_{n}, y_{n} ight) \ K_{2}=fleft(x_{n}+frac{1}{3} h, y_{n}+frac{1}{3} h K_{1} ight) \ K_{3}=fleft(x_{n}+frac{2}{3} h, y_{n}-frac{2}{3} h K_{1}+h K_{2} ight) \ K_{4}=fleft(x_{n}+h, y_{n}+h K_{1}-h K_{2}+h K_{3} ight) end{array} ight.]

四级四阶显式Gill公式

[left{egin{array}{l}y_{n+1}=y_{n}+frac{h}{6}left(K_{1}+(2-sqrt{2}) K_{2}+(2+sqrt{2}) K_{3}+K_{4} ight) \ K_{1}=fleft(x_{n}, y_{n} ight) \ K_{2}=fleft(x_{n}+frac{1}{2} h, y_{n}+frac{1}{2} h K_{1} ight) \ K_{3}=fleft(x_{n}+frac{1}{2} h, y_{n}+frac{sqrt{2}-1}{2} h K_{1}+left(1-frac{sqrt{2}}{2} ight) h K_{2} ight) \ K_{4}=fleft(x_{n}+h, y_{n}-frac{sqrt{2}}{2} h K_{2}+left(1+frac{sqrt{2}}{2} ight) h K_{3} ight)end{array} ight. ]

三个例题

1 分别用 Euler 法，改进的 Euler 法和经典 4 阶龙格-库塔法计算下列初值问题，并绘图比较:

(left{egin{array}{ll}y^{prime}=-y(1+x y) & (0 leq x leq 1) \ y(0)=1 & end{array} quadleft( ext { 精确解 }: quad y(x)=left(2 e^{x}-x-1 ight)^{-1} ight) ight.)

matlab代码

clear all,
close all
f=@(x,y)-y*(1+x*y);
h=0.1;
%% Euler method
x= [0:h:1];
N=size(x,2)-1
y1=[1,zeros(1,N)];
for n=1:N
    y1(n+1)=y1(n)+h*f(x(n),y1(n));
end
%% Improved Euler method
y2=[1,zeros(1,N)];
for n=1:N
    y2(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 formula
y3=[1,zeros(1,N)];
for n=1:N
    K1=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
%% 绘图
y=(2*exp(x)-x-1).^(-1);  %  Exact solution
plot(x,y,'k',x,y1,'xr',x,y2,'ob',x,y3,'*r','Markersize',10,'LineWidth',1.5)
legend('Exact','Euler','Improved Euler','Standard fourth-order Kutta')

运行结果

2.分别用 Euler 法，改进的 Euler 法和经典 4 阶龙格-库塔法计算下列初值问题，并绘图比较:

(2) (left{egin{array}{l}frac{mathrm{d} x}{mathrm{d} t}=x+y, x(0)=1 \ frac{mathrm{d} y}{mathrm{d} t}=-x+y, y(0)=2end{array} ight.)
(精确解: (left{egin{array}{l}x=e^{t} cos t+2 e^{t} sin t \ y=-e^{t} sin t+2 e^{t} cos tend{array} ight))

clear all,
close all
h=0.15; %定义步长
t=0:h:10; %给定参数t的范围
N=size(t,2)-1; 
%% Euler method
Y1(:,1)=[1;2];%赋初值
for n=1:N
    Y1(:,(n+1))=Y1(:,n)+h*F1(t(n),Y1(:,n));
end
x1=Y1(1,:);
y1=Y1(2,:);
%% Improved Euler method
Y2(:,1)=[1;2];
for n=1:N
    Y2(:,(n+1))=Y2(:,n)+h*F1(t(n),Y2(:,n));
    Y2(:,(n+1))=Y2(:,n)+h/2*(F1(t(n),Y2(:,n))+F1(t(n+1),Y2(:,n+1)));
end
x2=Y2(1,:);
y2=Y2(2,:);
%% Standard fourth-order explicit Kutta formula
Y3(:,1)=[1;2];
for n=1:N
    K1=F1(t(n),Y3(:,n));
    K2=F1(t(n)+1/2*h,Y3(:,n)+1/2*h*K1);
    K3=F1(t(n)+h,Y3(:,n)-h*K1+2*h*K2);
    Y3(:,n+1)=Y3(:,n)+h/6*(K1+4*K2+K3);
end
x3=Y3(1,:);
y3=Y3(2,:);
%% 精确解
x=exp(t).*cos(t)+2*exp(t).*sin(t);
y=-exp(t).*sin(t)+2*exp(t).*cos(t);
%% 绘图比较
figure
set(gcf,'position',[0.15 0.2 0.7 0.6])
subplot(1,2,1)
plot(t,x,'k',t,x1,'xr',t,x2,'ob',t,x3,'*r','Markersize',10,'LineWidth',1.5)
legend('Exact','Euler','Improved Euler','Standard fourth-order Kutta')
xlabel('t')
ylabel('x')

subplot(1,2,2)
plot(t,y,'k',t,y1,'xr',t,y2,'ob',t,y3,'*r','Markersize',10,'LineWidth',1.5)
legend('Exact','Euler','Improved Euler','Standard fourth-order Kutta')
xlabel('t')
ylabel('y')

函数脚本

function F1=f(t,Y)
%定义所求微分方程
x=Y(1);
y=Y(2);
f1=x+y;
f2=-x+y;
F1=[f1;f2];
end

运行结果

3.分别用 Euler 法，改进的 Euler 法和经典 4 阶龙格-库塔法计算下列初值问题，并绘图比较:

(3) (left{egin{array}{l}y^{prime prime}=5 e^{2 x} sin x-2 y+2 y^{prime}, quad x in[0,20] \ y(0)=-2, y^{prime}(0)=-3end{array} ight.)

clear all,
close all
h=0.5; %定义步长
x=0:h:20; %给定参数t的范围
N=size(x,2)-1; 
%% Euler method
Y1(:,1)=[-3;-2];%赋初值
for n=1:N
    Y1(:,(n+1))=Y1(:,n)+h*F2(x(n),Y1(:,n));
end
y1=Y1(2,:);
%% Improved Euler method
Y2(:,1)=[-3;-2];
for n=1:N
    Y2(:,(n+1))=Y2(:,n)+h*F2(x(n),Y2(:,n));
    Y2(:,(n+1))=Y2(:,n)+h/2*(F2(x(n),Y2(:,n))+F2(x(n+1),Y2(:,n+1)));
end
y2=Y2(2,:);
%% Standard fourth-order explicit Kutta formula
Y3(:,1)=[-3;-2];
for n=1:N
    K1=F2(x(n),Y3(:,n));
    K2=F2(x(n)+1/2*h,Y3(:,n)+1/2*h*K1);
    K3=F2(x(n)+h,Y3(:,n)-h*K1+2*h*K2);
    Y3(:,n+1)=Y3(:,n)+h/6*(K1+4*K2+K3);
end
y3=Y3(2,:);
%% 绘图比较
plot(x,y1,'xr',x,y2,'ob',x,y3,'*r','Markersize',10,'LineWidth',1.5)
legend('Euler','Improved Euler','Standard fourth-order Kutta')
xlabel('x')
ylabel('y')

函数脚本

function F2=f(x,Y)
%定义所求微分方程
z=Y(1);
y=Y(2);
f1=5*exp(2*x).*sin(x)-2*y+2*z;
f2=z;
F2=[f1;f2]

运行结果