这# 欧拉法、改进的欧拉法、龙格-库塔法求解初值问题
简介
通过求解简单的初值问题：
{
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和u(x)在区间[x0,∞)上进行离散化，然后构造递推公式，再进一步得到u(x)u(x) u(x)u(x)u(x)u(x)在这些位置的近似取值。
- 取定步长h，令 x n = x 0 + n h ( n = ± 1 , ± 2 , ⋯ ) x_n=x_0+nh(n=±1,±2,⋯) xn=x0+nh(n=±1,±2,⋯)
- 得到离散的位置： x 1 , x 2 , ⋯ , x n , x_1,x_2,⋯,x_n, x1,x2,⋯,xn,
- u(x)在这些点精确取值为： u ( x 1 ) , u ( x 2 ) , ⋯ , u ( x n ) u(x_1),u(x_2),⋯,u(x_n) u(x1),u(x2),⋯,u(xn)
- 利用数值解法得到的这些点的近似取值， u 1 , u 2 , ⋯ , u n u_1,u_2,\cdots,u_n u1,u2,⋯,un
欧拉法
欧拉法的核心就是将导数近似为差商。
将导数近似为向前差商，则有：
d
u
d
x
∣
x
=
x
n
≈
u
(
x
n
+
1
)
−
u
(
x
n
)
h
\left.\frac{d u}{d x}\right|_{x=x_{n}} \approx \frac{u\left(x_{n+1}\right)-u\left(x_{n}\right)} {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\left(x_{n+1}\right)=y\left(x_{n}\right)+h f\left(x_{n} \| u\left(x_{n}\right)\right)
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}+h f\left(x_{n}, u_{n}\right)
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.\frac{d u}{d x}\right|_{x=x_{n}} \approx \frac{u\left(x_{n}\right)-u\left(x_{n-1}\right)}{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}-h f\left(x_{n}, u_{n}\right)
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\left(x_{n+1}\right)=u\left(x_{n}\right)+\int_{x_{n}}^{x_{n+1}} f(x, u) d x
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}}^{x_{n+1}} f(x, u) d x \approx h f\left(x_{n+1}, u\left(x_{n+1}\right)\right)
∫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}}^{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]
∫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}\left[f\left(x_{n}, u_{n}\right)+f\left(x_{n+1}, u_{n+1}\right)\right]
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
)
]
\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.
{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
)
\left\{\begin{array}{l} u_{n+1}=u_{n}+h K_{1} \\ K_{1}=f\left(x_{n}, u_{n}\right) \end{array}\right.
{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
)
\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.
⎩⎨⎧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
\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.
⎩⎪⎨⎪⎧un+1=un+h∑i=1rRiKiK1=f(xn,un)Ki=f(xn+aih,un+∑j=1i−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 ) \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. ⎩⎪⎪⎨⎪⎪⎧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 ) \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. ⎩⎪⎪⎨⎪⎪⎧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 ) \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. ⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧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 ) \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. ⎩⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎧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−22hK2+(1+22)hK3)
三个例题
1 分别用 Euler 法，改进的 Euler 法和经典 4 阶龙格-库塔法计算下列初值问题，并绘图比较:
{ y ′ = − y ( 1 + x y ) ( 0 ≤ x ≤ 1 ) y ( 0 ) = 1 ( 精确解 : y ( x ) = ( 2 e x − x − 1 ) − 1 ) \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. {y′=−y(1+xy)y(0)=1(0≤x≤1)( 精确解 :y(x)=(2ex−x−1)−1)
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)
{
d
x
d
t
=
x
+
y
,
x
(
0
)
=
1
d
y
d
t
=
−
x
+
y
,
y
(
0
)
=
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.
{dtdx=x+y,x(0)=1dtdy=−x+y,y(0)=2
(精确解:
{
x
=
e
t
cos
t
+
2
e
t
sin
t
y
=
−
e
t
sin
t
+
2
e
t
cos
t
)
\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)
{x=etcost+2etsinty=−etsint+2etcost)
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) { y ′ ′ = 5 e 2 x sin x − 2 y + 2 y ′ , x ∈ [ 0 , 20 ] y ( 0 ) = − 2 , y ′ ( 0 ) = − 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. {y′′=5e2xsinx−2y+2y′,x∈[0,20]y(0)=−2,y′(0)=−3
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]
3834
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
