Java编辑欧拉定理 欧拉方法求初值问题

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

简介

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

\[\begin{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}\right|_{x=x_{n}} \approx \frac{u\left(x_{n+1}\right)-u\left(x_{n}\right)} {h} \]

代入(1)式，有：

\[u\left(x_{n+1}\right)=y\left(x_{n}\right)+h f\left(x_{n} \| u\left(x_{n}\right)\right) \]

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

\[u_{n+1}=u_{n}+h f\left(x_{n}, u_{n}\right) \]

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

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

类似的，就可以得到：

\[u_{n-1}=u_{n}-h f\left(x_{n}, u_{n}\right) \]

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

改进的欧拉法

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

\[u\left(x_{n+1}\right)=u\left(x_{n}\right)+\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 f\left(x_{n+1}, u\left(x_{n+1}\right)\right) \]

代入上式得：

\[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[f\left(x_{n}, u\left(x_{n}\right)\right)+f\left(x_{n+1}, u\left(x_{n+1}\right)\right)\right] \]

得到：

\[u_{n+1}=u_{n}+\frac{h}{2}\left[f\left(x_{n}, u_{n}\right)+f\left(x_{n+1}, u_{n+1}\right)\right] \]

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

\[\left\{\begin{array}{l} \bar{u}_{n+1}=u_{n}+h f\left(x_{n}, u_{n}\right) \\ u_{n+1}=u_{n}+\frac{h}{2}\left[f\left(x_{n}, u_{n}\right)+f\left(x_{n+1}, \bar{u}_{n+1}\right)\right] \end{array}\right.\]

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

\[\left\{\begin{array}{l} u_{n+1}=u_{n}+h K_{1} \\ K_{1}=f\left(x_{n}, u_{n}\right) \end{array}\right.\]

\[\left\{\begin{array}{l} u_{n+1}=u_{n}+\frac{h}{2}\left(K_{1}+K_{2}\right) \\ K_{1}=f\left(x_{n}, u_{n}\right) \\ K_{2}=f\left(x_{n}+h, u_{n}+K_{1}\right) \end{array}\right.\]

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

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

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

标准四阶显式Kutta公式

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

三级三阶显式公式

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

四级四阶显式Kutta公式

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

四级四阶显式Gill公式

\[\left\{\begin{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}\right) \\ K_{1}=f\left(x_{n}, y_{n}\right) \\ K_{2}=f\left(x_{n}+\frac{1}{2} h, y_{n}+\frac{1}{2} h K_{1}\right) \\ K_{3}=f\left(x_{n}+\frac{1}{2} h, y_{n}+\frac{\sqrt{2}-1}{2} h K_{1}+\left(1-\frac{\sqrt{2}}{2}\right) h K_{2}\right) \\ K_{4}=f\left(x_{n}+h, y_{n}-\frac{\sqrt{2}}{2} h K_{2}+\left(1+\frac{\sqrt{2}}{2}\right) h K_{3}\right)\end{array}\right. \]

三个例题

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

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

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\{\begin{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)=2\end{array}\right.\)
(精确解: \(\left\{\begin{array}{l}x=e^{t} \cos t+2 e^{t} \sin t \\ y=-e^{t} \sin t+2 e^{t} \cos t\end{array}\right)\)

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\{\begin{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)=-3\end{array}\right.\)

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]

运行结果