常微分方程初边值问题的数值解法
2023年11月30日
#analysis
文章目录
存在惟一解
一阶常微分方程初值问题的一般形式为:
{
d
y
d
x
=
f
(
x
,
y
)
,
a
≤
x
≤
b
y
(
a
)
=
α
\begin{cases} \frac{\mathrm d y}{\mathrm dx} =f(x,y) ,&a\le x\le b \\ \\ y(a)= \alpha \end{cases}
⎩
⎨
⎧dxdy=f(x,y),y(a)=αa≤x≤b
其中
f
{f}
f 是
x
{x}
x 和
y
{y}
y 的已知函数,
α
{ \alpha }
α 为给定的初值。
Lipschitz利普希兹条件
若函数
f
(
x
,
y
)
{f(x,y)}
f(x,y) 在区域
{
a
≤
x
≤
b
,
m
<
y
<
M
}
{ \lbrace a\le x\le b,m<y<M \rbrace}
{a≤x≤b,m<y<M} 上连续,且关于
y
{y}
y 满足Lipschitz条件:
∣
f
(
x
,
y
)
−
f
(
x
,
y
ˉ
)
∣
≤
L
∣
y
−
y
ˉ
∣
,
∀
y
,
y
ˉ
|f(x,y)-f(x,\bar y)|\le L|y-\bar y| \,\,,\,\, \forall y,\bar y
∣f(x,y)−f(x,yˉ)∣≤L∣y−yˉ∣,∀y,yˉ
其中
L
>
0
{L>0}
L>0 为Lipschiitz常数,则初值问题(在
a
{a }
a 的某邻域内)存在唯一解。此时Lipschitz常数不必小于
1
{1}
1 。
[!example]-
d y d x = 1 + x sin ( x y ) , x ∈ [ 0 , 2 ] \frac{\mathrm d y}{\mathrm dx}=1+x\sin(xy) \,\,,\,\, x\in[0,2] dxdy=1+xsin(xy),x∈[0,2]
y ( 0 ) = 1 y(0)=1 y(0)=1
解:对任意 y {y} y , y ˉ {\bar y} yˉ ,对变量 y {y} y 应用微分中值定理,存在 y {y} y ,使得
f ( x , y ) − f ( x , y ˉ ) y − y ˉ = ∂ ∂ y f ( x , y ) = x 2 cos ( x y ) \frac{f(x,y)-f(x,\bar y)}{y-\bar y}= \frac{\partial }{\partial y}f(x,y)=x^2\cos(xy) y−yˉf(x,y)−f(x,yˉ)=∂y∂f(x,y)=x2cos(xy)
于是有
∣ f ( x , y ) − f ( x , y ˉ ) ∣ = ∣ x 2 cos ( x y ) ∣ ∣ y − y ˉ ∣ ≤ 4 ∣ y − y ˉ ∣ |f(x,y)-f(x,\bar y)|=|x^2\cos(xy)||y-\bar y|\le 4|y-\bar y| ∣f(x,y)−f(x,yˉ)∣=∣x2cos(xy)∣∣y−yˉ∣≤4∣y−yˉ∣
因而此时 f ( x , y ) {f(x,y)} f(x,y) 关于 y {y} y 满足Lipschitz条件,且常数 L = 4 {L=4} L=4 。因此从理论上讲,初值问题存在惟一解
差分公式的格式
Euler公式
y
n
+
1
=
y
n
+
h
f
(
x
n
,
y
n
)
y_{n+1}=y_n+hf(x_n,y_n)
yn+1=yn+hf(xn,yn)
matlab实现
%% 欧拉公式例子 例9.1.1
[x,y] = odeEuler(@f,1,0,1,10);
[x' y']
function r = f(x,y)
r = y-2*x./y; % 导数公式
end
%% 欧拉公式求常微分方程初边值问题
% 输入函数,初值,起点,终点,区间数
% 输出每个点与对应的函数值
function [x,y] = odeEuler(f,y0,a,b,n)
h = (b-a)/n;
x = linspace(a,b,n+1);
y(1)=y0;
for i = 1:n
y(i+1) = y(i)+h*f(x(i),y(i));
end
end
梯形公式
y
n
+
1
=
y
n
+
h
2
(
f
(
x
n
,
y
n
)
+
f
(
x
n
+
1
,
y
n
+
1
)
)
y_{n+1}=y_n+ \frac{h}{2}(f(x_n,y_n)+f(x_{n+1},y_{n+1}))
yn+1=yn+2h(f(xn,yn)+f(xn+1,yn+1))
matlab实现,由于方程右边用到了左边的值,这里通过不动点迭代出
y
n
+
1
{y_{n+1}}
yn+1
%% 例9.1.2
[x,y] = odeTrapz(@f,1,0,1,10);
[x' y']
function r = f(x,y)
r = y-2*x./y;
end
%% 梯形公式求常微分方程初边值问题
% 输入函数,初值,起点,终点,区间数
% 输出每个点与对应的函数值
function [x,y] = odeTrapz(f,y0,a,b,n)
h = (b-a)/n;
x = linspace(a,b,n+1);
y(1)=y0;
for i = 1:n
t(1) = y(i)
for j = 1:10000 % 不动点迭代法
t(j+1) = y(i)+h/2*f(x(i),y(i))+h/2*f(x(i+1),t(j));
if abs(t(j+1)-t(j))<1e-12
y(i+1) = t(j+1);
break;
end
end
end
end
Euler中点公式
y
n
+
1
=
y
n
−
1
+
2
h
f
(
x
n
,
y
n
)
y_{n+1}=y_{n-1}+2hf(x_n,y_n)
yn+1=yn−1+2hf(xn,yn)
课本上没提,这里不写了,接手这篇笔记的人可以补上。
改进Euler方法(预估-矫正公式)
{
y
ˉ
n
+
1
=
y
n
+
h
f
(
x
n
,
y
n
)
预估
y
n
+
1
=
y
n
+
h
2
(
f
(
x
n
,
y
n
)
+
f
(
x
n
+
1
,
y
ˉ
n
+
1
)
)
校正
y
0
=
α
,
n
=
0
,
1
,
2
,
⋯
,
N
−
1
\begin{cases} \bar y_{n+1}=y_n+hf(x_n,y_n) & 预估\\ y_{n+1}=y_n+ \frac{h}{2}(f(x_n,y_n)+f(x_{n+1},\bar y_{n+1})) & 校正\\ y_0= \alpha ,n=0,1,2, \cdots ,N-1 \end{cases}
⎩
⎨
⎧yˉn+1=yn+hf(xn,yn)yn+1=yn+2h(f(xn,yn)+f(xn+1,yˉn+1))y0=α,n=0,1,2,⋯,N−1预估校正
等价的平均化形式比上面的少算一次
f
(
x
,
y
)
{f(x,y)}
f(x,y)
{
y
p
=
y
n
+
h
f
(
x
n
,
y
n
)
y
c
=
y
n
+
h
f
(
x
n
+
1
,
y
p
)
y
n
+
1
=
0.5
(
y
p
+
y
c
)
\begin{cases} y_p=y_n+hf(x_n,y_n) \\ y_c = y_n+hf(x_{n+1},y_p) \\ y_{n+1}= 0.5(y_p+y_c) \end{cases}
⎩
⎨
⎧yp=yn+hf(xn,yn)yc=yn+hf(xn+1,yp)yn+1=0.5(yp+yc)
等价二阶二段龙格库塔格式
{
y
n
+
1
=
y
n
+
h
2
[
k
1
+
k
2
]
k
1
=
f
(
x
n
,
y
n
)
k
2
=
f
(
x
n
+
h
,
y
n
+
h
k
1
)
y
0
=
α
,
n
=
0
,
1
,
2
,
⋯
,
N
−
1
\begin{cases} y_{n+1}=y_n+ \frac{h}{2}[k_1+k_2]\\ k_1=f(x_n,y_n)\\ k_2=f(x_n+h,y_n+hk_1)\\ y_0= \alpha ,n=0,1,2, \cdots ,N-1 \end{cases}
⎩
⎨
⎧yn+1=yn+2h[k1+k2]k1=f(xn,yn)k2=f(xn+h,yn+hk1)y0=α,n=0,1,2,⋯,N−1
改进欧拉公式精度高于欧拉公式,低于梯形公式,但比起梯形公式是显式格式,计算量较小,相对折中。
matlab实现
[x,y] = odeIEuler(@f,1,0,1,10);
[x' y']
function r = f(x,y)
r = y-2*x./y;
end
%% 改进Euler方法 平均化形式实现
% 输入函数,初值,起点,终点,区间数
% 输出每个点与对应的函数值
function [x,y] = odeIEuler(f,y0,a,b,n)
h = (b-a)/n;
x = linspace(a,b,n+1);
y(1) = y0;
for i = 1:n
yp = y(i)+h*f(x(i),y(i));
yc = y(i)+h*f(x(i+1),yp);
y(i+1) = 0.5*(yp+yc);
end
end
局部截断误差 y ( x n + 1 ) − y n + 1 {y(x_{n+1})-y_{n+1}} y(xn+1)−yn+1
局部截断误差的阶数
若单步差分公式的局部截断误差为
O
(
h
P
+
1
)
{O(h^{P+1})}
O(hP+1) ,则称该公式为
P
{P}
P 阶方法。
Taylor公式
一元Taylor公式:
y
(
x
n
+
1
)
=
y
(
x
n
+
h
)
=
y
(
x
n
)
+
y
′
(
x
n
)
h
+
y
′
′
(
x
n
)
2
!
h
2
+
y
′
′
′
(
x
n
)
3
!
h
3
+
O
(
h
4
)
y(x_{n+1})=y(x_n+h)=y(x_n)+y'(x_n)h+ \frac{y''(x_n)}{2!}h^2+ \frac{y'''(x_n)}{3!}h^3+O(h^4)
y(xn+1)=y(xn+h)=y(xn)+y′(xn)h+2!y′′(xn)h2+3!y′′′(xn)h3+O(h4)
二元Taylor公式:
f
(
x
n
+
h
,
y
n
+
k
)
=
f
(
x
n
,
y
n
)
+
∂
f
(
x
n
,
y
n
)
∂
x
h
+
∂
f
(
x
n
,
y
n
)
∂
y
k
+
1
2
!
(
∂
2
f
(
x
n
,
y
n
)
∂
x
2
h
2
+
2
∂
2
f
(
x
n
,
y
n
)
∂
x
∂
y
h
k
+
∂
2
f
(
x
n
,
y
n
)
∂
y
2
k
2
)
+
1
m
!
(
h
∂
∂
x
+
k
∂
∂
y
)
m
f
(
x
n
,
y
n
)
+
⋯
\begin{align*} f(x_n+h,y_n+k)=&f(x_n,y_n)+ \frac{\partial f(x_n,y_n)}{\partial x}h+ \frac{\partial f(x_n,y_n)}{\partial y}k \\ \\ &+ \frac{1}{2!} \bigg( \frac{\partial ^2f(x_n,y_n)}{\partial x^2} h^2+2 \frac{\partial ^2f(x_n,y_n)}{\partial x\partial y}hk + \frac{\partial ^2f(x_n,y_n)}{\partial y^2}k^2 \bigg) \\ \\ &+ \frac{1}{m!} \bigg( h \frac{\partial }{\partial x}+k \frac{\partial }{\partial y} \bigg)^mf(x_n,y_n)+ \cdots \end{align*}
f(xn+h,yn+k)=f(xn,yn)+∂x∂f(xn,yn)h+∂y∂f(xn,yn)k+2!1(∂x2∂2f(xn,yn)h2+2∂x∂y∂2f(xn,yn)hk+∂y2∂2f(xn,yn)k2)+m!1(h∂x∂+k∂y∂)mf(xn,yn)+⋯
y
′
(
x
n
)
,
y
′
′
(
x
n
)
,
y
′
′
′
(
x
n
)
{y'(x_n),y''(x_n),y'''(x_n)}
y′(xn),y′′(xn),y′′′(xn) 的写法
因
y
′
(
x
)
=
f
(
x
,
y
(
x
)
)
{y'(x)=f(x,y(x))}
y′(x)=f(x,y(x)) ,所以若
y
(
x
n
)
=
y
n
{y(x_n)=y_n}
y(xn)=yn ,则有
y
′
(
x
n
)
=
f
(
x
n
,
y
(
x
n
)
)
=
f
(
x
n
,
y
n
)
=
f
n
y
′
′
(
x
n
)
=
∂
f
(
x
n
,
y
n
)
∂
x
+
∂
f
(
x
n
,
y
n
)
∂
y
y
′
(
x
n
)
=
∂
f
n
∂
x
+
∂
f
n
∂
y
f
n
y
′
′
′
(
x
n
)
=
∂
2
f
n
∂
x
2
+
2
∂
2
f
n
∂
x
∂
y
f
n
+
∂
2
f
n
∂
y
2
f
n
2
+
∂
f
n
∂
x
∂
f
n
∂
y
+
(
∂
f
n
∂
y
)
2
f
n
\begin{align*} y'(x_n)=&f(x_n,y(x_n))=f(x_n,y_n)=f_n \\ \\ y''(x_n)=& \frac{\partial f(x_n,y_n)}{\partial x}+ \frac{\partial f(x_n,y_n)}{\partial y}y'(x_n)= \frac{\partial f_n}{\partial x}+ \frac{\partial f_n}{\partial y}f_n \\ \\ y'''(x_n)=& \frac{\partial ^2f_n}{\partial x^2}+ 2 \frac{\partial ^2f_n}{\partial x\partial y}f_n+ \frac{\partial ^2f_n}{\partial y^2}f_n^2+ \frac{\partial f_n}{\partial x} \frac{\partial f_n}{\partial y}+ (\frac{\partial f_n}{\partial y})^2f_n \end{align*}
y′(xn)=y′′(xn)=y′′′(xn)=f(xn,y(xn))=f(xn,yn)=fn∂x∂f(xn,yn)+∂y∂f(xn,yn)y′(xn)=∂x∂fn+∂y∂fnfn∂x2∂2fn+2∂x∂y∂2fnfn+∂y2∂2fnfn2+∂x∂fn∂y∂fn+(∂y∂fn)2fn
关于改进Euler方法的局部截断误差分析
改进Euler方法格式:
{
y
n
+
1
=
y
n
+
h
2
[
k
1
+
k
2
]
k
1
=
f
(
x
n
,
y
n
)
k
2
=
f
(
x
n
+
h
,
y
n
+
h
k
1
)
y
0
=
α
,
n
=
0
,
1
,
2
,
⋯
,
N
−
1
\begin{cases} y_{n+1}=y_n+ \frac{h}{2}[k_1+k_2]\\ k_1=f(x_n,y_n)\\ k_2=f(x_n+h,y_n+hk_1)\\ y_0= \alpha ,n=0,1,2, \cdots ,N-1 \end{cases}
⎩
⎨
⎧yn+1=yn+2h[k1+k2]k1=f(xn,yn)k2=f(xn+h,yn+hk1)y0=α,n=0,1,2,⋯,N−1
假设前
n
{n}
n 步的计算没有截断误差,即当
y
n
=
y
(
x
n
)
{y_n=y(x_n)}
yn=y(xn) 时
k
1
=
f
(
x
n
,
y
n
)
=
f
n
k_1=f(x_n,y_n)=f_n
k1=f(xn,yn)=fn
应用二元Taylor展开
k
2
=
f
(
x
n
+
h
,
y
n
+
h
k
1
)
=
f
n
+
∂
f
n
∂
x
h
+
∂
f
n
∂
y
h
k
1
+
1
2
(
∂
2
f
n
∂
x
2
h
2
+
2
∂
2
f
n
∂
x
∂
y
h
2
k
1
+
∂
2
f
n
∂
y
2
h
2
k
1
2
)
+
O
(
h
3
)
\begin{align*} k_2=& f(x_n+h,y_n+hk_1) \\ \\ =&f_n+ \frac{\partial f_n}{\partial x}h+ \frac{\partial f_n}{\partial y}hk_1+ \frac{1}{2}( \frac{\partial^2 f_n}{\partial x^2} h^2+2 \frac{\partial^2 f_n}{\partial x\partial y}h^2k_1+ \frac{\partial^2 f_n}{\partial y^2}h^2k_1^2 ) + O(h^3) \end{align*}
k2==f(xn+h,yn+hk1)fn+∂x∂fnh+∂y∂fnhk1+21(∂x2∂2fnh2+2∂x∂y∂2fnh2k1+∂y2∂2fnh2k12)+O(h3)
将
k
1
{k_1}
k1 和
k
2
{k_2}
k2 代入
y
n
+
1
{y_{n+1}}
yn+1 ,然后整合按照
h
{h}
h 的升幂排列
y
n
+
1
=
y
n
+
h
f
n
+
h
2
2
(
∂
f
n
∂
x
+
∂
f
n
∂
y
f
n
)
+
h
3
4
(
∂
2
f
n
∂
x
2
+
2
∂
2
f
n
∂
x
∂
y
f
n
+
∂
2
f
n
∂
y
2
f
n
2
)
+
O
(
h
4
)
y_{n+1}=y_n+hf_n+ \frac{h^2}{2}(\frac{\partial f_n}{\partial x}+ \frac{\partial f_n}{\partial y}f_n)+ \frac{h^3}{4}( \frac{\partial^2 f_n}{\partial x^2}+ 2 \frac{\partial^2 f_n}{\partial x\partial y}f_n+ \frac{\partial^2 f_n}{\partial y^2}f_n^2)+ O(h^4)
yn+1=yn+hfn+2h2(∂x∂fn+∂y∂fnfn)+4h3(∂x2∂2fn+2∂x∂y∂2fnfn+∂y2∂2fnfn2)+O(h4)
应用一元Taylor展开
y
(
x
n
+
1
)
=
y
(
x
n
+
h
)
=
y
(
x
n
)
+
y
′
(
x
n
)
h
+
y
′
′
(
x
n
)
2
!
h
2
+
y
′
′
′
(
x
n
)
3
!
h
3
+
O
(
h
4
)
y(x_{n+1})=y(x_n+h)=y(x_n)+y'(x_n)h+ \frac{y''(x_n)}{2!}h^2+\frac{y'''(x_n)}{3!}h^3+O(h^4)
y(xn+1)=y(xn+h)=y(xn)+y′(xn)h+2!y′′(xn)h2+3!y′′′(xn)h3+O(h4)
将
y
′
(
x
n
)
,
y
′
′
(
x
n
)
,
y
′
′
′
(
x
n
)
{y'(x_n),y''(x_n),y'''(x_n)}
y′(xn),y′′(xn),y′′′(xn) 代入
y
(
x
n
+
1
)
{y(x_{n+1})}
y(xn+1) ,并按照
h
{h}
h 的升幂排列,有
y
(
x
n
+
1
)
=
y
n
+
h
f
n
+
h
2
2
(
∂
f
n
∂
x
+
∂
f
n
∂
y
f
n
)
+
h
3
6
(
∂
2
f
n
∂
x
2
+
2
∂
2
f
n
∂
x
∂
y
f
n
+
∂
2
f
n
∂
y
2
f
n
2
+
∂
f
n
∂
x
∂
f
n
∂
y
+
(
∂
f
n
∂
y
)
2
f
n
)
+
O
(
h
4
)
\begin{align*} y(x_{n+1})=&y_n+hf_n+\frac{h^2}{2}(\frac{\partial f_n}{\partial x}+ \frac{\partial f_n}{\partial y}f_n) \\ \\ &+\frac{h^3}{6}( \frac{\partial^2 f_n}{\partial x^2}+ 2 \frac{\partial^2 f_n}{\partial x\partial y}f_n+ \frac{\partial^2 f_n}{\partial y^2}f_n^2+ \frac{\partial f_n}{\partial x} \frac{\partial f_n}{\partial y}+(\frac{\partial f_n}{\partial y})^2f_n)+O(h^4) \end{align*}
y(xn+1)=yn+hfn+2h2(∂x∂fn+∂y∂fnfn)+6h3(∂x2∂2fn+2∂x∂y∂2fnfn+∂y2∂2fnfn2+∂x∂fn∂y∂fn+(∂y∂fn)2fn)+O(h4)
当
y
n
=
y
(
x
n
)
{y_n=y(x_n)}
yn=y(xn) 时,
y
n
+
1
{y_{n+1}}
yn+1 的表达式与精确解
y
(
x
n
+
1
)
{y(x_{n+1})}
y(xn+1) 的前三项完全相同,因此想间可得其局部截断误差
y
(
x
n
+
1
)
−
y
n
+
1
=
O
(
h
3
)
{y(x_{n+1})-y_{n+1}=O(h^3)}
y(xn+1)−yn+1=O(h3) 。
因此改进Euler方法是
2
{2}
2 阶方法。
[!example]-
已知求解常微分方程初值问题
{ y ′ = f ( x , y ) , x ∈ [ a , b ] y ( a ) = α \begin{cases} y'=f(x,y) ,&x\in[a,b]\\ y( a )= \alpha \end{cases} {y′=f(x,y),y(a)=αx∈[a,b]
的 差分公式
{ y n + 1 = y n + h 4 ( k 1 + 3 k 2 ) k 1 = f ( x n , y n ) k 2 = f ( x n + 2 3 h , y n + 2 3 h k 1 ) y 0 = α \begin{cases} y_{n+1}=y_n+ \frac{h}{4}(k_1+3k_2)\\ k_1=f(x_n,y_n) \\ k_2=f(x_n+ \frac{2}{3}h, y_n+ \frac{2}{3}hk_1)\\ y_0= \alpha \end{cases} ⎩ ⎨ ⎧yn+1=yn+4h(k1+3k2)k1=f(xn,yn)k2=f(xn+32h,yn+32hk1)y0=α
请问此差分公式是几阶方法?并写出截断误差的主项。
解:
k 1 = f n k_1=f_n k1=fn
k 2 = f n + ∂ f n ∂ x 2 3 h + ∂ f n ∂ y 2 3 h f n + 1 2 ( ∂ 2 f n ∂ x 2 4 9 h 2 + 2 ∂ 2 f n ∂ x ∂ y 4 9 h 2 f n + ∂ 2 f n ∂ y 2 4 9 h 2 f n 2 ) + O ( h 3 ) \begin{align*} k_2=&f_n+ \frac{\partial f_n}{\partial x} \frac{2}{3}h+ \frac{\partial f_n}{\partial y}\frac{2}{3}hf_n \\ \\ &+\frac{1}{2} ( \frac{\partial^2 f_n}{\partial x^2} \frac{4}{9} h^2+ 2 \frac{\partial^2 f_n}{\partial x\partial y} \frac{4}{9}h^2f_n + \frac{\partial^2 f_n}{\partial y^2} \frac{4}{9}h^2f_n^2 )+O(h^3) \end{align*} k2=fn+∂x∂fn32h+∂y∂fn32hfn+21(∂x2∂2fn94h2+2∂x∂y∂2fn94h2fn+∂y2∂2fn94h2fn2)+O(h3)
将 k 1 {k_1} k1 和 k 2 {k_2} k2 代入 y n + 1 {y_{n+1}} yn+1 ,然后整合按照 h {h} h 的升幂排列
y n + 1 = y n + f n h + h 2 2 ( ∂ f n ∂ x + ∂ f n ∂ y f n ) + h 3 6 ( ∂ 2 f n ∂ x 2 + 2 ∂ 2 f n ∂ x ∂ y f n + ∂ 2 f n ∂ y 2 f n 2 ) + O ( h 4 ) \begin{align*} y_{n+1}=&y_n+ f_nh+ \frac{h^2}{2}(\frac{\partial f_n}{\partial x}+ \frac{\partial f_n}{\partial y}f_n) \\ \\ &+ \frac{h^3}{6}(\frac{\partial^2 f_n}{\partial x^2}+ 2 \frac{\partial^2 f_n}{\partial x\partial y}f_n+ \frac{\partial^2 f_n}{\partial y^2}f_n^2 )+O(h^4) \end{align*} yn+1=yn+fnh+2h2(∂x∂fn+∂y∂fnfn)+6h3(∂x2∂2fn+2∂x∂y∂2fnfn+∂y2∂2fnfn2)+O(h4)
精确解
y ( x n + 1 ) = y ( x n ) + y ′ ( x n ) h + y ′ ′ ( x n ) 2 ! h 2 + y ′ ′ ′ ( x n ) 3 ! h 3 + O ( h 4 ) = y n + f n h + h 2 2 ( ∂ f n ∂ x + ∂ f n ∂ y f n ) + h 3 6 ( ∂ 2 f n ∂ x 2 + 2 ∂ 2 f n ∂ x ∂ y f n + ∂ 2 f n ∂ y 2 f n 2 + ∂ f n ∂ x ∂ f n ∂ y + ( ∂ f n ∂ y ) 2 f n ) + O ( h 4 ) \begin{align*} y(x_{n+1})=&y(x_n)+y'(x_n)h+ \frac{y''(x_n)}{2!}h^2+ \frac{y'''(x_n)}{3!}h^3+O(h^4) \\ \\ =& y_n+ f_nh+ \frac{h^2}{2}(\frac{\partial f_n}{\partial x}+ \frac{\partial f_n}{\partial y}f_n) \\ \\ &+ \frac{h^3}{6}(\frac{\partial^2 f_n}{\partial x^2}+ 2 \frac{\partial^2 f_n}{\partial x\partial y}f_n+ \frac{\partial^2 f_n}{\partial y^2}f_n^2 + \frac{\partial f_n}{\partial x} \frac{\partial f_n}{\partial y}+ (\frac{\partial f_n}{\partial y})^2f_n )+O(h^4) \end{align*} y(xn+1)==y(xn)+y′(xn)h+2!y′′(xn)h2+3!y′′′(xn)h3+O(h4)yn+fnh+2h2(∂x∂fn+∂y∂fnfn)+6h3(∂x2∂2fn+2∂x∂y∂2fnfn+∂y2∂2fnfn2+∂x∂fn∂y∂fn+(∂y∂fn)2fn)+O(h4)
y ( x n + 1 ) − y n + 1 = 1 6 ( ∂ f n ∂ x ∂ f n ∂ y + ( ∂ f n ∂ y ) 2 f n ) h 3 + O ( h 4 ) = O ( h 3 ) \begin{align*} y(x_{n+1})-y_{n+1}= \frac{1}{6}( \frac{\partial f_n}{\partial x} \frac{\partial f_n}{\partial y}+ (\frac{\partial f_n}{\partial y})^2f_n)h^3+ O(h^4)=O(h^3) \end{align*} y(xn+1)−yn+1=61(∂x∂fn∂y∂fn+(∂y∂fn)2fn)h3+O(h4)=O(h3)
所以是二阶方法。
阶段误差为
1 6 ( ∂ f n ∂ x ∂ f n ∂ y + ( ∂ f n ∂ y ) 2 f n ) h 3 \frac{1}{6}( \frac{\partial f_n}{\partial x} \frac{\partial f_n}{\partial y}+ (\frac{\partial f_n}{\partial y})^2f_n)h^3 61(∂x∂fn∂y∂fn+(∂y∂fn)2fn)h3
龙格-库塔(Runge-Kutta)公式
基本思想是设法计算f(x,y)在某些节点上的函数值,然后对这些函数值做线性组合,构造含待定参数的近似计算公式,再把近似计算公式和真实解的泰勒展开式相比较,并确定参数的取值一时的前面的若干项吻合,从而获得一定精度的计算公式。
改进欧拉方法是二段二阶龙格库塔公式中的一种。
高阶龙格库塔公式一步计算
f
(
x
,
y
)
{f(x,y)}
f(x,y) 的次数大于阶数,所以一般不使用。
常用龙格库塔公式如下:
三段三阶龙格库塔公式
{
y
n
+
1
=
y
n
+
1
6
(
k
1
+
4
k
2
+
k
3
)
k
1
=
h
f
(
x
n
,
y
n
)
k
2
=
h
f
(
x
n
+
0.5
h
,
y
n
+
0.5
k
1
)
k
3
=
h
f
(
x
n
+
h
,
y
n
−
k
1
+
2
k
2
)
\begin{cases} y_{n+1}=y_n+ \frac{1}{6}(k_1+4k_2+k_3) \\ \\ k_1=hf(x_n,y_n) \\ \\ k_2=hf(x_n+0.5h,y_n+0.5k_1) \\ \\ k_3=hf(x_n+h,y_n-k_1+2k_2) \end{cases}
⎩
⎨
⎧yn+1=yn+61(k1+4k2+k3)k1=hf(xn,yn)k2=hf(xn+0.5h,yn+0.5k1)k3=hf(xn+h,yn−k1+2k2)
标准四段四阶龙格库塔公式
{
y
n
+
1
=
y
n
+
1
6
(
k
1
+
2
k
2
+
2
k
3
+
k
4
)
k
1
=
h
f
(
x
n
,
y
n
)
k
2
=
h
f
(
x
n
+
0.5
h
,
y
n
+
0.5
k
1
)
k
3
=
h
f
(
x
n
+
0.5
h
,
y
n
+
0.5
k
2
)
k
4
=
h
f
(
x
n
+
h
,
y
n
+
k
3
)
\begin{cases} y_{n+1}=y_n+ \frac{1}{6}(k_1+2k_2+2k_3+k_4) \\ \\ k_1=hf(x_n,y_n) \\ \\ k_2=hf(x_n+0.5h,y_n+0.5k_1) \\ \\ k_3=hf(x_n+0.5h,y_n+0.5k_2) \\ \\ k_4 = hf(x_n+h,y_n+k_3) \end{cases}
⎩
⎨
⎧yn+1=yn+61(k1+2k2+2k3+k4)k1=hf(xn,yn)k2=hf(xn+0.5h,yn+0.5k1)k3=hf(xn+0.5h,yn+0.5k2)k4=hf(xn+h,yn+k3)
matlab实现
[x,y] = ode4RK(@f,1,0,1,10);
[x' y']
function r = f(x,y)
r = y-2*x./y;
end
%% 标准四阶四段龙格库塔公式
% 输入函数,初值,起点,终点,区间数
% 输出每个点与对应的函数值
function [x,y] = ode4RK(f,y0,a,b,n)
h = (b-a)/n;
x = linspace(a,b,n+1);
y(1) = y0;
for i = 1:n
k1 = h*f(x(i),y(i));
k2 = h*f(x(i)+0.5*h,y(i)+0.5*k1);
k3 = h*f(x(i)+0.5*h,y(i)+0.5*k2);
k4 = h*f(x(i)+h,y(i)+k3);
y(i+1) = y(i)+1/6*(k1+2*k2+2*k3+k4);
end
end