【数学建模学习手册】第二章：微分方程（一）

本专栏内容为：数学建模原理 记录学习数学建模

????博主****个人主页：小小unicorn
⏩专栏分类：数学建模
????代码仓库：小小unicorn的代码仓库????
????????????关注我带你学习编程知识

微分方程

• 理论基础
• • 导数与微分
  • 一阶线性微分方程的解
  • 二阶常系数线性微分方程的解
  • 利用Python求函数的微分与积分
• 使用Scipy和Sympy解微分方程
• • 使用sympy求解微分方程解析解
  • • 例2.1 使用sympy解下面这个微分方程：
    • 例2.2 使用sympy解下面这个常微分方程组：
  • 使用scipy求解微分方程数值解
  • • 例2.3 使用scipy解下面这个微分方程的数值解：
    • 例2.4 使用scipy解下面这个微分方程的数值解：
    • 例2.5 使用scipy解下面这个高阶微分方程的数值解：
    • 例2.6 使用scipy解下面这个高阶微分方程的数值解：
    • 例2.7 使用scipy解下面这个微分方程组的数值解：
    • 例2.8 使用 scipy 解下面这个微分方程组的数值解：
    • 例2.9 使用scipy求解洛伦兹系统的数值解，参数与初始值自设：
• 偏微分方程的数值求解
• • 理论基础
  • 应用案例：
  • • 例2.10
    • 例2.11
    • 例2.12
    • 例2.13
    • 例2.14
    • 例2.15
    • 例2.16

理论基础

导数与微分

微分和导数其实是紧密相关的概念。我们通常将导数理解为函数在某一点处切线的斜率。而微分则描述的是当我们对自变量x施加一个非常小的增量d????时，函数值相应的变化量与之间的关系。当d????非常小的时候，函数的变化量就接近于在该点处切线的变化量d????.因此，我们可以用这种方式来理解微分：

在下图中，我们展示了函数、导数和微分之间的关系。微分实际上描述的是点????处切线的斜率；导数则描述的是割线MN的斜率。但当d????足够小的时候，切线的斜率和割线的斜率就会非常接近，这就是微分的核心概念。而微分方程，就是描述函数与其导数之间关系的方程。

相对于求微分，我们还有求积分的概念。积分本质上是根据已知的导数反推出原函数，这就是不定积分。而定积分则是在反推出原函数后，还需要计算该函数在特定区间内的值的差异。通常情况下，我们可以通过查阅常见函数的导数表来进行微分和不定积分的计算。

注意：割线斜率等于切线斜率的前提是dx非常小，这是一种极限思想的体现。虽然它们之间存在一个无穷小量PN的差距，但当我们在考虑dx时，这种差异就可以忽略不计了。这就是微分和积分的基本思想。

一阶线性微分方程的解

一阶线性微分方程描述的是怎么一回事呢？它是指形如下方的方程：

这里的$y$是一个未知数，而$P$和$Q$是已知的函数。我们的目标是找出y的解，即他的通解方式。为了解这个方程，我们通常会使用分离变量积分法和常数变易法。首先，我们先尝试解一个特殊情况的齐次方程，即当$Q(x)=0$时：

通过分离变量，我们得到：

接着，对两边进行不定积分，我们可以得到解的通式为：

其中C是一个常数，但在一般情况下，$Q(x)$不一定为0，所以我们需要将C换成一个函数$C(x)$,然后对$y$求导，并将其带入原方程中求的$C(x)$的通解。这就是常数变易法。进一步可推导出方程的通解为：

其中$C$为常数。

注意：这里的定积分符号用于求原函数。这就是为什么我们在高中学习的积分符号应该按照这种方式书写的原因。齐次方程指的是方程右边等于0的情况，而非齐次方程则是方程右边不恒等于0的情况。解非齐次方程更具有一般性，但很多非齐次方程的解也是基于齐次方程的解进行拓展的。

二阶常系数线性微分方程的解

二阶常系数线性微分方程可以表示为：

这个方程关联了二阶导数、一阶导数和函数本身。解决这个方程的一般策略是先考虑对应的齐次方程，即让$C(c)$为0：

解这种二阶常系数齐次线性微分方程时，我们通常使用特征根法。这个方法的关键是求解特征方程：

这个齐次方程的解的形式取决于特征方程的根。根据特征方程的不同实根、相同实根、或共轭复根，齐次微分方程的解会有不同的形式：

注意：这里为什么二次方程的根与齐次方程的解之间会有联系，这正是数学之美的体现之一。如果想检验这个方程的解是否正确，实际上并不难，可以使用 Vieta 定理将 ???? 和 ???? 代入，将两个方程统一起来，再通过换元法将其降为一阶微分方程进行验证。

对于一般的二阶非齐次线性微分方程，我们可以根据右侧????(????)的形式推导出一个特解。非齐次方程的通解等于齐次方程的通解加上非齐次方程的特解。求微分方程的特解有时需要观察法，但幸运的是，存在两种特殊形式：

其中P_m(x)是一个m次多项式，Q_n(x)是一个n次多项式，这两种方式的特解分别为：

其中????的取值取决于特征方程根的个数：如果有两个不同的实根，则k=2；如果有两个相同的实根，则k=1；如果没有实根，则k=0。通过上述形式，我们可以解出二阶线性微分方程。

特征根法和“特解+通解”的策略不仅适用于二阶线性微分方程，也适用于一般的高阶线性微分方程。只要特征方程是多项式，它至少满足韦达定理。在后续的差分方程中，特征根法同样会发挥重要作用。

利用Python求函数的微分与积分

在Python中，我们可以使用Numpy和SciPy这两个库来进行函数的微分和积分计算。下面将通过具体示例来说明如何使用这些库来求解函数的微分和积分。 假设我们需要计算函数f(x) = cos(2πx) * exp(-x) + 1.2在区间[0, 0.7]上的定积分。我们可以使用SciPy库中的quad函数来完成这个任务：

import numpy as np
from scipy.integrate import quad
# 定义函数
def f(x):
    return np.cos(2 * np.pi * x) * np.exp(-x) + 1.2
# 计算定积分
integral, error = quad(f, 0, 0.7)
print(f'定积分的结果是：{integral}')
# 定积分的结果是：0.7951866427656943 

除了使用SciPy库中的quad函数求解定积分外，我们还可以使用数值积分的方法来近似计算。一种常见的数值积分方法是梯形法则。

下面我们将通过一个示例来说明如何使用梯形法则来近似计算函数的定积分。 假设我们需要计算函数f(x) = cos(2πx) * exp(-x) + 1.2在区间[0, 0.7]上的定积分。我们可以使用梯形法则来近似求解：

h=x[1]-x[0]
xn=0.7
s=0
for i in range(1000):
    xn1=xn+h
    yn=np.cos(2*np.pi*xn)*np.exp(-xn)+1.2
    yn1=np.cos(2*np.pi*xn1)*np.exp(-xn1)+1.2
    s0=(yn+yn1)*h/2
    s+=s0
    xn=xn1
s
# 24.31183595181452

对于函数的微分，我们可以使用Numpy库中的gradient函数来近似求解。例如，我们想要求解函数f(x) = x^2在点x = 1处的导数：

#计算导数
import numpy as np
# 定义x的取值范围和步长
x = np.linspace(0, 2, 100)
y = x**2
# 计算导数
dydx = np.gradient(y, x)
# 在x=1处的导数值
derivative_at_1 = dydx[np.argmin(abs(x - 1))]
print(f'在x=1处的导数值是：{derivative_at_1}')
# 在x=1处的导数值是：1.9797979797979792

使用Scipy和Sympy解微分方程

前面我们见过了求微分方程解析解的一些方法，我们知道，微分方程的解本质上是通过给定函数与微分之间的关系求解出函数的表达式。但是事实上，大多数微分方程是没有解析解的，也就是无法求解出函数的具体解析式。这是不是意味着这样的微分方程不可解呢？也不尽然。

在上一章中我们已经知道，以前我们难以求解的超越方程也是可能给出数值解的，那么微分方程是否也会存在数值解呢？

使用sympy求解微分方程解析解

我们此前介绍的一阶、二阶常系数线性微分方程通解的形式就是一种解析解，但在科学与工程实际中我们遇到的微分方程形式会比这些基本形式更为复杂，条件也更多。事实上多数情况下，大多数微分方程其实是求不出解析解的，只能在不同取值条件下求一个数值解。那么如何编写算法去求数值解才能使精度尽可能提高呢？数值解会随着初始条件而变化，怎么变化呢？函数值又与自变量之间怎么变化呢？

在回答这些问题之前，请让我们先了解一番：如何使用python求解微分方程的解析解呢？但凡涉及到符号运算，通常都是使用sympy库实现。

Sympy是一个数学符号运算库。能解决积分、微分方程等各种数学运算方法，用起来也是很简单，效果可以和Matlab媲美。其中内置的Sympy.dsolve方法是解微分方程解析解的一种良好方式，而对于有初始值的微分方程问题，我们通常在求出其通解形式后通过解方程组的方法得到参数。

这个方法通过声明符号变量的方式求得最优解。

例如，我们看下面这个例子：

例2.1 使用sympy解下面这个微分方程：

若使用sympy，我们首先要声明两个符号变量，其中变量y是变量x的函数。代码如下：

from sympy import *
y = symbols('y', cls=Function)
x = symbols('x')
eq = Eq(y(x).diff(x,2)+2*y(x).diff(x,1)+y(x), x*x)
## y''+4y'+29y=0
print(dsolve(eq, y(x)))

这段代码通过sympy中的symbols类创建两个实例化的符号变量x和y，在y中我们通过cls参数声明y是一个scipy.Function对象（也就是说，y是一个函数）。

表达微分方程解析解的方法是通过创建一个Eq对象，这个对象分别存储方程左右两边。其中，y(x).diff(x,2)表明y是x的函数，然后需要取函数对x的2阶导数。最后，若想求解函数y的解析式，只需要调用dsolve(eq,y(x))函数即可。代码返回结果：

可以看到，代码能够给出完整的解析式。之所以还保留了参数C1和C2是因为在求解过程中没有给微分方程指定初值。

我们再来看一个例子，这个例子是使用sympy解一个常微分方程组：

例2.2 使用sympy解下面这个常微分方程组：

这个方程组里面的????₁,????₂,????₃都是关于????的函数，所以需要声明四个符号变量。不同的是，在这里每个函数都指定了初始值，并且三个函数的导数高度相关，该怎么描述这种相关呢？我们来看下面的例子：

t=symbols('t')
x1,x2,x3=symbols('x1,x2,x3',cls=Function)
eq=[x1(t).diff(t)-2*x1(t)+3*x2(t)-3*x3(t),
    x2(t).diff(t)-4*x1(t)+5*x2(t)-3*x3(t),
    x3(t).diff(t)-4*x1(t)+4*x2(t)-2*x3(t)]
con={x1(0):1, x2(0):2, x3(0):3}
s=dsolve(eq,ics=con)
print(s)

sympy当中内置的symbols工具是可以通过字符串批量创建变量的，这为我们带来了很大的方便。如果需要求解的是一个方程组，则使用列表将每一个方程表达出来即可。这里我们采取了不创建对象的方式，而是直接将方程组移项使每个方程右侧都为0。通过字典的方式保存函数的初始值，并利用ics参数传入dsolve从而得到方程的解。

结果返回的是一个Eq对象构成的列表，每个对象代表了一个函数的解析式。对于这个例子，大家可以发现：它是一个线性的微分方程组，而针对线性方程我们还可以使用矩阵的形式去表示。所以，这个问题还有第二种写法：

x=Matrix([x1(t),x2(t),x3(t)])
A=Matrix([[2,-3,3],[4,-5,3],[4,-4,2]])
eq=x.diff(t)-A*x
s=dsolve(eq,ics={x1(0):1, x2(0):2, x3(0):3})
print(s)

通过sympy中内置的符号矩阵Matrix对象构造函数向量和系数矩阵，通过对方程组矩阵化也可以得出一样的结果。返回值同上。使用sympy中的符号函数绘图得到结果如下：

from sympy.plotting import plot
from sympy import *
t=Symbol('t')
plot(2*exp(2*t) - exp(-t), line_color='red')
plot(2*exp(2*t) - exp(-t) + exp(-2*t), line_color='blue')
plot(2*exp(2*t) + exp(-2*t), line_color='green')

求解图如下：



sympy通过plotting下面的plot功能可以进行一些符号函数的绘图，但每一次调用都会创建一个独立的图窗，难以在同一张图上绘制多个函数的曲线。若要绘制多个函数则需要使用matplotlib来完成。

使用scipy求解微分方程数值解

微分方程的数值解是什么样子的呢？虽然大多数微分方程没有解析解，但解析式也并不是唯一可以表示函数的形式。函数的表示还可以用列表法和作图法来表示，而微分方程的数值解也正是像列表一样针对自变量数组中的每一个取值给出相对精确的因变量值。

Python求解微分方程数值解可以使用scipy库中的integrate包。在这当中有两个重要的函数：odeint和solve_ivp。但本质上，从底层来讲求解微分方程数值解的核心原理都是Euler 法和Runge-Kutta 法。关于这两个方法，我们会在后面进行进一步探讨。

我们先来了解一下odeint的用法吧。odeint()函数需要至少三个变量，第一个是微分方程函数，第二个是微分方程初值，第三个是微分的自变量。为了具体了解它的用法，我们通过一个例子来分析：

例2.3 使用scipy解下面这个微分方程的数值解：

首先需要通过def语句或者lambda表达式定义微分方程的表达式，然后定义微分方程的初值。代码如下：

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
dy=lambda y,x: 1/(1+x**2)-2*y**2 # y'=1/(1+x^2)-2y^2
'''
def dy(y,x):
    return 1/(1+x**2)-2*y**2
'''
x=np.arange(0,10.5,0.1) #从0开始，每次增加0.1，到10.5为止（取不到10.5）
sol=odeint(dy,0,x) # odeint输入：微分方程dy，y的首项（y(0)等于多少），自变量列表
print("x={}\n对应的数值解y={}".format(x,sol.T))
plt.plot(x,sol)
plt.show()

这里odeint函数传入的三个参数分别是函数表达式，函数的初值与自变量。自变量是一个数组，通过numpy.arange生成一个范围在[0, 10.5)的等差数列，公差为0.1。返回的结果sol是针对数组x中每个值的对应函数值，可以通过matplotlib.pyplot绘图得到函数的结果。函数的图像如图所示：

我们再来看一个例子，这个例子是一个不可积函数的积分问题：

例2.4 使用scipy解下面这个微分方程的数值解：

仿照例2.3中的代码，这个问题可以改写为：

def dy_dt(y,t):
    return np.sin(t**2)
y0=[1]
t = np.arange(-10,10,0.01)
y=odeint(dy_dt,y0,t)
plt.plot(t, y)
plt.show()

得到的结果必然是一个奇函数，图像为：

刚刚两个例子都是讲述了一阶微分方程如何求解，那么二阶及以上的高阶微分方程如何求解呢？事实上，Python求解微分方程数值解的时候是无法直接求解高阶微分方程的，必须通过换元降次的方法实现低阶化，把一个高阶微分方程替换成若干个一阶微分方程组成的微分方程组才能求解。

具体的，我们可以看下面这个例子：

例2.5 使用scipy解下面这个高阶微分方程的数值解：

这很显然是个二阶微分方程，并且不是常系数所以不能直接给出解析解。为了给这个方程做降次，
式子就可以代换为：

对于微分方程组，我们传入[y,u]两个函数的原函数值，返回的函数值为[y’,u’]。所以，只需要对每个微分表达式给出解析形式就可以了。代码如下：

# odeint是通过把二阶微分转化为一个方程组的形式求解高阶方程的
# y''=20(1-y^2)y'-y
def fvdp(y,t):
    '''
    要把y看出一个向量，y = [dy0,dy1,dy2,...]分别表示y的n阶导，那么
    y[0]就是需要求解的函数，y[1]表示一阶导，y[2]表示二阶导，以此类推
    '''
    dy1 = y[1]      # y[1]=dy/dt，一阶导                     y[0]表示原函数
    dy2 = 20*(1-y[0]**2) * y[1] - y[0]                    # y[1]表示一阶微分
    # y[0]是最初始，也就是需要求解的函数
    # 注意返回的顺序是[一阶导， 二阶导]，这就形成了一阶微分方程组
    return [dy1, dy2] 
    
# 求解的是一个二阶微分方程，所以输入的时候同时输入原函数y和微分y'
# y[0]表示原函数， y[1]表示一阶微分
# dy1表示一阶微分， dy2表示的是二阶微分
# 可以发现，dy1和y[1]表示的是同一个东西
# 把y''分离变量分离出来： dy2=20*(1-y[0]**2)*y[1]-y[0]
def solve_second_order_ode():
    '''
    求解二阶ODE
    '''
    x = np.arange(0,0.25,0.01)#给x规定范围
    y0 = [0.0, 2.0] # 初值条件
    # 初值[3.0, -5.0]表示y(0)=3,y'(0)=-5
    # 返回y，其中y[:,0]是y[0]的值，就是最终解，y[:,1]是y'(x)的值
    y = odeint(fvdp, y0, x)
    
    y1, = plt.plot(x,y[:,0],label='y')
    y1_1, = plt.plot(x,y[:,1],label='y‘')             
    plt.legend(handles=[y1,y1_1])   #创建图例
    
    plt.show()
solve_second_order_ode()

定义函数fvdp，传入y的原函数值和一阶导数值（列表传入），返回y的一阶导数值和二阶导数值。初值条件y(0)=0和y'(0)=2传入odeint函数中，自变量是取值[0, 0.25)的一个等距数组。解得的y其实包含两列，第一列是函数值，第二列是导数值。结果的图像如下。

图2.1.5展示的是原函数$y(x)$与一阶导数$y^{\prime }(x)$的图像。从图像中可以看到，原函数$y(x)$呈现出一种振荡衰减的趋势，随着$x$的增加，$y(x)$的振幅逐渐减小，最终趋于稳定。这是因为二阶微分方程中的非线性项起到了阻尼作用，当的绝对值接近$1$时，该项的值变小，从而减弱了$y$的增长速率，导致振荡的衰减。
同时，一阶导数的图像显示出与原函数相似的振荡衰减模式，但相比之下，其变化更加剧烈。这是因为直接受到非线性阻尼项的影响，而则是间接受到影响。

总的来说，这个微分方程组描述了一个非线性阻尼振荡系统，其解的行为随着初始条件和时间的变化而发生变化。在这个例子中，初始条件(0)=2导致了一个振荡衰减的解，这种解在物理学和工程学中很常见，用于描述许多实际系统的动态行为。

我们再来看一个更高阶函数的求解的案例。

例2.6 使用scipy解下面这个高阶微分方程的数值解：

这个案例当然可以和上面一样如法炮制，输入[y, y', y'']返回[y', y'', y''']。这里再次介绍一个案例是想引出Python求微分方程数值解的另一个函数solve_ivp的用法。

首先，仍然是通过换元法对函数进行定义：

def f(t,y):
    dy1 = y[1]
    dy2 = y[2]
    dy3 = -y[0]+dy1-dy2-np.cos(t)
    return [dy1,dy2,dy3]

Solve_ivp函数的用法与odeint非常类似，只不过比odeint多了两个参数。一个是t_span参数，表示自变量的取值范围；另一个是method参数，可以选择多种不同的数值求解算法。常见的内置方法包括RK45, RK23, DOP853, Radau, BDF等多种方法，通常使用RK45多一些。它的使用方法与odeint对比起来很类似，对这个问题进行代码实现如下：

def solve_high_order_ode():
    '''
    求解高阶ODE
    '''
    t = np.linspace(0,6,1000)
    tspan = (0.0, 6.0)
    y0 = [0.0, pi, 0.0]
    # 初值[0,1，0]表示y(0)=0,y'(0)=1，y''(0)=0
    # 返回y， 其中y[:,0]是y[0]的值 ，就是最终解 ，y[:,1]是y'(x)的值
    y = odeint(f, y0, t)
    y_ = solve_ivp(f,t_span=tspan, y0=y0, t_eval=t)
    plt.subplot(211)
    l1, = plt.plot(t