【数值分析】常微分方程的数值解，欧拉公式，梯形公式，龙格库塔公式，matlab实现_数值解和公式解

常微分方程初边值问题的数值解法

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$ 是 $x$ 和 $y$ 的已知函数， $\alpha$ 为给定的初值。

Lipschitz利普希兹条件
若函数 $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$ 满足Lipschitz条件：
$$|f(x,y)-f(x,\bar{y})|\le L|y-\bar{y}|,\forall y,\bar{y}$$
其中 $L>0$ 为Lipschiitz常数，则初值问题（在 $a$ 的某邻域内）存在唯一解。此时Lipschitz常数不必小于 $1$ 。

[!example]-
$$\frac{\mathrm{d}y}{\mathrm{d}x}=1+x\sin (xy),x\in [0,2]$$
$$y(0)=1$$
解：对任意 $y$ ， $\bar{y}$ ，对变量 $y$ 应用微分中值定理，存在 $y$ ，使得
$$\frac{f(x,y)-f(x,\bar{y})}{y-\bar{y}}=\frac{\partial }{\partial y}f(x,y)=x^2\cos (xy)$$
于是有
$$|f(x,y)-f(x,\bar{y})|=|x^2\cos (xy)||y-\bar{y}|\le 4|y-\bar{y}|$$
因而此时 $f(x,y)$ 关于 $y$ 满足Lipschitz条件，且常数 $L=4$ 。因此从理论上讲，初值问题存在惟一解

---

差分公式的格式

Euler公式

$$y_{n+1}=y_n+hf(x_n,y_n)$$
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
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

梯形公式

$$y_{n+1}=y_n+\frac{h}{2}(f(x_n,y_n)+f(x_{n+1},y_{n+1}))$$
matlab实现，由于方程右边用到了左边的值，这里通过不动点迭代出 $y_{n+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
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25

Euler中点公式

$$y_{n+1}=y_{n-1}+2hf(x_n,y_n)$$
课本上没提，这里不写了，接手这篇笔记的人可以补上。

改进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)$
$$\{\begin{array}{l}{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{array}$$
等价二阶二段龙格库塔格式
$$\{\begin{array}{l}{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+h{k}_1)\\ y_0=\alpha ,n=0,1,2,\cdots ,N-1\end{array}$$
改进欧拉公式精度高于欧拉公式，低于梯形公式，但比起梯形公式是显式格式，计算量较小，相对折中。
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
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19

---

局部截断误差 $y(x_{n+1})-y_{n+1}$

局部截断误差的阶数
若单步差分公式的局部截断误差为 $O(h^{P+1})$ ，则称该公式为 $P$ 阶方法。
Taylor公式
一元Taylor公式：
$$y(x_{n+1})=y(x_n+h)=y(x_n)+y^{\prime }(x_n)h+\frac{y^{\prime \prime }(x_n)}{2!}{h}^2+\frac{y^{\prime \prime \prime }(x_n)}{3!}{h}^3+O(h^4)$$
二元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^{\prime }(x_n),y^{\prime \prime }(x_n),y^{\prime \prime \prime }(x_n)$ 的写法
因 $y^{\prime }(x)=f(x,y(x))$ ，所以若 $y(x_n)=y_n$ ，则有
$$\begin{array}{rl}{y}^{\prime }(x_n)=& f(x_n,y(x_n))=f(x_n,y_n)=f_n\\ & \\ y^{\prime \prime }(x_n)=& \frac{\partial f(x_n,y_n)}{\partial x}+\frac{\partial f(x_n,y_n)}{\partial y}{y}^{\prime }(x_n)=\frac{\partial f_n}{\partial x}+\frac{\partial f_n}{\partial y}{f}_n\\ & \\ y^{\prime \prime \prime }(x_n)=& \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}{)}^2{f}_n\end{array}$$
关于改进Euler方法的局部截断误差分析
改进Euler方法格式：
$$\{\begin{array}{l}{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+h{k}_1)\\ y_0=\alpha ,n=0,1,2,\cdots ,N-1\end{array}$$
假设前 $n$ 步的计算没有截断误差，即当 $y_n=y(x_n)$ 时
$$k_1=f(x_n,y_n)=f_n$$
应用二元Taylor展开
$$\begin{array}{rl}{k}_2=& f(x_n+h,y_n+h{k}_1)\\ & \\ =& f_n+\frac{\partial f_n}{\partial x}h+\frac{\partial f_n}{\partial y}h{k}_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}^2{k}_1+\frac{{\partial }^2{f}_n}{\partial y^2}{h}^2{k}_1^2)+O(h^3)\end{array}$$
将 $k_1$ 和 $k_2$ 代入 $y_{n+1}$ ，然后整合按照 $h$ 的升幂排列
$$y_{n+1}=y_n+h{f}_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)$$
应用一元Taylor展开
$$y(x_{n+1})=y(x_n+h)=y(x_n)+y^{\prime }(x_n)h+\frac{y^{\prime \prime }(x_n)}{2!}{h}^2+\frac{y^{\prime \prime \prime }(x_n)}{3!}{h}^3+O(h^4)$$
将 $y^{\prime }(x_n),y^{\prime \prime }(x_n),y^{\prime \prime \prime }(x_n)$ 代入 $y(x_{n+1})$ ，并按照 $h$ 的升幂排列，有
$$\begin{array}{rl}y(x_{n+1})=& y_n+h{f}_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}{)}^2{f}_n)+O(h^4)\end{array}$$
当 $y_n=y(x_n)$ 时， $y_{n+1}$ 的表达式与精确解 $y(x_{n+1})$ 的前三项完全相同，因此想间可得其局部截断误差 $y(x_{n+1})-y_{n+1}=O(h^3)$ 。
因此改进Euler方法是 $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]
的 差分公式
$$\{\begin{array}{l}{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}h{k}_1)\\ y_0=\alpha \end{array}$$
请问此差分公式是几阶方法？并写出截断误差的主项。
解：
$$k_1=f_n$$
$$\begin{array}{rl}{k}_2=& f_n+\frac{\partial f_n}{\partial x}\frac{2}{3}h+\frac{\partial f_n}{\partial y}\frac{2}{3}h{f}_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}^2{f}_n+\frac{{\partial }^2{f}_n}{\partial y^2}\frac{4}{9}{h}^2{f}_n^2)+O(h^3)\end{array}$$
将 $k_1$ 和 $k_2$ 代入 $y_{n+1}$ ，然后整合按照 $h$ 的升幂排列
$$\begin{array}{rl}{y}_{n+1}=& y_n+f_{n}h+\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{array}$$
精确解
$$\begin{array}{rl}y(x_{n+1})=& y(x_n)+y^{\prime }(x_n)h+\frac{y^{\prime \prime }(x_n)}{2!}{h}^2+\frac{y^{\prime \prime \prime }(x_n)}{3!}{h}^3+O(h^4)\\ & \\ =& y_n+f_{n}h+\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}{)}^2{f}_n)+O(h^4)\end{array}$$
$$\begin{array}{r}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}{)}^2{f}_n)h^3+O(h^4)=O(h^3)\end{array}$$
所以是二阶方法。
阶段误差为
$$\frac{1}{6}(\frac{\partial f_n}{\partial x}\frac{\partial f_n}{\partial y}+(\frac{\partial f_n}{\partial y}{)}^2{f}_n)h^3$$

---

龙格-库塔（Runge-Kutta）公式

基本思想是设法计算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
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21

---

下链

---