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

1.两点边值问题形式

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

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

2.求解思路

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

1. 求解区域的划分
   将区间[0,1]等距划分为n份，节点记为$x_0,x_1,x_2...x_n$
   
2. 选取差分格式
   以二阶中心差分格式 $\frac{x_{i+1}+x_{i-1}-2x_i}{h^2}$ 代替题目中的 $u_{xx}$
   即可以得到 $\frac{x_{i+1}+x_{i-1}-2x_i}{h^2}=16{\pi }^{2}sin(4\pi x)$ （式2.1），其中 i =0,1,2…n
3. 边界处理
   这里已经给出了具体的边值 $u(0)=0,u(1)=0$ ，和方程是相容的，不需要额外的处理.
4. 得出代数方程组并求解
   将上面的（式2.1）写成代数方程组的格式，即：
   $$\{\begin{array}{l}{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\end{array}$$
   其中$a=-\frac{1}{h^2},b=\frac{2}{h^2}$
   将方程组左端转化为成矩阵形式，即为一个类三对角矩阵，将其记为A：
   $$A=\left[\begin{array}{cccccc}1& 0& 0& 0& \cdots & 0\\ a& b& a& 0& \cdots & 0\\ 0& a& b& a& \cdots & 0\\ ⋮& ⋮& ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & a& b& a\\ 0& 0& \cdots & 0& 0& 1\end{array}\right]$$
   方程右端记为矩阵b:
   $$b={\left[\begin{array}{ccc}f(x_0),f(x_1),f(x_2)& \cdots & f(x_n)\end{array}\right]}^T$$
   即求解 $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)=sin(4\pi 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)$ ，所以理论上来说网格每加密一倍，与真解的误差大致会缩小到原来的$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$。

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