【数学知识】常微分方程的数值求解及python实现1_python 微分方程数值解

引言

在工程和科学技术的实际问题中，常需要求解微分方程。不管是本科生参加数学建模竞赛还是研究生进行自己的课题研究，常微分方程都是经常出现的。
然而，只有简单的和典型的微分方程可以求出解析解，而在实际问题中的微分方程往往无法求出解析解。

常微分方程认识

定义与常见形式

凡含有参数，未知函数和未知函数导数 (或微分) 的方程，称为微分方程，有时简称为方程，未知函数是一元函数的微分方程称作常微分方程，未知函数是多元函数的微分方程称作偏微分方程。
常微分方程定义：
$$F(x,y,y^{\prime },y^{\prime \prime },...,y^{(n)})=0$$
在高等数学中我们见过以下常微分方程：

• 初值问题(1)：
  { y ′ = f ( x , y ) a ≤ x ≤ b y ( a ) = y 0 \left\{

  $$\begin{array}{l} y'=f(x,y) \quad a\le x \le b \\ y(a) = y_0 \end{array}$$ \right. {y′=f(x,y)a≤x≤by(a)=y0​​

• 初值问题(2)：
  { y ′ ′ = f ( x , y , y ′ ) a ≤ x ≤ b y ( a ) = y 0 , y ′ ( a ) = α \left\{

  $$\begin{array}{l} y''=f(x,y,y') \quad a\le x \le b \\ y(a) = y_0 , y'(a)=\alpha \end{array}$$ \right. {y′′=f(x,y,y′)a≤x≤by(a)=y0​,y′(a)=α​

• 边值问题
  { y ′ ′ = f ( x , y , y ′ ) a ≤ x ≤ b y ( a ) = y 0 , y ′ ( b ) = y n \left\{

  $$\begin{array}{l} y''=f(x,y,y') \quad a\le x \le b \\ y(a) = y_0 , y'(b)=y_n \end{array}$$ \right. {y′′=f(x,y,y′)a≤x≤by(a)=y0​,y′(b)=yn​​

初值问题(1)的解存在的条件

• 存在和唯一性定理：如果连续函数f(x,y)关于y满足Lipschitz条件，即：
  $\exists$正数$L$，使得$\forall x\in [a,b]$，均有
  $|f(x,y_1)-f(x,y_2)|\le L|y_1-y_2|$
  则初值问题(1)的解存在且唯一

初值问题(1)的数值解

对于初值问题(1)：
{ y ′ = f ( x , y ) a ≤ x ≤ b y ( a ) = y 0 \left\{

$$\begin{array}{l} y'=f(x,y) \quad a\le x \le b \\ y(a) = y_0 \end{array}$$ \right. {y′=f(x,y)a≤x≤by(a)=y0​​
要求它的数值解，就是求未知函数$y(x)$在区间$[a,b]$上的一系列离散点(节点)
$$a=x_0\le x_1\le x_2\le \cdots \le x_n=b$$
上函数值$y(x_k)$的近似值$y_k(k=1,2,\cdots ,n)$，这就是初值问题(1)的数值解。
求它的数值解的关键在于$y^{\prime }(x)$的数值计算问题或者它的等价的积分方程$y(x)=y_0+{\int }_a^{x}f(t,y(t))dx$的数值计算问题。

求解微分方程的数值方法分为：数值微分和数值积分两种方法。

实际应用中，通常取求解区间$[a,b]$的等分点作为离散点，即
$$x_k=a+kh,k=0,1,\cdots ,n,whereh=\frac{n}{b-a}$$

推导初值问题的数值方法的途径：

• Talor展开
• 利用差商离散导数
• 利用数值积分方法

求初值问题数值解的方法是步进法，即从已知的初值$y_0$出发，通过一定的计算求$y_1$，然后由$y_1$或$y_0$和$y_1$求出$y_2$，依次计算到$y_n$，这时，可以分为

• 单步法：只利用$y_k$来计算$y_{k+1}$
• 多步法：计算$y_{k+1}$时，除了利用$y_k$，还需要用到已算出的若干$y_{k-j}$($j=1,2,\cdots ,l-1$)，并称为$l$步法

常微分方程数值求解方法

欧拉(尤拉) (Euler) 公式推导

方法一：利用Taylor展开

因为
$$x_{k+1}-x_k=h$$

故
$$y(x_{k+1})=y(x_k+h)$$

可将$y(x_{k+1})=y(x_k+h)$按照1阶泰勒公式展开，得
$$y(x_{k+1})=y(x_k)+\frac{1}{1!}{y}^{\prime }(x_k)(x_{k+1}-x_k)$$

也就是

$$y(x_{k+1})=y(x_k)+h{y}^{\prime }(x_k)$$
又根据初值条件
{ y ′ = f ( x , y ) a ≤ x ≤ b y ( a ) = y 0 \left\{

$$\begin{array}{l} y'=f(x,y) \quad a\le x \le b \\ y(a) = y_0 \end{array}$$ \right. {y′=f(x,y)a≤x≤by(a)=y0​​

可以得到数值解序列 { y k } \{y_k\} {yk​}的计算递推公式：
{ y 0 = y ( a ) y k + 1 = y k + h f ( x k , y k ) , k = 0 , 1 , ⋯   , n − 1 \left\{

$$\begin{array}{l} y_0 = y(a) \\ y_{k+1}=y_k+hf(x_k,y_k), k = 0,1,\cdots , n-1 \end{array}$$ \right. {y0​=y(a)yk+1​=yk​+hf(xk​,yk​),k=0,1,⋯,n−1​

方法二：利用差商离散导数

首先，导数的定义为：
$$f^{\prime }(x)=\frac{f(x)-f(x-\Delta x)}{\Delta x}$$

我们对$x_k$处的导数用差商来近似代替，如向前差商：
$$f^{\prime }(x_k)\approx \frac{y(x_{k+1})-y(x_k)}{h}$$

则微分方程初值问题可以化为：
{ y ( x k + 1 ) − y ( x k ) h ≈ f ( x k , y ( x k ) ) , k = 0 , 1 , ⋯   , n − 1 y ( a ) = y 0

$$\begin{cases} \frac{y(x_{k+1})-y(x_k)}{h} \approx f(x_k,y(x_k)) , k = 0,1,\cdots , n-1 \\ y(a) = y_0 \end{cases}$$ {hy(xk+1​)−y(xk​)​≈f(xk​,y(xk​)),k=0,1,⋯,n−1y(a)=y0​​

将近似号改为等号，精确解$y(x_k)$改为近似序列$y_k$，满足：
{ y 0 = y ( a ) y k + 1 = y k + h f ( x k , y k ) , k = 0 , 1 , ⋯   , n − 1 \left\{

$$\begin{array}{l} y_0 = y(a) \\ y_{k+1}=y_k+hf(x_k,y_k), k = 0,1,\cdots , n-1 \end{array}$$ \right. {y0​=y(a)yk+1​=yk​+hf(xk​,yk​),k=0,1,⋯,n−1​

第三种方法——利用数值积分方法这里就不介绍了
总的来说，上面推导的结果，就是欧拉(尤拉) (Euler) 公式，这里重写一遍：
{ y 0 = y ( a ) y k + 1 = y k + h f ( x k , y k ) , k = 0 , 1 , ⋯   , n − 1 \left\{

$$\begin{array}{l} y_0 = y(a) \\ y_{k+1}=y_k+hf(x_k,y_k), k = 0,1,\cdots , n-1 \end{array}$$ \right. {y0​=y(a)yk+1​=yk​+hf(xk​,yk​),k=0,1,⋯,n−1​

欧拉(尤拉) (Euler) 公式及其变形/改进

上一节详细推到了欧拉(尤拉) (Euler) 公式，这里介绍它的几种不同的形式：

• 显式欧拉公式：原始的欧拉公式
  { y 0 = y ( a ) y k + 1 = y k + h f ( x k , y k ) , k = 0 , 1 , ⋯   , n − 1 \left\{

  $$\begin{array}{l} y_0 = y(a) \\ y_{k+1}=y_k+hf(x_k,y_k), k = 0,1,\cdots , n-1 \end{array}$$ \right. {y0​=y(a)yk+1​=yk​+hf(xk​,yk​),k=0,1,⋯,n−1​

• 隐式欧拉公式：使用向后差商代替导数得到
  { y 0 = y ( a ) y k + 1 = y k + h f ( x k + 1 , y k + 1 ) , k = 0 , 1 , ⋯   , n − 1 \left\{

  $$\begin{array}{l} y_0 = y(a) \\ y_{k+1}=y_k+hf(x_{k+1},y_{k+1}), k = 0,1,\cdots , n-1 \end{array}$$ \right. {y0​=y(a)yk+1​=yk​+hf(xk+1​,yk+1​),k=0,1,⋯,n−1​

• 梯形公式：将显式欧拉公式和隐式欧拉公式作算术平均
  { y 0 = y ( a ) y k + 1 = y k + h 2 [ f ( x k , y k ) + f ( x k + 1 , y k + 1 ) ] , k = 0 , 1 , ⋯   , n − 1 \left\{

  $$\begin{array}{l} y_0 = y(a) \\ y_{k+1}=y_k+\frac{h}{2}[f(x_{k},y_{k})+f(x_{k+1},y_{k+1})], k = 0,1,\cdots , n-1 \end{array}$$ \right. {y0​=y(a)yk+1​=yk​+2h​[f(xk​,yk​)+f(xk+1​,yk+1​)],k=0,1,⋯,n−1​

• 预估-校正Euler方法：

即首先计算出初步的近似值$\overline{y_{k+1}}$:
$$\overline{y_{k+1}}=y_k+hf(x_k,y_k)$$

然后用该预估值替代梯形公式右端中的$y_{k+1}$进行计算，从而得到最后的校正值$y_{k+1}$：
$$y_{k+1}=y_k+\frac{h}{2}[f(x_k,y_k)+f(x_{k+1},\overline{y_{k+1}})]$$

欧拉(尤拉) (Euler) 公式Python实现

导入包，代码如下：

from pylab import *
import warnings
warnings.filterwarnings('ignore')
import matplotlib.pyplot as plt
%matplotlib inline
1
2
3
4
5

定义函数euler_xianshi(f,a=0,b=1,ya=1,h=0.1,print_result=True)使用显式欧拉公式进行微分方程数值计算代码如下：

def euler_xianshi(f,a=0,b=1,ya=1,h=0.1,print_result=True):
    '''显式欧拉公式, Explicit Euler formula
    f：需要求解的微分方程，y'=f(x,y)
        f的格式：
        def f(x,y):
            ...
            return dy
    a：求解区间起始值
    b：求解区间终止值
    ya：起始条件，ya=y(a)
    h：求解步长（区间[a,b]n等分）
    print_result：显示显式欧拉公式的每一步计算结果
    Returns
    ----------
    res：返回显式欧拉公式的求解结果
    '''
    res = []
    xi = a 
    yi = ya
    while xi<=b: # 在求解区间范围
        y = yi + h*f(xi,yi)
        if print_result:
            print('xi:{:.2f}, yi:{:.6f}'.format(xi,yi))
        res.append(y)
        xi, yi = xi+h, y
    
    return res

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
26
27
28

例1：求解下列微分方程

{ y ′ = x − y + 1 0 ≤ x ≤ 1 y ( 0 ) = 1 \left\{

$$\begin{array}{l} y'= x-y+1 \quad 0\le x \le 1 \\ y(0) = 1 \end{array}$$ \right. {y′=x−y+10≤x≤1y(0)=1​

• 第1步：定义函数表示求解的微分方程

def f(x,y):
    '''
    求解的微分方程：y'=x-y+1
    '''
    return x-y+1
1
2
3
4
5

• 第2步：应用前面定义的euler_xianshi函数进求解

a=0 # 求解区间起点
b=1 # 求解区间终点
y0=1    # 初值条件
h=0.05  # 步长
y = euler_xianshi(f,a,b,y0,h,print_result=True)
1
2
3
4
5

• 第3步：计算真实值

import math
import numpy as np
def fun_y(x):
    return x+np.exp(-x)

x_true = np.arange(a,b,0.1*h)
y_true = fun_y(x_true)
1
2
3
4
5
6
7

• 第4步：比较数值计算结果和真实值

x = np.arange(a,b,h)
plt.plot(x_true,y_true,label="true", color="red")
plt.scatter(x,y,label="predict")
plt.legend()
1
2
3
4

结果如下：

【作者简介】陈艺荣，男，目前在华南理工大学电子与信息学院广东省人体数据科学工程技术研究中心攻读博士，担任IEEE Access、IEEE Photonics Journal的审稿人。两次获得美国大学生数学建模竞赛(MCM)一等奖，获得2017年全国大学生数学建模竞赛(广东赛区)一等奖、2018年广东省大学生电子设计竞赛一等奖等科技竞赛奖项，主持一项2017-2019年国家级大学生创新训练项目获得优秀结题，参与两项广东大学生科技创新培育专项资金、一项2018-2019年国家级大学生创新训练项目获得良好结题，发表SCI论文4篇，授权实用新型专利8项，受理发明专利13项。
我的主页
我的Github
我的CSDN博客
我的Linkedin