【微分方程数值解】有限差分法（二）两点边值问题数值算例（附python代码）

1.两点边值问题形式

一般的两点边值问题形式为：
−uxx=f(x)-u_{xx}=f(x)−uxx=f(x)
u(a)=α,u(b)=βu(a)=\alpha,u(b)=\betau(a)=α,u(b)=β

下给出一个具体两点边值的问题用以分析：
−uxx=16π2sin(4πx)-u_{xx}=16\pi^2sin(4\pi x)−uxx=16π2sin(4πx)
u(0)=0,u(1)=0u(0)=0,u(1)=0u(0)=0,u(1)=0
其中问题的真解为：u(x)=sin(4πx)u(x)=sin(4πx)u(x)=sin(4πx)

2.求解思路

按照我们上一节中有限差分法步骤：

求解区域的划分
将区间[0,1]等距划分为n份，节点记为x0,x1,x2...xnx_0,x_1,x_2...x_nx0,x1,x2...xn
选取差分格式
以二阶中心差分格式 xi+1+xi−1−2xih2\frac{x_{i+1}+x_{i-1}-2x_{i}}{h^2}h2xi+1+xi−1−2xi 代替题目中的 uxxu_{xx}uxx
即可以得到 xi+1+xi−1−2xih2=16π2sin(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)=0u(0)=0,u(1)=0u(0)=0,u(1)=0 ，和方程是相容的，不需要额外的处理.
得出代数方程组并求解
将上面的（式2.1）写成代数方程组的格式，即：
{x0=0ax0+bx1+ax2=f(x1)ax1+bx2+ax3=f(x2)⋮axn−2+bxn−1+axn=f(xn−1)xn=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=−1h2,b=2h2a=-\frac{1}{h^2},b=\frac{2}{h^2}a=−h21,b=h22
将方程组左端转化为成矩阵形式，即为一个类三对角矩阵，将其记为A：
A=[1000⋯0aba0⋯00aba⋯0⋮⋮⋮⋮⋱⋮00⋯aba00⋯001]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(x0),f(x1),f(x2)⋯f(xn)]Tb=\begin{bmatrix} {f(x_0)},{f(x_1)},{f(x_2)}&{\cdots}&{f(x_n)}\\ \end{bmatrix}^Tb=[f(x0),f(x1),f(x2)⋯f(xn)]T
即求解 Ax=bAx=bAx=b

3.求解代码与结果分析

先以将区间划分为20份为例，求出数值解

import numpy as np
#准备工作
def initial(L,R,NS):x = np.linspace(L,R,NS+1)h = (R-L) / NSreturn x,h#右端函数f
def f(x):return 16*(np.pi)**2*np.sin(4*np.pi*x)
#解方程
def solve(NS):#矩阵AA1 = 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)#矩阵bbr = f(x)uh = np.linalg.solve(A,br)return uh
L = 0.0
R = 1.0
NS = 20x,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-016.07510380e-01 -1.13114821e-14 -6.07510380e-01 -9.82972443e-01-9.82972443e-01 -6.07510380e-01 -3.28309853e-14  6.07510380e-019.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)=sin(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-015.87785252e-01  1.22464680e-16 -5.87785252e-01 -9.51056516e-01-9.51056516e-01 -5.87785252e-01 -2.44929360e-16  5.87785252e-019.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(h2)O(h^2)O(h2) ，所以理论上来说网格每加密一倍，与真解的误差大致会缩小到原来的1/41/41/4。下讨论网格加密时的变化：

#误差与网格关系讨论
N = [10,20,40,80]
print("接下来是网格密度加倍时误差变化情况")
for i in range (4):NS = N[i]x,h = initial(L,R,NS)a = -1/h**2b = 2/h**2uh = solve(NS)t_u = ture_u(x)ee,e0 = err(t_u,uh)print(ee)

结果：

接下来是网格密度加倍时误差变化情况
0.13569108849388545
0.03191592653467212
0.00826541696629457
0.002058706764599627

由结果观察可知，当网格每加密一倍时误差确实减少为原来的1/41/41/4。

至此，两点边值的相关问题已经分析完了，下一节我们将讨论均匀直棒热传导问题（抛物型问题）。博主在准备考研中，不过还是会尽快更新的~