2023年6月29日发(作者：)

常微分⽅程的解析解（⽅法归纳）以及基于Python的⼆阶微分⽅程边值问题的数值算例实现常微分⽅程的解析解(⽅法归纳)以及基于Python的微分⽅程数值解算例实现本⽂归纳常见常微分⽅程的解析解解法以及基于Python的微分⽅程数值解算例实现。⽂章⽬录常微分⽅程的解析解考虑常微分⽅程的解析解法，我们⼀般可以将其归纳为如下⼏类：1. 可分离变量的微分⽅程（⼀阶）2. ⼀阶齐次（⾮齐次）线性微分⽅程（⼀阶）3. ⼆阶常系数微分⽅程（⼆阶）4. ⾼阶常系数微分⽅程（n阶）可分离变量的微分⽅程（⼀阶）这类微分⽅程可以变形成如下形式：f(x)dx=g(y)dy两边同时积分即可解出函数，难点主要在于不定积分，是最简单的微分⽅程。某些⽅程看似不可分离变量，但是经过换元之后，其实还是可分离变量的，不要被这种⽅程迷惑。⼀阶齐次（⾮齐次）线性微分⽅程（⼀阶）形如dydx+P(x)y=Q(x)的⽅程叫做⼀阶线性微分⽅程，若Q(x)为0，则⽅程齐次，否则称为⾮齐次。解法： (直接套公式)y(x)=e−∫P(x)dx(∫e∫P(x)dxQ(x)dx+C)伯努利⽅程形如dy1dx+P(x)y=Q(x)yn,n∈R,n=的⽅程称为伯努利⽅程，这种⽅程可以通过以下步骤化为⼀阶线性微分⽅程：y−ndydx+P(x)y1−n=Q(x)1dy1−n1−n⋅dx+P(x)y1−n=Q(x)令

y1−n=u, ⽅程两边同时乘以

1−n，得到dudx+(1−n)P(x)u=(1−n)Q(x)即

dx+P′(x)u=Q′(x).这就将伯努利⽅程归结为可以套公式的⼀阶线性微分⽅程。du⼆阶常系数微分⽅程（⼆阶）形如y′′+py′+qy=f(x)的⽅程称为⼆阶常系数微分⽅程，若f(x)≡0，则⽅程称为齐次的，反之称为⾮齐次的。以下默认⽅程是⾮齐次的。求解此类⽅程分两步：1. 求出齐次通解2. 求出⾮齐次特解原⽅程的解=齐次通解+⾮齐次特解齐次通解的求法⾸先假设

f(x)≡0.⽤特征⽅程法，写出对应的特征⽅程并且求解：λ2+pλ+q=0解的情况分为以下三种：情况⼀：⽅程有两个不同的实数解假设两个实数解分别是

λ1、λ2, 此时⽅程的通解是Y(x)=C1eλ1x+C2eλ2x情况⼆：⽅程有⼀个⼆重解假设该解等于λ，此时⽅程的通解是Y(x)=(C1+C2x)eλx情况三：⽅程有⼀对共轭复解假设这对解是

α±iβ, 此时⽅程的通解是Y(x)=eαx(C1cos(βx)+C2sin(βx))⾮齐次特解的求法对于

f(x) 和特征根的情况，对特解的情况做如下归纳：1.

f(x)=Pm(x)，其中

Pm(x) 表⽰

x 的最⾼次数为

m 的多项式。若0不是⽅程特征解则⽅程有特解

y∗=Qm(x)若0是⽅程的单特征解则⽅程有特解

y∗=xQm(x)若0是⽅程的⼆重特征解则⽅程有特解

y∗=x2Qm(x)其中

Qm(x)=b0+b1x+…+bmxm,

bi(i=0,1,…,m)是需要带回原⽅程来确定的系数。2.

f(x)=eαxPm(x)若

α 不是⽅程特征解则⽅程有特解

y∗=eαxQm(x)若

α 是⽅程的单特征解则⽅程有特解

y∗=xeαxQm(x)若

α 是⽅程的⼆重特征解则⽅程有特解

y∗=x2eαxQm(x)3. $f(x)=e^{alpha x}(a_{1}cos(beta x)+a_{2}sin(beta x)) $若

α±iβ 不是特征解则⽅程有特解

y∗=eαx(A1cos(βx)+A2sin(βx))若

α±iβ 是特征解则⽅程有特解

y∗=xeαx(A1cos(βx)+A2sin(βx))其中

A1、A2 是需要带回原⽅程来确定的系数。⾼阶常系数微分⽅程（n阶）形如y(n)+p1y(n−1)+...+pn−1y′+pny=f(x)的⽅程叫做⾼阶常系数微分⽅程，若

f(x)≡0，则⽅程是齐次的，否则是⾮齐次的。下⾯默认⽅程是⾮齐次的。求解此类⽅程分两步：1. 求出齐次通解2. 求出⾮齐次特解原⽅程的解=齐次通解+⾮齐次特解齐次通解的求法(参考⼆阶常系数微分⽅程解法)⾮齐次特解的求法(参考⼆阶常系数微分⽅程解法)算例考虑带有第三类边界条件的⼆阶常系数微分⽅程边值问题y′′(x)+2y′(x)−3y(x)=3x+1,⎧⎪y′(1)−y(1)=1,⎨x∈[1,2]′⎪⎩y(2)−y(2)=−4.1. 求出上述两点边值问题的解析解;2. 通过有限差分⽅法算出其数值解;(取h=5)并计算误差、绘图、输出1.0,1.2,1.4,1.6,1.8,2.0处的数值解.1问题⼀:两点边值问题的解析解由于此⽅程是⾮齐次的,故求解此类⽅程分两步：1. 求出齐次通解2. 求出⾮齐次特解原⽅程的解=齐次通解+⾮齐次特解齐次通解的求法⾸先假设

y′′(x)+2y′(x)−3y(x)≡0. ⽤特征⽅程法，写出对应的特征⽅程λ2+2λ−3=0求解得到两个不同的实数特征根:λ1=1,λ2=−3.此时⽅程的齐次通解是y(x)=C1eλ1x+C2eλ2x⾮齐次特解的求法由于

λ1=0,λ2=0. 所以⾮齐次特解形式为y∗=ax+b将上式代⼊控制⽅程有2a−3ax−3b=3x+1求解得:

a=b=−1, 即⾮齐次特解为

y∗=−x−1.原⽅程的解=齐次通解+⾮齐次特解于是,原⽅程的全解为y(x)=C1ex+C2e−3x−x−1因为该问题给出的是第三类边界条件,故需要求解的导函数y′(x)=C1ex−3C2e−3x−1且有y(1)=C1e+C2e−3−2⎧⎪⎪⎪⎪y′(1)=C1e−3C2e−3−1⎨y(2)=C1e2+C2e−6−3⎪⎪⎪⎪y′(2)=C1e2−3C2e−6−1⎩将以上各式代⼊边界条件y′(1)−y(1)=−4C2e−3+1=1{y′(2)+y(2)=2C1e2−2C2e−6−4=−4解此⽅程组可得:

C1=C2=0.综上所述,原两点边值问题的解为y=−x−1常微分⽅程的数值解⼀般⼆阶线性常微分⽅程边值问题的差分法对⼀般的⼆阶微分⽅程边值问题⎧y′′(x)+p(x)y′(x)+q(x)y(x)=f(x),⎨a1y(a)+a2y′(a)=a⎩β1y(b)+β2y′(b)=βa

I=[a,b]N 等分, 即得到区间

I=[a,b] 的⼀个⽹格剖分:a=x0

xi=a+ih(i=0,1,⋯,N), 步长

h=N.2. 对式中的⼆阶导数仍⽤数值微分公式b−ay(xi+1)−2y(xi)+y(xi−1)h2y′′(xi)=−12y(4)(ξi),h2xi−1<ξi

y(xi) 的近似值

yi 代替

y(xi),

pi=p(xi),qi=q(xi),fi=f(xi) 便得到差分⽅程组2yi+yi+1)+2h(yi+1−yi−1)+q1yi=fi,⎧h2(yi−1a−2⎨a1y0+2h(−3y0+4y1−y2)=αβ2⎩β1yN+2h(yN−2−4yN−1+3yN)=β1pii=1,2,⋯,N−1其中

qi=q(xi),pi=p(xi),fi=f(xi),i=1,2,⋯,N−1,yi 是

y(xi) 的近似值. 整理得⎧(2hα1−3α2)y0+4α2y1−α2y2=2hα⎨(2−hpi)yi−1−2(2−h2qi)yi+(2+hpi)yi+1=2h2fi,⎩β2yN−2−4β2yN−1+(3β2+2hβ1)yN=2hβi=1,2,⋯,N−1特别地, 若

α1=1,α2=0,β1=1,β2=0, 则原问题中的边界条件是第⼀类边值条件:

y(a)=α,y(b)=β; 此时⽅程组为⎧−2(2−h2q1)y1+(2+hp1)y2=2h2f1−(2−hp1)α,⎨(2−hpi)yi−1−2(2−h2q1)yi+(2+hpi)yi+1=2h2fi,i=2,3,⋯,N−2⎩(2−hpN−1)yN−2−2(2−h2qN−1)yN−1=2h2fN−1−(2+hpN−1)β以上⽅程组是三对⾓⽅程组，⽤解三对⾓⽅程组的追赶法求解差分⽅程组，便得边值问题的差分解 .3. 讨论差分⽅程组的解是否收敛到原微分⽅程的解，估计误差. 这⾥就不再详细介绍.数值算例考虑带有第三类边界条件的⼆阶常系数微分⽅程边值问题y′′(x)+2y′(x)−3y(x)=3x+1,⎧⎪y′(1)−y(1)=1,⎨x∈[1,2]′⎪⎩y(2)−y(2)=−4.1. 求出上述两点边值问题的解析解;2. 通过有限差分⽅法算出其数值解;(取h=5)并计算误差、绘图、输出1.0,1.2,1.4,1.6,1.8,2.0处的数值解.1问题⼆:有限差分⽅法算出其数值解及误差对于原问题,取步长h=0.2,⽤有限差分求其近似解,并将结果与精确解y(x)=-x-1进⾏⽐较.解:因为2−1N=h=5,p(x)=2,q(x)=−3,f(x)=3x+1α1=−1,α2=1,α=1β1=1,β2=1,β=−4代⼊上述差分⽅程组公式得到差分格式⎧(−2h−3)y0+4y1−y2=2h⎨(2−2h)yi−1−2(2+3h2)yi+(2+2h)yi+1=2h2fi,⎩yN−2−4yN−1+(3+2h)yN=−8hi=1,2,⋯,N−1化简得⎧(−2h−3)y0+4y1+(−1)y2=2h⎨(2−2h)yi−1+(−4−6h2)yi+(2+2h)yi+1=2h2fi,⎩yN−2+(−4)yN−1+(3+2h)yN=−8hi=1,2,⋯,N−1写成矩阵形式⎡y0⎤⎡2h⎤−2h−34−12⎡⎤y12⎢⎥⎢2h2f1⎥2−2h−4−6h2+2h⎢⎢⎥y2⎥⎢2hf2⎥2⎢⎥⎢⎥⎢⎥2−2h−4−6h2+2h⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⋱⋱⋱⋮⎥⎢⋮⎢⎥⎢⎥2⎢⎥⎢yN−2⎥⎢2h2fN−2⎥2−2h−4−6h2+2h⎢⎥⎢⎥⎢⎥⎢2−2h−4−6h22+2h⎥⎢yN−1⎥⎢2h2fN−1⎥⎣1−43+2h⎦⎣yN⎦=⎣−8h⎦令⎡y0⎤⎡2h⎤−2h−34−12⎡⎤y2hf112⎢⎥⎢⎥2+2h⎢y2⎥⎢2h2f2⎥⎢2−2h−4−6h⎥⎢⎥⎢⎥⎢⎥2−2h−4−6h22+2h⎢⎥⎢⎥⎢⎥⎢⎥⎢⋮⎥⎢⎥⋱⋱⋱⋮⎢⎥⎢⎥⎢⎥22⎢⎥⎢⎥⎢⎥2−2h−4−6h2+2h⎢⎥⎢yN−2⎥⎢2hfN−2⎥⎢⎢yN−1⎥⎢2h2fN−1⎥2−2h−4−6h22+2h⎥⎣1−43+2h⎦,Y=⎣yN⎦,b=⎣−8h⎦A=得到代数⽅程AY=b使⽤解三对⾓⽅程组的追赶法求解此差分⽅程组,得到差分解

Y.⼆阶线性常系数微分⽅程边值问题的差分法Python代码先以将区间划分为5份(h=0.2)为例，求出数值解import numpy as npimport as plt# 准备⼯作def initial(L, R, NS): x = ce(L, R, NS + 1) h = (R - L) / NS return x, h# 右端函数fdef f(x): return 3 * x + 1# 解⽅程def solve(NS): # 矩阵A A1 = ([-2 * h - 3] + [b] * (NS - 1) + [3 + 2 * h]) # print("A1",A1) A2 = ([a] * (NS - 1) + [-4]) # print("A2", A2) A3 = ([4] + [c] * (NS - 1)) # print("A3", A3) A4 = ([-1] + [0] * (NS - 2)) # print("A4", A4) A5 = ([0] * (NS - 2) + [1]) # print("A5", A5) A = (A1) + (A2, -1) + (A3, 1) + (A4, 2) + (A5, -2) # print("A", A) # 矩阵b br = 2 * h ** 2 * f(x) + ( [2 * h - 2 * h ** 2 * f(x[0])] + [0] * (NS - 1) + [-8 * h + 2 * h ** 2 * f(x[NS])]) uh = (A, br) return uhL = 1.0R = 2.0NS = 5x, h = initial(L, R, NS)a = 2 - 2 * hb = -4 - 6 * h ** 2c = 2 + 2 * huh = solve(NS)print("h=%f时数值解n" % h, uh)结果:h=0.200000时数值解 [-1.48151709 -1.56543883 -1.6245972 -1.6531625 -1.64418896 -1.58929215]是不是解出数值解就完事了呢？当然不是。我们可以将问题的差分格式解与问题的真解进⾏⽐较，以得到解的可靠性。通过数学计算我们得到问题的真解为

u(x)=−x−1，现⽤范数来衡量误差的⼤⼩：# 真解函数def ture_u(x): ture_u = -x - 1 return ture_ut_u = ture_u(x)print("h=%f时真解n" % h, t_u)# 误差范数def err(ture_u, uh): ee = max(abs(ture_u - uh)) e0 = (sum((ture_u - uh) ** 2) * h) return ee, e0ee, e0 = err(t_u, uh)print('L_∞(最⼤模)范数下的误差', ee)print('L_2(平⽅和)范数下的误差', e0)结果:h=0.200000时真解 [-2. -2.2 -2.4 -2.6 -2.8 -3. ]L_∞(最⼤模)范数下的误差 1.4164L_2(平⽅和)范数下的误差 1.8784接下来绘图⽐较

h=0.2 时数值解与真解的差距：# 绘图⽐较()(x, uh, label='Numerical solution')(x, t_u, label='Exact solution')("solution")()()结果:哎呀呀! 根据输出显⽰的误差,发现此时数值解的求解效果令⼈难以接受,

L∞(最⼤模)范数下的误差⾼达1.41,从图像也能看出数值解与解析解不仅相差甚远,甚⾄连曲线⾛势都不太⼀致.(数值解为⼀条曲线,解析解为直线)此处⾸先认为解析解或者查分格式计算错了,休息了⼀晚上起来重新推导了⼀遍后并未发现明显错误. 后来在玩GTA的时候忽然想起导师之前提到过步长会影响数值解稳定性···,于是重新编写程序,修改步长后发现确实如此.将区间划分为

N=1280 份, 即

h=0.000781 时.结果:h=0.000781时数值解 [-1.99798816 -1.99876784 -1.99954751 ... -2.99297729 -2.99375427 -2.99453125]h=0.000781时真解 [-2. -2.00078125 -2.0015625 ... -2.9984375 -2.99921875 -3. ]L_∞(最⼤模)范数下的误差 0.166322L_2(平⽅和)范数下的误差 0.5910903绘图⽐较

h=0.000781 时数值解与真解的差距：可以看到此时已经具有较⾼精度, 若还未达到需求精度, 可通过缩短步长h 继续提⾼精度.最后，我们还可以从数学的⾓度分析所采⽤的差分格式的⼀些性质。因为差分格式的误差为O(h2)， 所以理论上来说⽹格每加密⼀倍，与真解的误差⼤致会缩⼩到原来的4. 下讨论⽹格加密时的变化：# 误差与⽹格关系讨论N = [5, 10, 20, 40, 80, 160, 320, 640, 1280, 2560]print("⽹格密度加倍时误差变化情况")for i in range(len(N)): NS = N[i] x, h = initial(L, R, NS) a = 2 - 2 * h b = -4 - 6 * h ** 2 c = 2 + 2 * h uh = solve(NS) t_u = ture_u(x) ee, e0 = err(t_u, uh) print('步长h=%f' % h, end='t') print('最⼤模误差=%f' % ee)1结果：⽹格密度加倍时误差变化情况步长h=0.200000 最⼤模误差=1.410708步长h=0.100000 最⼤模误差=0.701417步长h=0.050000 最⼤模误差=0.350182步长h=0.025000 最⼤模误差=0.175023步长h=0.012500 最⼤模误差=0.087503步长h=0.006250 最⼤模误差=0.043750步长h=0.003125 最⼤模误差=0.021875步长h=0.001563 最⼤模误差=0.010938步长h=0.000781 最⼤模误差=0.005469步长h=0.000391 最⼤模误差=0.002734完整代码# 开发者: Leo 刘# 开发环境: macOs Big Sur# 开发时间: 2021/10/5 6:51 下午# 邮箱 : 517093978@# @Software: PyCharm# ----------------------------------------------------------------------------------------------------------import numpy as npimport as plt# 准备⼯作def initial(L, R, NS): x = ce(L, R, NS + 1) h = (R - L) / NS return x, h# 右端函数fdef f(x): return 3 * x + 1# 解⽅程def solve(NS): # 矩阵A A1 = ([-2 * h - 3] + [b] * (NS - 1) + [3 + 2 * h]) # print("A1",A1) A2 = ([a] * (NS - 1) + [-4]) # print("A2", A2) A3 = ([4] + [c] * (NS - 1)) # print("A3", A3) A4 = ([-1] + [0] * (NS - 2)) # print("A4", A4) A5 = ([0] * (NS - 2) + [1]) # print("A5", A5) A = (A1) + (A2, -1) + (A3, 1) + (A4, 2) + (A5, -2) # print("A", A) # 矩阵b br = 2 * h ** 2 * f(x) + ( [2 * h - 2 * h ** 2 * f(x[0])] + [0] * (NS - 1) + [-8 * h + 2 * h ** 2 * f(x[NS])]) uh = (A, br) return uhL = 1.0R = 2.0NS = 1280x, h = initial(L, R, NS)a = 2 - 2 * hb = -4 - 6 * h ** 2c = 2 + 2 * huh = solve(NS)print("h=%f时数值解n" % h, uh)# 真解函数def ture_u(x): ture_u = -x - 1 return ture_ut_u = ture_u(x)print("h=%f时真解n" % h, t_u)# 误差范数def err(ture_u, uh): ee = max(abs(ture_u - uh)) e0 = (sum((ture_u - uh) ** 2) * h) return ee, e0ee, e0 = err(t_u, uh)print('L_∞(最⼤模)范数下的误差', ee)print('L_2(平⽅和)范数下的误差', e0)# 绘图⽐较()(x, uh, label='Numerical solution')(x, t_u, label='Exact solution')("solution")()()# 误差与⽹格关系讨论N = [5, 10, 20, 40, 80, 160, 320, 640, 1280, 2560]print("⽹格密度加倍时误差变化情况")for i in range(len(N)): NS = N[i] x, h = initial(L, R, NS) a = 2 - 2 * h b = -4 - 6 * h ** 2 c = 2 + 2 * h uh = solve(NS) t_u = ture_u(x) ee, e0 = err(t_u, uh) print('步长h=%f' % h, end='t') print('最⼤模误差=%f' % ee)

数值算例来源: 《微分⽅程数值解》-