本文归纳常见常微分方程的解析解解法以及基于Python的微分方程数值解算例实现。
- 常微分方程的解析解(方法归纳)以及基于Python的微分方程数值解算例实现
- 常微分方程的解析解
- 可分离变量的微分方程（一阶）
- 一阶齐次（非齐次）线性微分方程（一阶）
- 二阶常系数微分方程（二阶）
- 高阶常系数微分方程（n阶）
- 算例
- 常微分方程的数值解
- 一般二阶线性常微分方程边值问题的差分法
- 数值算例
- 完整代码
考虑常微分方程的解析解法，我们一般可以将其归纳为如下几类：
- 可分离变量的微分方程（一阶）
- 一阶齐次（非齐次）线性微分方程（一阶）
- 二阶常系数微分方程（二阶）
- 高阶常系数微分方程（ n n n阶）
可分离变量的微分方程（一阶）
这类微分方程可以变形成如下形式：
f
(
x
)
d
x
=
g
(
y
)
d
y
f(x)dx=g(y)dy
f(x)dx=g(y)dy
两边同时积分即可解出函数，难点主要在于不定积分，是最简单的微分方程。
某些方程看似不可分离变量，但是经过换元之后，其实还是可分离变量的，不要被这种方程迷惑。
一阶齐次（非齐次）线性微分方程（一阶）
形如
d
y
d
x
+
P
(
x
)
y
=
Q
(
x
)
frac{dy}{dx}+P(x)y=Q(x)
dxdy+P(x)y=Q(x)
的方程叫做一阶线性微分方程，若
Q
(
x
)
Q(x)
Q(x)为0，则方程齐次，否则称为非齐次。
解法： (直接套公式)
y
(
x
)
=
e
−
∫
P
(
x
)
d
x
(
∫
e
∫
P
(
x
)
d
x
Q
(
x
)
d
x
+
C
)
y(x)=e^{-int{P(x)}dx}(int{e^{int{P(x)dx}}Q(x)}dx+C)
y(x)=e−∫P(x)dx(∫e∫P(x)dxQ(x)dx+C)
伯努利方程
形如
d
y
d
x
+
P
(
x
)
y
=
Q
(
x
)
y
n
,
n
∈
R
,
n
≠
1
frac{dy}{dx}+P(x)y=Q(x)y^{n},ninmathbb{R},nne1
dxdy+P(x)y=Q(x)yn,n∈R,n=1
的方程称为伯努利方程，这种方程可以通过以下步骤化为一阶线性微分方程：
y
−
n
d
y
d
x
+
P
(
x
)
y
1
−
n
=
Q
(
x
)
y^{-n}frac{dy}{dx}+P(x)y^{1-n}=Q(x)
y−ndxdy+P(x)y1−n=Q(x)
1
1
−
n
⋅
d
y
1
−
n
d
x
+
P
(
x
)
y
1
−
n
=
Q
(
x
)
frac{1}{1-n}·frac{dy^{1-n}}{dx}+P(x)y^{1-n}=Q(x)
1−n1⋅dxdy1−n+P(x)y1−n=Q(x)
令
y
1
−
n
=
u
y^{1-n}=u
y1−n=u, 方程两边同时乘以
1
−
n
1−n
1−n，得到
d
u
d
x
+
(
1
−
n
)
P
(
x
)
u
=
(
1
−
n
)
Q
(
x
)
frac{du}{dx}+(1-n)P(x)u=(1-n)Q(x)
dxdu+(1−n)P(x)u=(1−n)Q(x)
即
d
u
d
x
+
P
′
(
x
)
u
=
Q
′
(
x
)
frac{du}{dx}+P'(x)u=Q'(x)
dxdu+P′(x)u=Q′(x).
这就将伯努利方程归结为可以套公式的一阶线性微分方程。
二阶常系数微分方程（二阶）
形如
y
′
′
+
p
y
′
+
q
y
=
f
(
x
)
y''+py'+qy=f(x)
y′′+py′+qy=f(x)
的方程称为二阶常系数微分方程，若
f
(
x
)
≡
0
f(x)equiv0
f(x)≡0，则方程称为齐次的，反之称为非齐次的。以下默认方程是非齐次的。
求解此类方程分两步：
- 求出齐次通解
- 求出非齐次特解
原方程的解=齐次通解+非齐次特解
- 齐次通解的求法
首先假设
f
(
x
)
≡
0
f(x)equiv0
f(x)≡0.用特征方程法，写出对应的特征方程并且求解：
λ
2
+
p
λ
+
q
=
0
lambda^2+plambda+q=0
λ2+pλ+q=0
解的情况分为以下三种：
情况一：方程有两个不同的实数解
假设两个实数解分别是
λ
1
、
λ
2
lambda_1、lambda_2
λ1、λ2, 此时方程的通解是
Y
(
x
)
=
C
1
e
λ
1
x
+
C
2
e
λ
2
x
Y(x)=C_{1}e^{lambda_1x}+C_{2}e^{lambda_2x}
Y(x)=C1eλ1x+C2eλ2x
情况二：方程有一个二重解
假设该解等于
λ
lambda
λ，此时方程的通解是
Y
(
x
)
=
(
C
1
+
C
2
x
)
e
λ
x
Y(x)=(C_{1}+C_{2}x)e^{lambda x}
Y(x)=(C1+C2x)eλx
情况三：方程有一对共轭复解
假设这对解是
α
±
i
β
alphapm ibeta
α±iβ, 此时方程的通解是
Y
(
x
)
=
e
α
x
(
C
1
c
o
s
(
β
x
)
+
C
2
s
i
n
(
β
x
)
)
Y(x)=e^{alpha x}(C_{1}cos(beta x)+C_{2}sin(beta x))
Y(x)=eαx(C1cos(βx)+C2sin(βx))
- 非齐次特解的求法
对于 f ( x ) f(x) f(x) 和特征根的情况，对特解的情况做如下归纳：
-
f ( x ) = P m ( x ) f(x)=P_{m}(x) f(x)=Pm(x)，其中 P m ( x ) P_{m}(x) Pm(x) 表示 x x x 的最高次数为 m m m 的多项式。
若0不是方程特征解
则方程有特解 y ∗ = Q m ( x ) y^{*}=Q_{m}(x) y∗=Qm(x)
若0是方程的单特征解
则方程有特解 y ∗ = x Q m ( x ) y^{*}=xQ_{m}(x) y∗=xQm(x)
若0是方程的二重特征解
则方程有特解 y ∗ = x 2 Q m ( x ) y^{*}=x^{2}Q_{m}(x) y∗=x2Qm(x)
其中 Q m ( x ) = b 0 + b 1 x + … + b m x m Q_{m}(x)=b_{0}+b_{1}x+…+b_{m}x^{m} Qm(x)=b0+b1x+…+bmxm, b i ( i = 0 , 1 , … , m ) b_{i}(i = 0,1,…,m) bi(i=0,1,…,m)是需要带回原方程来确定的系数。
-
f ( x ) = e α x P m ( x ) f(x)=e^{alpha x}P_{m}(x) f(x)=eαxPm(x)
若 α alpha α 不是方程特征解
则方程有特解 y ∗ = e α x Q m ( x ) y^{*}=e^{alpha x}Q_{m}(x) y∗=eαxQm(x)
若 α alpha α 是方程的单特征解
则方程有特解 y ∗ = x e α x Q m ( x ) y^{*}=xe^{alpha x}Q_{m}(x) y∗=xeαxQm(x)
若 α alpha α 是方程的二重特征解
则方程有特解 y ∗ = x 2 e α x Q m ( x ) y^{*}=x^2e^{alpha x}Q_{m}(x) y∗=x2eαxQm(x)
-
$f(x)=e^{alpha x}(a_{1}cos(beta x)+a_{2}sin(beta x)) $
若 α ± i β alphapm ibeta α±iβ 不是特征解
则方程有特解 y ∗ = e α x ( A 1 c o s ( β x ) + A 2 s i n ( β x ) ) y^{*}=e^{alpha x}(A_{1}cos(beta x)+A_{2}sin(beta x)) y∗=eαx(A1cos(βx)+A2sin(βx))
若 α ± i β alphapm ibeta α±iβ 是特征解
则方程有特解 y ∗ = x e α x ( A 1 c o s ( β x ) + A 2 s i n ( β x ) ) y^{*}=xe^{alpha x}(A_{1}cos(beta x)+A_{2}sin(beta x)) y∗=xeαx(A1cos(βx)+A2sin(βx))
其中 A 1 、 A 2 A_{1}、A_{2} A1、A2 是需要带回原方程来确定的系数。
高阶常系数微分方程（n阶）
形如
y
(
n
)
+
p
1
y
(
n
−
1
)
+
.
.
.
+
p
n
−
1
y
′
+
p
n
y
=
f
(
x
)
y^{(n)}+p_{1}y^{(n-1)}+...+p_{n-1}y'+p_{n}y=f(x)
y(n)+p1y(n−1)+...+pn−1y′+pny=f(x)
的方程叫做高阶常系数微分方程，若
f
(
x
)
≡
0
f(x)equiv0
f(x)≡0，则方程是齐次的，否则是非齐次的。下面默认方程是非齐次的。
求解此类方程分两步：
- 求出齐次通解
- 求出非齐次特解
原方程的解=齐次通解+非齐次特解
- 齐次通解的求法(参考二阶常系数微分方程解法)
- 非齐次特解的求法(参考二阶常系数微分方程解法)
算例
考虑带有第三类边界条件的二阶常系数微分方程边值问题
{
y
′
′
(
x
)
+
2
y
′
(
x
)
−
3
y
(
x
)
=
3
x
+
1
,
x
∈
[
1
,
2
]
y
′
(
1
)
−
y
(
1
)
=
1
,
y
′
(
2
)
−
y
(
2
)
=
−
4.
left{ begin{aligned} y''(x)+2y'(x)-3y(x) & = 3x+1, & xin [1,2]\ y'(1)-y(1) & = 1,\ y'(2)-y(2) & = -4. end{aligned} right.
⎩⎪⎨⎪⎧y′′(x)+2y′(x)−3y(x)y′(1)−y(1)y′(2)−y(2)=3x+1,=1,=−4.x∈[1,2]
- 求出上述两点边值问题的解析解;
- 通过有限差分方法算出其数值解;(取 h = 1 5 h=frac{1}{5} h=51)并计算误差、绘图、输出 1.0 , 1.2 , 1.4 , 1.6 , 1.8 , 2.0 1.0,1.2,1.4,1.6,1.8,2.0 1.0,1.2,1.4,1.6,1.8,2.0处的数值解.
问题一:两点边值问题的解析解
由于此方程是非齐次的,故求解此类方程分两步：
- 求出齐次通解
- 求出非齐次特解
原方程的解=齐次通解+非齐次特解
- 齐次通解的求法
首先假设
y
′
′
(
x
)
+
2
y
′
(
x
)
−
3
y
(
x
)
≡
0
y''(x)+2y'(x)-3y(x)equiv0
y′′(x)+2y′(x)−3y(x)≡0. 用特征方程法，写出对应的特征方程
λ
2
+
2
λ
−
3
=
0
lambda^2+2lambda-3=0
λ2+2λ−3=0
求解得到两个不同的实数特征根:
λ
1
=
1
,
λ
2
=
−
3
lambda_1=1,lambda_2=-3
λ1=1,λ2=−3.
此时方程的齐次通解是
y
(
x
)
=
C
1
e
λ
1
x
+
C
2
e
λ
2
x
y(x)=C_{1}e^{lambda_1x}+C_{2}e^{lambda_2x}
y(x)=C1eλ1x+C2eλ2x
- 非齐次特解的求法
由于
λ
1
≠
0
,
λ
2
≠
0
lambda_1ne0,lambda_2ne0
λ1=0,λ2=0. 所以非齐次特解形式为
y
∗
=
a
x
+
b
y^{*}=ax+b
y∗=ax+b
将上式代入控制方程有
2
a
−
3
a
x
−
3
b
=
3
x
+
1
2a-3ax-3b=3x+1
2a−3ax−3b=3x+1
求解得:
a
=
b
=
−
1
a=b=-1
a=b=−1, 即非齐次特解为
y
∗
=
−
x
−
1
y^{*}=-x-1
y∗=−x−1.
原方程的解=齐次通解+非齐次特解
于是,原方程的全解为
y
(
x
)
=
C
1
e
x
+
C
2
e
−
3
x
−
x
−
1
y(x)=C_{1}e^{x}+C_{2}e^{-3x}-x-1
y(x)=C1ex+C2e−3x−x−1
因为该问题给出的是第三类边界条件,故需要求解的导函数
y
′
(
x
)
=
C
1
e
x
−
3
C
2
e
−
3
x
−
1
y'(x)=C_{1}e^{x}-3C_{2}e^{-3x}-1
y′(x)=C1ex−3C2e−3x−1
且有
{
y
(
1
)
=
C
1
e
+
C
2
e
−
3
−
2
y
′
(
1
)
=
C
1
e
−
3
C
2
e
−
3
−
1
y
(
2
)
=
C
1
e
2
+
C
2
e
−
6
−
3
y
′
(
2
)
=
C
1
e
2
−
3
C
2
e
−
6
−
1
left{ begin{aligned} y(1)&=C_1e+C_2e^{-3}-2\ y'(1)&=C_1e-3C_2e^{-3}-1\ y(2)&=C_1e^2+C_2e^{-6}-3\ y'(2)&=C_1e^2-3C_2e^{-6}-1 end{aligned} right.
⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧y(1)y′(1)y(2)y′(2)=C1e+C2e−3−2=C1e−3C2e−3−1=C1e2+C2e−6−3=C1e2−3C2e−6−1
将以上各式代入边界条件
{
y
′
(
1
)
−
y
(
1
)
=
−
4
C
2
e
−
3
+
1
=
1
y
′
(
2
)
+
y
(
2
)
=
2
C
1
e
2
−
2
C
2
e
−
6
−
4
=
−
4
left{ begin{aligned} y'(1)-y(1)&=-4C_2e^{-3}+1=1\ y'(2)+y(2)&=2C_1e^2-2C_2e^{-6}-4=-4 end{aligned} right.
{y′(1)−y(1)y′(2)+y(2)=−4C2e−3+1=1=2C1e2−2C2e−6−4=−4
解此方程组可得:
C
1
=
C
2
=
0
C_1=C_2=0
C1=C2=0.
综上所述,原两点边值问题的解为
y
=
−
x
−
1
y=-x-1
y=−x−1
对一般的二阶微分方程边值问题
{
y
′
′
(
x
)
+
p
(
x
)
y
′
(
x
)
+
q
(
x
)
y
(
x
)
=
f
(
x
)
,
a
<
x
<
b
a
1
y
(
a
)
+
a
2
y
′
(
a
)
=
a
β
1
y
(
b
)
+
β
2
y
′
(
b
)
=
β
left{begin{array}{l} y^{prime prime}(x)+p(x) y^{prime}(x)+q(x) y(x)=f(x), quad a
为求解的近似值, 类似于前面的做法,
-
把区间 I = [ a , b ] N I=[a, b] N I=[a,b]N 等分, 即得到区间 I = [ a , b ] I=[a, b] I=[a,b] 的一个网格剖分:
a = x 0 < x 1 < ⋯ < x N − 1 < x N = b a=x_{0} -
对式中的二阶导数仍用数值微分公式
y ′ ′ ( x i ) = y ( x i + 1 ) − 2 y ( x i ) + y ( x i − 1 ) h 2 − h 2 12 y ( 4 ) ( ξ i ) , x i − 1 < ξ i < x i y^{prime prime}left(x_{i}right)=frac{yleft(x_{i+1}right)-2 yleft(x_{i}right)+yleft(x_{i-1}right)}{h^{2}}-frac{h^{2}}{12} y^{(4)}left(xi_{i}right), quad x_{i-1}
{ y ′ ( x i ) = y ( x i + 1 ) − y ( x i − 1 ) 2 h − h 2 6 y ′ ′ ′ ( ξ i ) , x i − 1 < ξ i < x i , i = 1 , 2 , ⋯ , N − 1 y ′ ( x 0 ) = − 3 y ( x 0 ) + 4 y ( x 1 ) − y ( x 2 ) 2 h + h 2 3 y ′ ′ ′ ( ξ 0 ) , x 0 < ξ 0 < x 2 , y ′ ( x N ) = y ( x N − 2 ) − 4 y ( x N − 1 ) + 3 y ( x N ) 2 h + h 2 3 y ′ ′ ′ ( ξ N ) , x N − 2 < ξ N < x N left{begin{array}{l} y^{prime}left(x_{i}right)=frac{yleft(x_{i+1}right)-yleft(x_{i-1}right)}{2 h}-frac{h^{2}}{6} y^{prime prime prime}left(xi_{i}right), quad x_{i-1}
{ 1 h 2 ( y i − 1 − 2 y i + y i + 1 ) + p i 2 h ( y i + 1 − y i − 1 ) + q 1 y i = f i , i = 1 , 2 , ⋯ , N − 1 a 1 y 0 + a 2 2 h ( − 3 y 0 + 4 y 1 − y 2 ) = α β 1 y N + β 2 2 h ( y N − 2 − 4 y N − 1 + 3 y N ) = β left{begin{array}{l} frac{1}{h^{2}}left(y_{i-1}-2 y_{i}+y_{i+1}right)+frac{p_{i}}{2 h}left(y_{i+1}-y_{i-1}right)+q_{1} y_{i}=f_{i}, quad i=1,2, cdots, N-1 \ a_{1} y_{0}+frac{a_{2}}{2 h}left(-3 y_{0}+4 y_{1}-y_{2}right)=alpha \ beta_{1} y_{N}+frac{beta_{2}}{2 h}left(y_{N-2}-4 y_{N-1}+3 y_{N}right)=beta end{array}right. ⎩⎨⎧h21(yi−1−2yi+yi+1)+2hpi(yi+1−yi−1)+q1yi=fi,i=1,2,⋯,N−1a1y0+2ha2(−3y0+4y1−y2)=αβ1yN+2hβ2(yN−2−4yN−1+3yN)=β
其中 q i = q ( x i ) , p i = p ( x i ) , f i = f ( x i ) , i = 1 , 2 , ⋯ , N − 1 , y i q_{i}=qleft(x_{i}right), p_{i}=pleft(x_{i}right), f_{i}=fleft(x_{i}right), i=1,2, cdots, N-1, y_{i} qi=q(xi),pi=p(xi),fi=f(xi),i=1,2,⋯,N−1,yi 是 y ( x i ) yleft(x_{i}right) y(xi) 的近似值. 整理得
{ ( 2 h α 1 − 3 α 2 ) y 0 + 4 α 2 y 1 − α 2 y 2 = 2 h α ( 2 − h p i ) y i − 1 − 2 ( 2 − h 2 q i ) y i + ( 2 + h p i ) y i + 1 = 2 h 2 f i , i = 1 , 2 , ⋯ , N − 1 β 2 y N − 2 − 4 β 2 y N − 1 + ( 3 β 2 + 2 h β 1 ) y N = 2 h β left{begin{array}{l}left(2 h alpha_{1}-3 alpha_{2}right) y_{0}+4 alpha_{2} y_{1}-alpha_{2} y_{2}=2 h alpha \ left(2-h p_{i}right) y_{i-1}-2left(2-h^{2} q_{i}right) y_{i}+left(2+h p_{i}right) y_{i+1}=2 h^{2} f_{i}, quad i=1,2, cdots, N-1 \ beta_{2} y_{N-2}-4 beta_{2} y_{N-1}+left(3 beta_{2}+2 h beta_{1}right) y_{N}=2 h betaend{array}right. ⎩⎨⎧(2hα1−3α2)y0+4α2y1−α2y2=2hα(2−hpi)yi−1−2(2−h2qi)yi+(2+hpi)yi+1=2h2fi,i=1,2,⋯,N−1β2yN−2−4β2yN−1+(3β2+2hβ1)yN=2hβ
特别地, 若 α 1 = 1 , α 2 = 0 , β 1 = 1 , β 2 = 0 alpha_{1}=1, alpha_{2}=0, beta_{1}=1, beta_{2}=0 α1=1,α2=0,β1=1,β2=0, 则原问题中的边界条件是第一类边值条件: y ( a ) = α , y ( b ) = β y(a)=alpha, y(b)=beta y(a)=α,y(b)=β; 此时方程组为
{ − 2 ( 2 − h 2 q 1 ) y 1 + ( 2 + h p 1 ) y 2 = 2 h 2 f 1 − ( 2 − h p 1 ) α , ( 2 − h p i ) y i − 1 − 2 ( 2 − h 2 q 1 ) y i + ( 2 + h p i ) y i + 1 = 2 h 2 f i , i = 2 , 3 , ⋯ , N − 2 ( 2 − h p N − 1 ) y N − 2 − 2 ( 2 − h 2 q N − 1 ) y N − 1 = 2 h 2 f N − 1 − ( 2 + h p N − 1 ) β left{begin{array}{l} -2left(2-h^{2} q_{1}right) y_{1}+left(2+h p_{1}right) y_{2}=2 h^{2} f_{1}-left(2-h p_{1}right) alpha, \ left(2-h p_{i}right) y_{i-1}-2left(2-h^{2} q_{1}right) y_{i}+left(2+h p_{i}right) y_{i+1}=2 h^{2} f_{i}, quad i=2,3, cdots, N-2 \ left(2-h p_{N-1}right) y_{N-2}-2left(2-h^{2} q_{N-1}right) y_{N-1}=2 h^{2} f_{N-1}-left(2+h p_{N-1}right) beta end{array}right. ⎩⎨⎧−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)β
以上方程组是三对角方程组，用解三对角方程组的追赶法求解差分方程组，便得边值问题的差分解 . -
讨论差分方程组的解是否收敛到原微分方程的解，估计误差. 这里就不再详细介绍.
数值算例
考虑带有第三类边界条件的二阶常系数微分方程边值问题
{
y
′
′
(
x
)
+
2
y
′
(
x
)
−
3
y
(
x
)
=
3
x
+
1
,
x
∈
[
1
,
2
]
y
′
(
1
)
−
y
(
1
)
=
1
,
y
′
(
2
)
−
y
(
2
)
=
−
4.
left{ begin{aligned} y''(x)+2y'(x)-3y(x) & = 3x+1, & xin [1,2]\ y'(1)-y(1) & = 1,\ y'(2)-y(2) & = -4. end{aligned} right.
⎩⎪⎨⎪⎧y′′(x)+2y′(x)−3y(x)y′(1)−y(1)y′(2)−y(2)=3x+1,=1,=−4.x∈[1,2]
- 求出上述两点边值问题的解析解;
- 通过有限差分方法算出其数值解;(取 h = 1 5 h=frac{1}{5} h=51)并计算误差、绘图、输出 1.0 , 1.2 , 1.4 , 1.6 , 1.8 , 2.0 1.0,1.2,1.4,1.6,1.8,2.0 1.0,1.2,1.4,1.6,1.8,2.0处的数值解.
问题二:有限差分方法算出其数值解及误差
对于原问题,取步长h=0.2,用有限差分求其近似解,并将结果与精确解y(x)=-x-1进行比较.
解:
因为
N = 2 − 1 h = 5 , p ( x ) = 2 , q ( x ) = − 3 , f ( x ) = 3 x + 1 N=frac{2-1}{h}=5,p(x)=2,q(x)=-3,f(x)=3x+1 N=h2−1=5,p(x)=2,q(x)=−3,f(x)=3x+1
α 1 = − 1 , α 2 = 1 , α = 1 alpha_1=-1,alpha_2=1,alpha=1 α1=−1,α2=1,α=1
β 1 = 1 , β 2 = 1 , β = − 4 beta_1=1,beta_2=1,beta=-4 β1=1,β2=1,β=−4
代入上述差分方程组公式得到差分格式
{ ( − 2 h − 3 ) y 0 + 4 y 1 − y 2 = 2 h ( 2 − 2 h ) y i − 1 − 2 ( 2 + 3 h 2 ) y i + ( 2 + 2 h ) y i + 1 = 2 h 2 f i , i = 1 , 2 , ⋯ , N − 1 y N − 2 − 4 y N − 1 + ( 3 + 2 h ) y N = − 8 h left{begin{array}{l} left(-2h-3right) y_{0}+4 y_{1}- y_{2}=2h \ left(2-2hright) y_{i-1}-2left(2+3h^{2} right) y_{i}+left(2+2hright) y_{i+1}=2h^2f_i, quad i=1,2, cdots, N-1 \ y_{N-2}-4 y_{N-1}+left(3 + 2hright) y_{N}=-8hend{array}right. ⎩⎨⎧(−2h−3)y0+4y1−y2=2h(2−2h)yi−1−2(2+3h2)yi+(2+2h)yi+1=2h2fi,i=1,2,⋯,N−1yN−2−4yN−1+(3+2h)yN=−8h
化简得
{ ( − 2 h − 3 ) y 0 + 4 y 1 + ( − 1 ) y 2 = 2 h ( 2 − 2 h ) y i − 1 + ( − 4 − 6 h 2 ) y i + ( 2 + 2 h ) y i + 1 = 2 h 2 f i , i = 1 , 2 , ⋯ , N − 1 y N − 2 + ( − 4 ) y N − 1 + ( 3 + 2 h ) y N = − 8 h left{begin{array}{l} left(-2h-3right) y_{0}+4 y_{1}+(-1) y_{2}=2h \ left(2-2hright) y_{i-1}+left(-4-6h^{2} right) y_{i}+left(2+2hright) y_{i+1}=2h^2f_i, quad i=1,2, cdots, N-1 \ y_{N-2}+(-4) y_{N-1}+left(3 + 2hright) y_{N}=-8hend{array}right. ⎩⎨⎧(−2h−3)y0+4y1+(−1)y2=2h(2−2h)yi−1+(−4−6h2)yi+(2+2h)yi+1=2h2fi,i=1,2,⋯,N−1yN−2+(−4)yN−1+(3+2h)yN=−8h
写成矩阵形式
[ − 2 h − 3 4 − 1 2 − 2 h − 4 − 6 h 2 2 + 2 h 2 − 2 h − 4 − 6 h 2 2 + 2 h ⋱ ⋱ ⋱ 2 − 2 h − 4 − 6 h 2 2 + 2 h 2 − 2 h − 4 − 6 h 2 2 + 2 h 1 − 4 3 + 2 h ] [ y 0 y 1 y 2 ⋮ y N − 2 y N − 1 y N ] = [ 2 h 2 h 2 f 1 2 h 2 f 2 ⋮ 2 h 2 f N − 2 2 h 2 f N − 1 − 8 h ] left[ begin{matrix} -2h-3 & 4 & -1 & & & & \ 2-2h & -4-6h^2 & 2+2h & & & & \ & 2-2h & -4-6h^2 & 2+2h & & & \ & & ddots & ddots & ddots & & \ & & & 2-2h & -4-6h^2 & 2+2h \ & & & & 2-2h & -4-6h^2 & 2+2h \ & & & & 1 & -4 & 3+2h end{matrix} right] left[ begin{matrix} y_0\ y_1 \ y_2 \ vdots \ y_{N-2} \ y_{N-1} \ y_{N} end{matrix} right]=left[ begin{matrix} 2h \ 2h^2f_1 \ 2h^2f_2 \ vdots \ 2h^2f_{N-2} \ 2h^2f_{N-1} \ -8h end{matrix} right] ⎣⎢⎢⎢⎢⎢⎢⎢⎢⎡−2h−32−2h4−4−6h22−2h−12+2h−4−6h2⋱2+2h⋱2−2h⋱−4−6h22−2h12+2h−4−6h2−42+2h3+2h⎦⎥⎥⎥⎥⎥⎥⎥⎥⎤⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡y0y1y2⋮yN−2yN−1yN⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤=⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡2h2h2f12h2f2⋮2h2fN−22h2fN−1−8h⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤
令
A = [ − 2 h − 3 4 − 1 2 − 2 h − 4 − 6 h 2 2 + 2 h 2 − 2 h − 4 − 6 h 2 2 + 2 h ⋱ ⋱ ⋱ 2 − 2 h − 4 − 6 h 2 2 + 2 h 2 − 2 h − 4 − 6 h 2 2 + 2 h 1 − 4 3 + 2 h ] , Y = [ y 0 y 1 y 2 ⋮ y N − 2 y N − 1 y N ] , b = [ 2 h 2 h 2 f 1 2 h 2 f 2 ⋮ 2 h 2 f N − 2 2 h 2 f N − 1 − 8 h ] A=left[ begin{matrix} -2h-3 & 4 & -1 & & & & \ 2-2h & -4-6h^2 & 2+2h & & & & \ & 2-2h & -4-6h^2 & 2+2h & & & \ & & ddots & ddots & ddots & & \ & & & 2-2h & -4-6h^2 & 2+2h \ & & & & 2-2h & -4-6h^2 & 2+2h \ & & & & 1 & -4 & 3+2h end{matrix} right] , Y=left[ begin{matrix} y_0\ y_1 \ y_2 \ vdots \ y_{N-2} \ y_{N-1} \ y_{N} end{matrix}right] , b=left[ begin{matrix} 2h \ 2h^2f_1 \ 2h^2f_2 \ vdots \ 2h^2f_{N-2} \ 2h^2f_{N-1} \ -8h end{matrix} right] A=⎣⎢⎢⎢⎢⎢⎢⎢⎢⎡−2h−32−2h4−4−6h22−2h−12+2h−4−6h2⋱2+2h⋱2−2h⋱−4−6h22−2h12+2h−4−6h2−42+2h3+2h⎦⎥⎥⎥⎥⎥⎥⎥⎥⎤,Y=⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡y0y1y2⋮yN−2yN−1yN⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤,b=⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡2h2h2f12h2f2⋮2h2fN−22h2fN−1−8h⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤
得到代数方程
A Y = b AY=b AY=b
使用解三对角方程组的追赶法求解此差分方程组,得到差分解 Y Y Y.
二阶线性常系数微分方程边值问题的差分法Python代码
先以将区间划分为5份( h = 0.2 h=0.2 h=0.2)为例，求出数值解
import numpy as np
import matplotlib.pyplot as plt
# 准备工作
def initial(L, R, NS):
x = np.linspace(L, R, NS + 1)
h = (R - L) / NS
return x, h
# 右端函数f
def f(x):
return 3 * x + 1
# 解方程
def solve(NS):
# 矩阵A
A1 = np.array([-2 * h - 3] + [b] * (NS - 1) + [3 + 2 * h])
# print("A1",A1)
A2 = np.array([a] * (NS - 1) + [-4])
# print("A2", A2)
A3 = np.array([4] + [c] * (NS - 1))
# print("A3", A3)
A4 = np.array([-1] + [0] * (NS - 2))
# print("A4", A4)
A5 = np.array([0] * (NS - 2) + [1])
# print("A5", A5)
A = np.diag(A1) + np.diag(A2, -1) + np.diag(A3, 1) + np.diag(A4, 2) + np.diag(A5, -2)
# print("A", A)
# 矩阵b
br = 2 * h ** 2 * f(x) + np.array(
[2 * h - 2 * h ** 2 * f(x[0])] + [0] * (NS - 1) + [-8 * h + 2 * h ** 2 * f(x[NS])])
uh = np.linalg.solve(A, br)
return uh
L = 1.0
R = 2.0
NS = 5
x, h = initial(L, R, NS)
a = 2 - 2 * h
b = -4 - 6 * h ** 2
c = 2 + 2 * h
uh = 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 u ( x ) = -x-1 u(x)=−x−1，现用范数来衡量误差的大小：
# 真解函数
def ture_u(x):
ture_u = -x - 1
return ture_u
t_u = ture_u(x)
print("h=%f时真解n" % h, t_u)
# 误差范数
def err(ture_u, uh):
ee = max(abs(ture_u - uh))
e0 = np.sqrt(sum((ture_u - uh) ** 2) * h)
return ee, e0
ee, 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.4107078471946164 L_2(平方和)范数下的误差 1.0483548020308784
接下来绘图比较 h = 0.2 h=0.2 h=0.2 时数值解与真解的差距：
# 绘图比较
plt.figure()
plt.plot(x, uh, label='Numerical solution')
plt.plot(x, t_u, label='Exact solution')
plt.title("solution")
plt.legend()
plt.show()
结果:
哎呀呀! 根据输出显示的误差,发现此时数值解的求解效果令人难以接受,
L
∞
L_∞
L∞(最大模)范数下的误差高达1.41,从图像也能看出数值解与解析解不仅相差甚远,甚至连曲线走势都不太一致.(数值解为一条曲线,解析解为直线)
此处首先认为解析解或者查分格式计算错了,休息了一晚上起来重新推导了一遍后并未发现明显错误. 后来在玩GTA的时候忽然想起导师之前提到过步长会影响数值解稳定性···,于是重新编写程序,修改步长后发现确实如此.
将区间划分为 N = 1280 N=1280 N=1280 份, 即 h = 0.000781 h=0.000781 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.005468750663166322 L_2(平方和)范数下的误差 0.0035976563635910903
绘图比较
h
=
0.000781
h=0.000781
h=0.000781 时数值解与真解的差距：
可以看到此时已经具有较高精度, 若还未达到需求精度, 可通过缩短步长
h
h
h 继续提高精度.
最后，我们还可以从数学的角度分析所采用的差分格式的一些性质。因为差分格式的误差为 O ( h 2 ) O(h^2) O(h2)， 所以理论上来说网格每加密一倍，与真解的误差大致会缩小到原来的 1 4 frac{1}{4} 41. 下讨论网格加密时的变化：
# 误差与网格关系讨论
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)
结果：
网格密度加倍时误差变化情况 步长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@qq.com
# @Software: PyCharm
# ----------------------------------------------------------------------------------------------------------
import numpy as np
import matplotlib.pyplot as plt
# 准备工作
def initial(L, R, NS):
x = np.linspace(L, R, NS + 1)
h = (R - L) / NS
return x, h
# 右端函数f
def f(x):
return 3 * x + 1
# 解方程
def solve(NS):
# 矩阵A
A1 = np.array([-2 * h - 3] + [b] * (NS - 1) + [3 + 2 * h])
# print("A1",A1)
A2 = np.array([a] * (NS - 1) + [-4])
# print("A2", A2)
A3 = np.array([4] + [c] * (NS - 1))
# print("A3", A3)
A4 = np.array([-1] + [0] * (NS - 2))
# print("A4", A4)
A5 = np.array([0] * (NS - 2) + [1])
# print("A5", A5)
A = np.diag(A1) + np.diag(A2, -1) + np.diag(A3, 1) + np.diag(A4, 2) + np.diag(A5, -2)
# print("A", A)
# 矩阵b
br = 2 * h ** 2 * f(x) + np.array(
[2 * h - 2 * h ** 2 * f(x[0])] + [0] * (NS - 1) + [-8 * h + 2 * h ** 2 * f(x[NS])])
uh = np.linalg.solve(A, br)
return uh
L = 1.0
R = 2.0
NS = 1280
x, h = initial(L, R, NS)
a = 2 - 2 * h
b = -4 - 6 * h ** 2
c = 2 + 2 * h
uh = solve(NS)
print("h=%f时数值解n" % h, uh)
# 真解函数
def ture_u(x):
ture_u = -x - 1
return ture_u
t_u = ture_u(x)
print("h=%f时真解n" % h, t_u)
# 误差范数
def err(ture_u, uh):
ee = max(abs(ture_u - uh))
e0 = np.sqrt(sum((ture_u - uh) ** 2) * h)
return ee, e0
ee, e0 = err(t_u, uh)
print('L_∞(最大模)范数下的误差', ee)
print('L_2(平方和)范数下的误差', e0)
# 绘图比较
plt.figure()
plt.plot(x, uh, label='Numerical solution')
plt.plot(x, t_u, label='Exact solution')
plt.title("solution")
plt.legend()
plt.show()
# 误差与网格关系讨论
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)
数值算例来源: 《微分方程数值解》-M.Ran
