歐拉法、改進的歐拉法、龍格-庫塔法求解初值問題

簡介

通過求解簡單的初值問題：

\[\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]

運行結果