1.两点边值问题形式
一般的两点边值问题形式为:
−
u
x
x
=
f
(
x
)
-u_{xx}=f(x)
−uxx=f(x)
u
(
a
)
=
α
,
u
(
b
)
=
β
u(a)=\alpha,u(b)=\beta
u(a)=α,u(b)=β
下给出一个具体两点边值的问题用以分析:
−
u
x
x
=
16
π
2
s
i
n
(
4
π
x
)
-u_{xx}=16\pi^2sin(4\pi x)
−uxx=16π2sin(4πx)
u
(
0
)
=
0
,
u
(
1
)
=
0
u(0)=0,u(1)=0
u(0)=0,u(1)=0
其中问题的真解为:
u
(
x
)
=
s
i
n
(
4
π
x
)
u(x)=sin(4πx)
u(x)=sin(4πx)
2.求解思路
按照我们上一节中有限差分法步骤:
- 求解区域的划分
将区间[0,1]等距划分为n份,节点记为 x 0 , x 1 , x 2 . . . x n x_0,x_1,x_2...x_n x0,x1,x2...xn
- 选取差分格式
以二阶中心差分格式 x i + 1 + x i − 1 − 2 x i h 2 \frac{x_{i+1}+x_{i-1}-2x_{i}}{h^2} h2xi+1+xi−1−2xi 代替题目中的 u x x u_{xx} uxx
即可以得到 x i + 1 + x i − 1 − 2 x i h 2 = 16 π 2 s i n ( 4 π x ) \frac{x_{i+1}+x_{i-1}-2x_{i}}{h^2}=16\pi^2sin(4\pi x) h2xi+1+xi−1−2xi=16π2sin(4πx) (式2.1),其中 i =0,1,2…n - 边界处理
这里已经给出了具体的边值 u ( 0 ) = 0 , u ( 1 ) = 0 u(0)=0,u(1)=0 u(0)=0,u(1)=0 ,和方程是相容的,不需要额外的处理. - 得出代数方程组并求解
将上面的(式2.1)写成代数方程组的格式,即:
{ x 0 = 0 a x 0 + b x 1 + a x 2 = f ( x 1 ) a x 1 + b x 2 + a x 3 = f ( x 2 ) ⋮ a x n − 2 + b x n − 1 + a x n = f ( x n − 1 ) x n = 0 \begin{cases} x_0=0\\ ax_0+bx_1+ax_2=f(x_1)\\ ax_1+bx_2+ax_3=f(x_2)\\ \vdots\\ ax_{n-2}+bx_{n-1}+ax_n=f(x_{n-1})\\ x_n=0 \end{cases} ⎩⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎧x0=0ax0+bx1+ax2=f(x1)ax1+bx2+ax3=f(x2)⋮axn−2+bxn−1+axn=f(xn−1)xn=0
其中 a = − 1 h 2 , b = 2 h 2 a=-\frac{1}{h^2},b=\frac{2}{h^2} a=−h21,b=h22
将方程组左端转化为成矩阵形式,即为一个类三对角矩阵,将其记为A:
A = [ 1 0 0 0 ⋯ 0 a b a 0 ⋯ 0 0 a b a ⋯ 0 ⋮ ⋮ ⋮ ⋮ ⋱ ⋮ 0 0 ⋯ a b a 0 0 ⋯ 0 0 1 ] A=\begin{bmatrix} {1}&{0}&{0}&{0}&{\cdots}&{0}\\ {a}&{b}&{a}&{0}&{\cdots}&{0}\\ {0}&{a}&{b}&{a}&{\cdots}&{0}\\ {\vdots}&{\vdots}&{\vdots}&{\vdots}&{\ddots}&{\vdots}\\ {0}&{0}&{\cdots}&{a}&{b}&{a}\\ {0}&{0}&{\cdots}&{0}&{0}&{1}\\ \end{bmatrix} A=⎣⎢⎢⎢⎢⎢⎢⎢⎡1a0⋮000ba⋮000ab⋮⋯⋯00a⋮a0⋯⋯⋯⋱b0000⋮a1⎦⎥⎥⎥⎥⎥⎥⎥⎤
方程右端记为矩阵b:
b = [ f ( x 0 ) , f ( x 1 ) , f ( x 2 ) ⋯ f ( x n ) ] T b=\begin{bmatrix} {f(x_0)},{f(x_1)},{f(x_2)}&{\cdots}&{f(x_n)}\\ \end{bmatrix}^T b=[f(x0),f(x1),f(x2)⋯f(xn)]T
即求解 A x = b Ax=b Ax=b
3.求解代码与结果分析
先以将区间划分为20份为例,求出数值解
import numpy as np
#准备工作
def initial(L,R,NS):
x = np.linspace(L,R,NS+1)
h = (R-L) / NS
return x,h
#右端函数f
def f(x):
return 16*(np.pi)**2*np.sin(4*np.pi*x)
#解方程
def solve(NS):
#矩阵A
A1 = np.array([1]+[b]*(NS-1)+[1])
A2 = np.array([a]*(NS-1)+[0])
A3 = np.array([0]+[a]*(NS-1))
A = np.diag(A1)+np.diag(A2,-1)+np.diag(A3,1)
#矩阵b
br = f(x)
uh = np.linalg.solve(A,br)
return uh
L = 0.0
R = 1.0
NS = 20
x,h = initial(L,R,NS)
a = -1/h**2
b = 2/h**2
uh = solve(NS)
print("N=20时数值解\n",uh)
结果:
[Out]:
N=20时数值解
[ 1.09917379e-14 6.07510380e-01 9.82972443e-01 9.82972443e-01
6.07510380e-01 -1.13114821e-14 -6.07510380e-01 -9.82972443e-01
-9.82972443e-01 -6.07510380e-01 -3.28309853e-14 6.07510380e-01
9.82972443e-01 9.82972443e-01 6.07510380e-01 -5.45035987e-14
-6.07510380e-01 -9.82972443e-01 -9.82972443e-01 -6.07510380e-01
-7.73553884e-14]
是不是解出数值解就完事了呢?当然不是。我们可以将问题的差分格式解与问题的真解进行比较,以得到解的可靠性。通过数学计算我们得到问题的真解为 u ( x ) = s i n ( 4 π x ) u(x)=sin(4πx) u(x)=sin(4πx),现用范数来衡量误差的大小:
#真解函数
def ture_u(x):
ture_u = np.sin(4*np.pi*x)
return ture_u
t_u = ture_u(x)
print("N=20时真解\n",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('最大模范数下的误差',ee)
print('平方和范数下的误差',e0)
结果:
[Out]:
N=20时真解
[ 0.00000000e+00 5.87785252e-01 9.51056516e-01 9.51056516e-01
5.87785252e-01 1.22464680e-16 -5.87785252e-01 -9.51056516e-01
-9.51056516e-01 -5.87785252e-01 -2.44929360e-16 5.87785252e-01
9.51056516e-01 9.51056516e-01 5.87785252e-01 3.67394040e-16
-5.87785252e-01 -9.51056516e-01 -9.51056516e-01 -5.87785252e-01
-4.89858720e-16]
最大模范数下的误差 0.03191592653467212
平方和范数下的误差 0.023729365914438843
接下来绘图比较NS=20时与真解的差距:
#绘图比较
import matplotlib.pyplot as plt
plt.figure()
plt.plot(x,uh)
plt.plot(x,t_u)
plt.title("solution")
plt.show()
结果:
最后,我们还可以从数学的角度分析所采用的差分格式的一些性质。因为差分格式的误差为 O ( h 2 ) O(h^2) O(h2) ,所以理论上来说网格每加密一倍,与真解的误差大致会缩小到原来的 1 / 4 1/4 1/4。下讨论网格加密时的变化:
#误差与网格关系讨论
N = [10,20,40,80]
print("接下来是网格密度加倍时误差变化情况")
for i in range (4):
NS = N[i]
x,h = initial(L,R,NS)
a = -1/h**2
b = 2/h**2
uh = solve(NS)
t_u = ture_u(x)
ee,e0 = err(t_u,uh)
print(ee)
结果:
接下来是网格密度加倍时误差变化情况
0.13569108849388545
0.03191592653467212
0.00826541696629457
0.002058706764599627
由结果观察可知,当网格每加密一倍时误差确实减少为原来的 1 / 4 1/4 1/4。
至此,两点边值的相关问题已经分析完了,下一节我们将讨论均匀直棒热传导问题(抛物型问题)。博主在准备考研中,不过还是会尽快更新的~