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

最新推荐文章于 2024-07-25 22:24:14 发布

小小unicorn

最新推荐文章于 2024-07-25 22:24:14 发布

阅读量1.3k

点赞数 33

分类专栏：数学建模文章标签：数学建模学习

本文链接：https://blog.csdn.net/weixin_72066135/article/details/140270023

版权

数学建模专栏收录该内容

4 篇文章 1 订阅

订阅专栏

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

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

理论基础

导数与微分

微分和导数其实是紧密相关的概念。我们通常将导数理解为函数在某一点处切线的斜率。而微分则描述的是当我们对自变量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^{'} (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,y[:,0],label='y(0) Initial Function')
    l2, = plt.plot(t,y[:,1],label='y(1) The first order of Initial Function')
    l3, = plt.plot(t,y[:,2],label='y(2) The second order of Initial Function')
    plt.legend(handles=[l1,l2,l3])
    plt.grid('on')
    plt.subplot(212)
    l4, = plt.plot(y_.t, y_.y[0,:],'r--', label='y(0) Initial Function')
    l5,= plt.plot(y_.t,y_.y[1,:],'g--', label='y(1) The first order of Initial Function')
    l6, = plt.plot(y_.t,y_.y[2,:],'b-', label='y(2) The second order of Initial Function')
    plt.legend(handles=[l4,l5,l6]) # 显示图例
    plt.grid('on')
    plt.show()
solve_high_order_ode()

这里通过matplotlib.pyplot中提供的绘图接口绘制了两个数值解的图像。由于没有设置method参数，这里默认solve_ivp使用RK45（4-5阶 Runge-Kutta 法）方法进行求解。所解得的结果如图所示：
在这里插入图片描述
图中上半部分是使用odeint求解得到的结果，下半部分是由solve_ivp得到的结果，二者大差不差。一般来说对于普通的微分方程odeint与solve_ivp得到的结果差异不会太大，但有些情况下函数的微分容易发散，就会导致求解结果出现比较大的差异。odeint内置的原理也是4-5阶 Runge-Kutta 法，版本比较早所以求解也相对较为稳定。但solve_ivp则是后来新增的方法，有可能出现不太稳定的现象，内置方法较多所以也更加灵活。
Python求解微分方程组的模式有两种：一是采用基于基本原理自己写相关函数，这样操作比较繁琐，但是对于整个的求解过程会比较清晰明了；第二就是利用python下面的ode求解器，熟悉相关的输入输出，就可以完成数值求解。基于这个demo，在不同方向领域可以套用不同的微分方程组模型，进行仿真求解。但无论是常微分方程组还是偏微分方程组，使用的都是同一套思路，就是用差分代替微分。

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

在这里插入图片描述
这个例子和例2.2很像，但不同的是这是也一个非线性方程组。那么，输入就需要以数组的形式传入[x, y]两个函数的函数值，返回它们的导数值。这里使用solve_ivp对这个方程组进行求解如下：

def fun(t, w):
    x = w[0]
    y = w[1]
    return [-x**3-y,-y**3+x]
# 初始条件
y0 = [1,0.5]
yy = solve_ivp(fun, (0,100), y0, method='RK45',t_eval = np.arange(0,100,0.2) )
t = yy.t
data = yy.y
plt.plot(t, data[0, :])
plt.plot(t, data[1, :])
plt.xlabel("时间s")
plt.show()

绘制图像如图所示：
在这里插入图片描述
可以看到二者处于相互干扰的震荡状态，振幅随着时间的推移逐渐收敛并趋于稳定。如果在微分方程组里面还出现了高阶微分，我们又应该怎么做呢？来看下面这个例子：

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

在这里插入图片描述
这个例子就比较有趣了，在方程组里面涉及到了两个函数的导数怎么求解呢？本质上还是要使用换元法。完全可以令u=x’, v=y’，然后带入到方程组中把一个二元二阶微分方程组变为四元一阶微分方程组。

代码实现如下所示：

def fun(t, w):
    x = w[0]
    y = w[1]
    dx = w[2]
    dy = w[3]
    # 求导以后[x,y,dx,dy]变为[dx,dy,d2x,d2y]
    # d2x为w[2]，d2y为w[5]
    return [dx,dy,-dy-3*x+np.cos(2*t),4*dx-3*y+np.sin(2*t)]
# 初始条件
y0 = [0,0,1/5,1/6]
yy = solve_ivp(fun, (0,100), y0, method='RK45',t_eval = np.arange(0,100,0.2) )
t = yy.t
data = yy.y
plt.figure(figsize=(12,6))
plt.subplot(1,2,1)
plt.plot(t, data[0, :])
plt.plot(t, data[1, :])
plt.legend(['x','y'])
plt.xlabel("时间s")
plt.subplot(1,2,2)
plt.plot(t, data[2, :])
plt.plot(t, data[3, :])
plt.legend(["x' ","y' "])
plt.xlabel("时间s")
plt.show()

得到图像如图所示：
在这里插入图片描述
可以看到，图像呈现出一定的周期规律但并不是简谐运动。这样的方程解往往在物理学中有着实际意义，例如，这样的方程可以描述物体同时出现平动和摆动的过程中，位移-速度-加速度与角度-角速度-角加速度之间存在的关系。这样的例子曾出现在2022年全国大学生数学建模竞赛A题中。最后一个例子是对蝴蝶效应的求解。

例2.9 使用scipy求解洛伦兹系统的数值解，参数与初始值自设：

在这里插入图片描述

import numpy as np
from scipy.integrate import odeint
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
def dmove(Point, t, sets):
    p, r, b = sets
    x, y, z = Point
    return np.array([p*(y-x), x*(r-z)-y, x*y-b*z])
t = np.arange(0, 30, 0.001)
P1 = odeint(dmove, (0., 1., 0.), t, args=([10., 28., 8/3],))
P2 = odeint(dmove, (0., 1.01, 0.), t, args=([10., 28., 8/3],)) 
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot(P1[:,0], P1[:,1], P1[:,2], label='P1 trajectory')
ax.plot(P2[:,0], P2[:,1], P2[:,2], label='P2 trajectory')
ax.set_xlabel('X Axis')
ax.set_ylabel('Y Axis')
ax.set_zlabel('Z Axis')
ax.legend()
plt.show()

Mpl_toolkits.mplot3d提供了进行三维曲线、曲面绘制的函数，这里使用里面提供的三维坐标系绘制洛伦兹系统中点的运动轨迹。我们这里基于不同的初值绘制了两个点的轨迹P1和P2，并展示在图中：
在这里插入图片描述
可以看到，曲线的形状呈现双螺旋状，有些像蝴蝶的翅膀。所以洛伦兹系统又被叫做“蝴蝶效应”。蝴蝶效应本质上就是指，即使给这个系统的初始值一点微小的变化，曲线的形状也会出现很大不同。仅仅是把 $y$ 改变了 $0.01$ ，曲线的密集程度与蝴蝶翅膀的大小也是有所不同的。这是个混沌系统里面的典型案例。

偏微分方程的数值求解

偏微分方程是针对多元函数来说的，它在物理学中有着很深刻的现实意义。但是，偏微分方程往往比常微分方程更难求解，并且Python也没有提供偏微分方程的专用工具包。怎么求解这些偏微分方程呢？我们要始终把握一个思想：就是把连续的问题离散化。这一节会通过一系列的物理案例来看到Python如何求解一些典型的偏微分方程。

理论基础

偏微分方程实际上就是由多元函数、自变量与多元函数的偏导数及高阶偏导数之间构成的方程。它在工程中很多地方都有深刻应用，比如波动力学、热学、电磁学等。我们常研究的就是二元函数的二阶偏微分方程，其基本形式为：
在这里插入图片描述
在方程中，如果𝐴𝐵𝐶三个常系数不全为0，定义判别式Δ=B² −4AC，当判别式大于0称其为双曲线式方程；若判别式等于0，则称其为抛物线式方程；若判别式小于0，则称其为椭圆式方程。

刚刚我们说到，二阶偏微分方程主要有三类：椭圆方程，抛物方程和双曲方程。双曲方程描述变量以一定速度沿某个方向传播，常用于描述振动与波动问题。椭圆方程描述变量以一定深度沿所有方向传播，常用于描述静电场、引力场等稳态问题。抛物方程描述变量沿下游传播，常用于描述热传导和扩散等瞬态问题。它们都在工程中有实际应用。

偏微分方程的定解问题通常很难求出解析解，只能通过数值计算方法对偏微分方程的近似求解。常用偏微分方程数值解法有包括有限差分方法、有限元方法、有限体方法、共轭梯度法，等等。在使用这些数值方法时通常先对问题的求解区域进行网格剖分，然后将定解问题离散为代数方程组，求出在离散网格点上的近似值。

偏微分方程的典型应用有很多。描述热源传热过程中温度变化的热传导方程本质上是一个抛物类微分方程。大名鼎鼎的韦东奕大神所着重研究的纳维-斯托克斯方程，所描述的是流体流速与流体密度、压力、外阻力之间的关系，在机械工程、能源工程等制造领域有着重要应用。还有电磁场中非常重要的麦克斯韦方程组，本质上也是偏微分方程。下面我们会根据一系列的具体应用来看如何去构建与求解微分方程。

应用案例：

在这里插入图片描述

例2.10

使用偏微分方程对RC电路进行建模，分析电容放电过程中电量随时间的变化。
在这里插入图片描述
$\left\{ \begin{align} IR + \frac{Q}{C} &= 0, \\ I &= \frac{\mathrm{d}Q}{\mathrm{d}t}, \end{align} \right. \tag{1}$

那么我们为什么还要来谈这个案例？我想通过这个案例来给大家讲述一下把一个连续问题离散化的方法。我们先从一阶微分方程的计算原理说起。将一阶微分方程离散化，实际上也就是把它写成迭代、差分的形式。对于一个连续的问题有
在这里插入图片描述
画出𝑦(𝑡)的图像如图所示：

在图中，如果求解区间为 $x_0, x_0+h]$ ，已知点 $x_{0},f(x_{0}))$ ，要求解 $f(x_0+h)$ 应该怎么做？根据微分方程 $f (x, y) = 0$ 可以很容易地解得每个点的导数值。在距离 $h$ 很小的情况下，我们说，切线与割线是可以逼近的，也就可以用切线的对应值逼近函数值。但在上图中，我们发现：如果使用 $x_0$ 点处的导数作为切线斜率，那么得到的估计值是比实际值要小的；但如果使用 $x_0+h$ 点的导数作为斜率，那么估计值比实际值又大一些。能不能取一个折中的方案呢？很简单，把二者进行一个平均就可以了： $y_{n+1} = y_{n} + \frac{f(y_{n}, t_{n}) + f(y_{n+1}, t_{n+1})}{2} h. \tag{2.3.5}$ 那么对于例2.10中的电容放电过程，这个问题的代码可以这样写：

import numpy as np
import matplotlib.pyplot as plt
rc = 2.0 #设置常数
dt = 0.5 #设置步长
n = 1000 #设置分割段数
t = 0.0 #设置初始时间
q = 1.0 #设置初始电量
#先定义三个空列表
qt=[] #用来盛放差分得到的q值
qt0=[] #用来盛放解析得到的q值
time = [] #用来盛放时间值
for i in range(n):
    t = t + dt
    q1 = q - q*dt/rc #qn+1的近似值
    q = q - 0.5*(q1*dt/rc + q*dt/rc) #差分递推关系
    q0 = np.exp(-t/rc) #解析关系
    qt.append(q) #差分得到的q值列表
    qt0.append(q0) #解析得到的q值列表
    time.append(t) #时间列表
plt.plot(time,qt,'o',label='Euler-Modify') #差分得到的电量随时间的变化
plt.plot(time,qt0,'r-',label='Analytical') #解析得到的电量随时间的变化
plt.xlabel('time')
plt.ylabel('charge')
plt.xlim(0,20)
plt.ylim(-0.2,1.0)
plt.legend(loc='upper right')
plt.show()

这个案例我们没有用任何包里面的微分方程求解器，纯手写的情况下解了这个微分方程。它的结果如图所示：
在这里插入图片描述
这个方法被称为Euler 法，是在求解常微分方程中的一种常见数值方法。

例2.11

一维热传导方程是一个典型的抛物型二阶偏微分方程。设𝑢(𝑥,𝑡)表示在时间𝑡,空间𝑥处的温度，则根据傅里叶定律(单位时间内流经单位面积的热量和该处温度的负梯度成正比),可以导出热传导方程：
在这里插入图片描述
其中𝜆称为热扩散率，𝑘,𝐶,𝑝分别为热导率，比热和质量密度，是由系统本身确定的常量。问题的形式为：

请求解这个问题的数值解。

一元函数的微分方程可以绘制曲线，那么二元函数的偏微分方程应该就可以绘制曲面。那么，怎么对这个问题进行离散化呢？对于一个一般的二阶抛物型偏微分方程：
在这里插入图片描述
对时间和空间的一阶偏导数微分是容易离散化的：

那么，对于空间的二阶微分，可以看作是一阶微分的再微分：

了解了这个原理，我们将例2.11中的边界条件和迭代规则进行翻译如下：

import numpy as np
import matplotlib.pyplot as plt
h = 0.1#空间步长
N =30#空间步数
dt = 0.0001#时间步长
M = 10000#时间的步数
A = dt/(h**2) #lambda*tau/h^2
U = np.zeros([N+1,M+1])#建立二维空数组
Space = np.arange(0,(N+1)*h,h)#建立空间等差数列，从0到3，公差是h
#边界条件
for k in np.arange(0,M+1):
    U[0,k] = 0.0
    U[N,k] = 0.0
#初始条件
for i in np.arange(0,N):
    U[i,0]=4*i*h*(3-i*h)
#递推关系
for k in np.arange(0,M):
    for i in np.arange(1,N):
        U[i,k+1]=A*U[i+1,k]+(1-2*A)*U[i,k]+A*U[i-1,k]

将解空间抽象为以时间为横坐标、空间为纵坐标的网格，翻译时间与空间的边界条件，对于网格内每个点使用二重循环遍历每个点，根据差分后的迭代方程进行演化。不同时刻的温度随空间坐标的变化图像如下：

plt.plot(Space,U[:,0], 'g-', label='t=0',linewidth=1.0)
plt.plot(Space,U[:,3000], 'b-', label='t=3/10',linewidth=1.0)
plt.plot(Space,U[:,6000], 'k-', label='t=6/10',linewidth=1.0)
plt.plot(Space,U[:,9000], 'r-', label='t=9/10',linewidth=1.0)
plt.plot(Space,U[:,10000], 'y-', label='t=1',linewidth=1.0)
plt.ylabel('u(x,t)', fontsize=20)
plt.xlabel('x', fontsize=20)
plt.xlim(0,3)
plt.ylim(-2,10)
plt.legend(loc='upper right')
plt.show()

在这里插入图片描述
在图中可以看到，随着时间的推进，温度分布呈现出一种动态变化的过程。在 $t = 0$ 时（绿线），我们看到起始的温度分布情况；随着时间的推移，图中的曲线显示出温度在不同位置的变化。

整体上，这些曲线描述了温度如何随着时间从初始状态演变到一个稳态分布。这种分析对于理解热传导、扩散过程，以及如何在时间上控制温度分布都是非常有用的。图形的y轴显示的是𝑢(𝑥,𝑡)即位置𝑥在时间𝑡的温度，而x轴表示空间坐标。曲线下方较深的颜色表示较低的温度值，而曲线顶部较浅的颜色表示较高的温度值。通过设置坐标轴的范围和图例的位置，这张图为观察者提供了清晰的数据解读。

在这里插入图片描述

#温度等高线随时空坐标的变化，温度越高，颜色越偏红
extent = [0,1,0,3] #时间和空间的取值范围
levels = np.arange(0,10,0.1)#温度等高线的变化范围0-10，变化间隔为0.1
plt.contourf(U,levels,origin='lower',extent=extent,cmap=plt.cm.jet)
plt.ylabel('x', fontsize=20)
plt.xlabel('t', fontsize=20)
plt.show()

在这里插入图片描述
从图中可以看到的是一个温度分布的热力图，其中颜色的变化表明了不同温度的区域。温度越高，颜色越偏向红色，温度较低的区域则显现为蓝色。这种热力图通常用于显示温度如何在空间内分布以及如何随时间变化。图中的 $\mathrm{y}$ 轴（标记为 $x$ ）代表空间坐标，而 $\mathrm{x}$ 轴（标记为 $t$ ）代表时间。热力图覆盖的范围是时间从 $0$ 到 $1$ ，空间从 $0$ 到 $3$ 。可以清晰地看到，在图的左侧（时间较早）温度整体较低，而在图的右侧（时间较晚）温度较高，这表示随着时间的推移，整体温度有所上升。
等温线的密集区表示温度变化较大的区域，而等温线的稀疏区则表示温度变化较小的区域。从热力图中我们可以推断，最高温区域集中在图的右上角，而最低温区域则在左下角。通过色彩的深浅变化，我们可以直观地看到温度在空间中如何变化以及时间对这种分布的影响。

例2.12

平流过程是大气运动中重要的过程。平流方程（Advection equation）描述某一物理量的平流作用而引起局地变化的物理过程，最简单的形式是一维平流方程。
在这里插入图片描述
式中𝑢为某物理量，𝑣为系统速度，𝑥为水平方向分量，𝑡为时间。该方程可以求得解析解：

考虑一维线性平流偏微分方程的数值解法，采用有限差分法求解。简单地，采用一阶迎风格式的差分方法（First-order Upwind)，一阶导数的差分表达式为：
在这里插入图片描述
代码实现如下：

import numpy as np
import matplotlib.pyplot as plt
# 初始条件函数 U(x,0)
def funcUx_0(x, p): 
    u0 = np.sin(2 * (x-p)**2)
    return u0
# 输入参数
v1 = 1.0  # 平流方程参数，系统速度
p = 0.25  # 初始条件函数 u(x,0) 中的参数
tc = 0  # 开始时间
te = 1.0  # 终止时间: (0, te)
xa = 0.0  # 空间范围: (xa, xb)
xb = np.pi
dt = 0.02  # 时间差分步长
nNodes = 100  # 空间网格数
# 初始化
nsteps = round(te/dt)
dx = (xb - xa) / nNodes
x = np.arange(xa-dx, xb+2*dx, dx)
ux_0 = funcUx_0(x, p)
u = ux_0.copy()  # u(j)
ujp = ux_0.copy()  # u(j+1)
# 时域差分
for i in range(nsteps):
    plt.clf()  # 清除当前 figure 的所有axes, 但是保留当前窗口
    # 计算 u(j+1)
    for j in range(nNodes + 2):
        ujp[j] = u[j] - (v1 * dt/dx) * (u[j] - u[j-1])
    # 更新边界条件
    u = ujp.copy()
    u[0] = u[nNodes + 1]
    u[nNodes+2] = u[1]
    tc += dt
    # 绘图
plt.plot(x, u, 'b-', label="v1= 1.0")
plt.axis((xa-0.1, xb + 0.1, -1.1, 1.1))
plt.xlabel("x")
plt.ylabel("U(x)")
plt.legend(loc=(0.05,0.05))
plt.show()

在这里插入图片描述
从图中可以发现，函数 $U (x)$ 显示了随空间 $x$ 变化的波动性质。这个波动可能代表一维平流方程在某一特定时间 $t$ 的数值解。图中蓝色的线表示速度 $v_{1} = 1.0$ 下的解，而波形的变化暗示了初态条件 $\sin \big[ 2(x-p)^{2} \big]$ 随时间的演化。
注意到曲线在x轴的不同位置出现了波峰和波谷，这表明了函数值随位置的变化并非均匀。由于是一阶迎风格式的数值解法，我们可能会观察到与理论解相比有一定程度的数值扩散或者数值耗散，这是由于一阶方法在数值传输过程中的固有特性。

在这个示例中，曲线的形状可能表示了经过一段时间演化后，平流作用在初始条件
𝑈(𝑥,0)上的效果。如果初始波形移动的速度是𝑣₁，那么这个图形可能表明波形随着时间的推进而向右移动。然而，由于是一阶迎风差分，我们也可以预期波形会有一定程度的变形，这在实际中表现为波峰变得不那么尖锐，以及波形整体变得更加平坦。这是数值方法的离散化误差导致的结果。

这个数值解可以帮助理解平流方程在数值模拟中的行为，特别是当理论解难以获得时，数值方法提供了一种有效的途径来近似解决实际问题。

例2.13

波动方程（wave equation）是典型的双曲偏微分方程，广泛应用于声学，电磁学，和流体力学等领域，描述自然界中的各种的波动现象，包括横波和纵波，例如声波、光波和水波。考虑如下二维波动方程的初边值问题：
在这里插入图片描述
式中：𝑢是振幅；𝑐为波的传播速率，𝑐可以是固定常数，或位置的函数 𝑐(𝑥,𝑦)也可以是振幅的函数 𝑐(𝑢)

考虑二维波动偏微分方程的数值解法，采用有限差分法求解。简单地，采用迎风法的三点差分格式，将上述的偏微分方程离散为差分方程：
在这里插入图片描述
可以得到迭代规则为：

代码实现：

import numpy as np
import matplotlib.pyplot as plt
# 模型参数
c = 1.0  # 波的传播速率
tc, te = 0.0, 1.0  # 时间范围，0<t<te
xa, xb = 0.0, 1.0  # 空间范围，xa<x<xb
ya, yb = 0.0, 1.0  # 空间范围，ya<y<yb
# 初始化
c2 = c*c  # 方程参数
dt = 0.01  # 时间步长
dx = dy = 0.02  # 空间步长
tNodes = round(te/dt)  # t轴 时序网格数
xNodes = round((xb-xa)/dx)  # $\mathrm{x}$轴 空间网格数
yNodes = round((yb-ya)/dy)  # $\mathrm{y}$轴 空间网格数
tZone = np.arange(0, (tNodes+1)*dt, dt)  # 建立空间网格
xZone = np.arange(0, (xNodes+1)*dx, dx)  # 建立空间网格
yZone = np.arange(0, (yNodes+1)*dy, dy)  # 建立空间网格
xx, yy = np.meshgrid(xZone, yZone)  # 生成网格点的坐标 xx,yy (二维数组)
# 步长比检验(r>1 则算法不稳定)
r = 4 * c2 * dt*dt / (dx*dx+dy*dy)
print("dt = {:.2f}, dx = {:.2f}, dy = {:.2f}, r = {:.2f}".format(dt,dx,dy,r))
assert r < 1.0, "Error: r>1, unstable step ratio of dt2/(dx2+dy2) ."
rx = c*c * dt**2/dx**2
ry = c*c * dt**2/dy**2
# 绘图
fig = plt.figure(figsize=(8,6))
ax1 = fig.add_subplot(111, projection='3d')
# 计算初始值
U = np.zeros([tNodes+1, xNodes+1, yNodes+1])  # 建立三维数组
U[0] = np.sin(6*np.pi*xx)+np.cos(4*np.pi*yy)  # U[0,:,:]
U[1] = np.sin(6*np.pi*xx)+np.cos(4*np.pi*yy)  # U[1,:,:]
surf = ax1.plot_surface(xx, yy, U[0,:,:], rstride=2, cstride=2, cmap=plt.cm.coolwarm)
# 有限差分法求解
for k in range(2,tNodes+1):
    for i in range(1,xNodes):
        for j in range(1,yNodes):
            U[k,i,j] = rx*(U[k-1,i-1,j]+U[k-1,i+1,j]) + ry*(U[k-1,i,j-1]+U[k-1,i,j+1])\
                     + 2*(1-rx-ry)*U[k-1,i,j] -U[k-2,i,j]
surf = ax1.plot_surface(xx, yy, U[k,:,:], rstride=2, cstride=2, cmap='rainbow')
ax1.set_xlim3d(0, 1.0)
ax1.set_ylim3d(0, 1.0)
ax1.set_zlim3d(-2, 2)
ax1.set_title("2D wave equationt (t= %.2f)" % (k*dt))
ax1.set_xlabel("x")
ax1.set_ylabel("y")
plt.show()

这个函数是一个三元函数，实际上是可以做出一个曲面随时间变化的动画的。大家可以尝试使用matplotlib提供的动画功能进行绘制，这里展示其中一个瞬时状态：
在这里插入图片描述

例2.14

热传导方程（heat equation）是典型的抛物形偏微分方程，也成为扩散方程。广泛应用于声学，电磁学，和流体力学等领域，描述自然界中的各种的波动现象，包括横波和纵波，例如声波、光波和水波。之前的例2.11我们已经看到了一维热传导方程的求解，现在考虑如下二维热传导方程的初边值问题：
在这里插入图片描述
类似的，同上一个例子中的有限差分法，我们将这个方程离散化可以得到：

事实上，这个方程还有一种矩阵形式：

其中

我们这次就用这个矩阵形式去进行偏微分方程的求解：

import numpy as np
import matplotlib.pyplot as plt
def showcontourf(zMat, xyRange, tNow):  # 绘制等温云图
    x = np.linspace(xyRange[0], xyRange[1], zMat.shape[1])
    y = np.linspace(xyRange[2], xyRange[3], zMat.shape[0])
    xx,yy = np.meshgrid(x,y)
    zMax = np.max(zMat)
    yMax, xMax = np.where(zMat==zMax)[0][0], np.where(zMat==zMax)[1][0]
    levels = np.arange(0,100,1)
    showText = "time = {:.1f} s\nhotpoint = {:.1f} C".format(tNow, zMax)
    plt.plot(x[xMax],y[yMax],'ro')  # 绘制最高温度点
    plt.contourf(xx, yy, zMat, 100, cmap=plt.cm.get_cmap('jet'), origin='lower', levels = levels)
    plt.annotate(showText, xy=(x[xMax],y[yMax]), xytext=(x[xMax],y[yMax]),fontsize=10)
    plt.colorbar()
    plt.xlabel('X')
    plt.ylabel('Y')
    plt.title('Temperature distribution of the plate')
    plt.show()
# 模型参数
uIni = 25  # 初始温度值
uBound = 25.0  # 边界条件
c = 1.0  # 热传导参数
qv = 50.0  # 热源功率
x_0, y0 = 0.0, 3.0  # 热源初始位置
vx, vy = 2.0, 1.0  # 热源移动速度
# 求解范围
tc, te = 0.0, 5.0  # 时间范围，0<t<te
xa, xb = 0.0, 16.0  # 空间范围，xa<x<xb
ya, yb = 0.0, 12.0  # 空间范围，ya<y<yb
# 初始化
dt = 0.002  # 时间步长
dx = dy = 0.1  # 空间步长
tNodes = round(te/dt)  # t轴 时序网格数
xNodes = round((xb-xa)/dx)  # $\mathrm{x}$轴 空间网格数
yNodes = round((yb-ya)/dy)  # $\mathrm{y}$轴 空间网格数
xyRange = np.array([xa, xb, ya, yb])
xZone = np.linspace(xa, xb, xNodes+1)  # 建立空间网格
yZone = np.linspace(ya, yb, yNodes+1)  # 建立空间网格
xx,yy = np.meshgrid(xZone, yZone)  # 生成网格点的坐标 xx,yy (二维数组)
# 计算 差分系数矩阵 A、B (三对角对称矩阵)，差分系数 rx,ry,ft
A = (-2) * np.eye(xNodes+1, k=0) + (1) * np.eye(xNodes+1, k=-1) + (1) * np.eye(xNodes+1, k=1)
B = (-2) * np.eye(yNodes+1, k=0) + (1) * np.eye(yNodes+1, k=-1) + (1) * np.eye(yNodes+1, k=1)
rx, ry, ft = c*dt/(dx*dx), c*dt/(dy*dy), qv*dt
# 计算 初始值
U = uIni * np.ones((yNodes+1, xNodes+1))  # 初始温度 u0
# 前向Euler 法一阶差分求解
for k in range(tNodes+1):
    t = k * dt  # 当前时间
    # 热源条件
    # (1) 恒定热源：Qv(x_0,y0,t) = qv, 功率、位置 恒定
    # Qv = qv
    # (2) 热源功率随时间变化 Qv(x_0,y0,t)=f(t)
    # Qv = qv*np.sin(t*np.pi) if t<2.0 else qv
    # (3) 热源位置随时间变化 Qv(x,y,t)=f(x(t),y(t))
    xt, yt = x_0+vx*t, y0+vy*t  # 热源位置变化
    Qv = qv * np.exp(-((xx-xt)**2+(yy-yt)**2))  # 热源方程
    # 边界条件
    U[:,0] = U[:,-1] = uBound
    U[0,:] = U[-1,:] = uBound
    # 差分求解
    U = U + rx * np.dot(U,A) + ry * np.dot(B,U) + Qv*dt
    if k % 100 == 0:
        print('t={:.2f}s\tTmax={:.1f}  Tmin={:.1f}'.format(t, np.max(U), np.min(U)))
showcontourf(U, xyRange, k*dt)  # 绘制等温云图

这段代码只需要把showcontourf放在循环体内即可实现动画效果，可以清晰看到热源在空间中的分布。我这里展示最后状态下的温度分布图：

在这里插入图片描述
最终状态下的温度分布图显示了热源随时间在平板中移动的情况，并将热量传递给周围材料。图像中的红点标记了最高温度点，即‘热点’，对应于最后时间步骤中热源的当前位置。颜色渐变代表温度分布，红色是最热的区域，蓝色是最冷的。等温线表示温度相等的水平线。

热点周围的温度梯度平滑，这表明热量通过平板均匀扩散。这样的模拟在许多应用中都很有用，例如在散热器设计、理解制造过程中的热梯度，甚至在诸如地热源传播热量的自然现象中。

例2.15

椭圆偏微分方程是一类重要的偏微分方程，主要用来描述物理的平衡稳定状态，如定常状态下的电磁场、引力场和反应扩散现象等，广泛应用于流体力学、弹性力学、电磁学、几何学和变分法中。考虑如下二维 Poisson 方程：
在这里插入图片描述
这个方程怎么解呢？考虑二维椭圆偏微分方程的数值解法，采用有限差分法求解。简单地，采用五点差分格式表示二阶导数的差分表达式，将上述的偏微分方程离散为差分方程：

椭圆型偏微分描述不随时间变化的均衡状态，没有初始条件，因此不能沿时间步长递推求解。对上式的差分方程，可以通过矩阵求逆方法求解，但当
ℎ
h较小时网格很多，矩阵求逆的内存占用和计算量极大。于是，可以使用迭代松弛法递推求得二维椭圆方程的数值解。假定𝑥和𝑦的间距都为ℎ，松弛系数为𝑤，则迭代过程可以表示为：

在这里插入图片描述
考虑一个特殊情况，也就是当𝑓(𝑥,𝑦)的情况下，迭代松弛法的代码如下：

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
# 求解范围
xa, xb = 0.0, 1.0  # 空间范围，xa<x<xb
ya, yb = 0.0, 1.0  # 空间范围，ya<y<yb
# 初始化
h = 0.01  # 空间步长, dx = dy = 0.01
w = 0.5  # 松弛因子
nodes = round((xb-xa)/h)  # $\mathrm{x}$轴 空间网格数
# 边值条件
u = np.zeros((nodes+1, nodes+1))
for i in range(nodes+1):
    u[i, 0] = 1.0 + np.sin(0.5*(i-50)/np.pi)
    u[i, -1] = -1.0 + 0.5*np.sin((i-50)/np.pi)
    u[0, i] = -1.0 + 0.5*np.sin((i-50)/np.pi)
    u[-1, i] = 1.0 + np.sin(0.5*(50-i)/np.pi)
# 迭代松弛法求解
for iter in range(100):
    for i in range(1, nodes):
        for j in range(1, nodes):
            u[i, j] = w/4 * (u[i-1, j] + u[i+1, j] + u[i, j-1] + u[i, j+1]) + (1-w) * u[i, j]
# 绘图
x = np.linspace(0, 1, nodes+1)
y = np.linspace(0, 1, nodes+1)
xx, yy = np.meshgrid(x, y)
fig = plt.figure(figsize=(8,6))
ax = fig.add_subplot(111, projection='3d')
surf = ax.plot_surface(xx, yy, u, cmap=plt.get_cmap('rainbow'))
fig.colorbar(surf, shrink=0.5)
ax.set_xlim3d(0, 1.0)
ax.set_ylim3d(0, 1.0)
ax.set_zlim3d(-2, 2.5)
ax.set_title("2D elliptic partial differential equation")
ax.set_xlabel("X")
ax.set_ylabel("Y")
plt.show()

在这里插入图片描述
从图中可以看到，解呈现出一系列波峰和波谷，这与边界条件的正弦函数有关。每个波峰和波谷都是颜色映射中的热点和冷点，分别代表着方程解的局部最大值和最小值。波峰出现在和轴的中间位置，而波谷则出现在四个角和边缘。颜色的变化代表了解的幅度，从红色的高值到蓝色的低值。整体上，解的形状表现出了椭圆方程特有的对称性和周期性。

sympy中的pdsolve给出了一些简单偏微分方程的解析解法，但sympy目前不支持二阶偏微分方程的求解。pdsolve的调用格式形如pdsolve(eq, func=None, hint=‘default’, dict=False, solvefun=None, **kwargs)，具体使用我们可以看到下面的例子：

例2.16

使用sympy中的pdsolve解下面这个偏微分方程：
在这里插入图片描述
可以给出如下代码

from sympy.solvers.pde import pdsolve
from sympy import Function, pprint, exp
from sympy.abc import x,y
f = Function('f')
eq = -2*f(x,y).diff(x) + 4*f(x,y).diff(y) + 5*f(x,y) - exp(x + 3*y)
pdsolve(eq)