二维Poisson方程五点差分格式及简单求解方法Python实现

二维Poisson方程简介

给出 二维 $Poisson$ 方程 $Dirichlet$ 边值问题:
{ − Δ u = f ( x , y ) ( x , y ) ∈ Ω u = φ ( x , y ) ( x , y ) ∈ Γ

$$\begin{cases} - \Delta u = f(x, y) \quad (x, y) \in \Omega \\ u = \varphi(x, y) \quad (x, y) \in \Gamma \end{cases}$$ {−Δu=f(x,y)(x,y)∈Ωu=φ(x,y)(x,y)∈Γ​
为简单起见， 假设 $$\Omega =[a,b]\times [c,d]$$

五点差分格式建立

将区间 $[a,b]$ 作 $m$ 等分, 记 $$h_1=(b-a)/m,x_i=a+i{h}_1,0\le i\le m_i$$; 将区间 $[c,d]$ 作 n 等分, 记 $$h_2=(d-c)/n,y_j=c+j{h}_2,0\le j\le n$$称 $h_1$ 为 $x$ 方向的步长, $h_2$为 $y$ 方向上的步长, 用两簇平行线
$$x=x_i,0\le i\le m$$ $$y=y_j,0\le j\le n$$

将区域$\Omega$剖分为 $mn$个小矩形, 称两簇直线的交点 $(x_i,y_i)$ 为网格结点.
记
Ω h = { ( x i , y j ) ∣ 0 ≤ i ≤ m ,   0 ≤ j ≤ n } \Omega_h = \left\lbrace (x_i, y_j)| 0 \leq i \leq m, \ 0 \leq j \leq n \right\rbrace Ωh​={(xi​,yj​)∣0≤i≤m, 0≤j≤n}

称属于 $\Omega$的结点
Ω ˚ h = { ( x i , y j ) ∣ 1 ≤ i ≤ m − 1 , 1 ≤ j ≤ n − 1 } \mathring{\Omega}_h = \left\lbrace (x_i, y_j) | 1 \leq i \leq m-1, 1 \leq j \leq n-1 \right\rbrace Ω˚h​={(xi​,yj​)∣1≤i≤m−1,1≤j≤n−1}
称为内结点, 称位于 ${\Gamma }_h={\Omega }_h$ \ ${\mathring{\Omega }}_h$上的结点为边界结点

记
ω = { ( i , j ) ∣ ( x i , x j ) ∈ Ω ˚ h } , γ = { ( i , j ) ∣ ( x i , y j ) ∈ Γ h } \omega = \left\lbrace (i, j) | (x_i, x_j) \in \mathring{\Omega}_h \right\rbrace, \qquad \gamma = \left\lbrace (i,j) | (x_i, y_j) \in \Gamma_h \right\rbrace ω={(i,j)∣(xi​,xj​)∈Ω˚h​},γ={(i,j)∣(xi​,yj​)∈Γh​}

记
S h = { v ∣ v = { v i j ∣ 0 ≤ i ≤ m ,   0 ≤ j ≤ n } 为 Ω h 上的网格函数 } S_h = \left\lbrace v | v = \left\lbrace v_{ij} | 0 \leq i \leq m, \ 0 \leq j \leq n\right\rbrace \text{为} \Omega_h \text{上的网格函数}\right\rbrace Sh​={v∣v={vij​∣0≤i≤m, 0≤j≤n}为Ωh​上的网格函数}

设 v = { v i j   ∣ 0 ≤ i ≤ m , 0 ≤ j ≤ n } ∈ S h v = \left\lbrace v_{ij\ } | 0 \leq i \leq m, 0 \leq j \leq n \right\rbrace \in S_h v={vij ​∣0≤i≤m,0≤j≤n}∈Sh​ 给出如下记号:

$$D_x{v}_{ij}=\frac{1}{h_1}(v_{i+1,j}-v_{ij}),D_{\bar{x}}{v}_{ij}=\frac{1}{h_1}(v_{ij}-v_{i-1,j})$$$$D_y{v}_{ij}=\frac{1}{h_2}(v_{i,j+1}-v_{ij}),D_{\bar{y}}{v}_{ij}=\frac{1}{h_2}(v_{ij}-v_{i,j-1})$$ $${\delta }_x^2{v}_{ij}=\frac{1}{h_1}(D_x{v}_{ij}-D_{\bar{x}}{v}_{ij}),{\delta }_y^2{v}_{ij}=\frac{1}{h_2}(D_y{v}_{ij}-D_{\bar{y}}{v}_{ij})$$

$$\begin{gathered} \Vert v{\Vert }_{\infty }=\underset{0\le i\le m0 \\ \le j\le n}{\max }|v_{ij}| \end{gathered}$$ 称为 $v$ 的无穷范数.

下面在结点处考虑上述椭圆型方程边值问题:

$$-\left[\frac{{\partial }^{2}u}{\partial x^2}(x_i,y_j)+\frac{{\partial }^{2}u}{\partial y^2}(x_i,y_j)\right]=f(x_i,y_j),(i,j)\in \omega$$ $$u(x_i,y_j)=\phi (x_i,y_j),(i,j)\in \gamma$$

定义 ${\Omega }_h$上的网格函数:
U = { U i j ∣ 0 ≤ i ≤ m ,   0 ≤ j ≤ n } , U = \left\lbrace U_{ij}| 0 \leq i \leq m , \ 0 \leq j \leq n \right\rbrace, U={Uij​∣0≤i≤m, 0≤j≤n},
其中
$$U_{ij}=u(x_i,y_j),0\le i\le m,0\le j\le n.$$

将原方程中的二阶导数项利用差分估计得:

$$\frac{{\partial }^{2}u}{\partial x^2}(x_i,y_j)=\frac{1}{h_1^2}\left[u(x_{i-1},y_j)-2u(x_i,y_j)+u(x_{i+1},y_j)\right]-\frac{h_1^2}{12}\frac{{\partial }^{4}u({\xi }_{ij},y_j)}{\partial x^4}$$$$={\delta }_x^2{U}_{ij}-\frac{h_1^2}{12}\frac{{\partial }^{4}u({\xi }_{ij},y_j)}{\partial x^4},x_{i-1}<{\xi }_{ij}<x_{i+1}$$

$$\frac{{\partial }^{2}u}{\partial y^2}(x_i,y_j)=\frac{1}{h_2^2}\left[u(x_i,y_{j-1})-2u(x_i,y_j)+u(x_i,y_{j+1})\right]-\frac{h_2^2}{12}\frac{{\partial }^{4}u(x_i,{\eta }_{ij})}{\partial y^4}$$
$$={\delta }_y^2{U}_{ij}-\frac{h_2^2}{12}\frac{{\partial }^{4}u(x_i,{\eta }_{ij})}{\partial y^4},y_{j-1}<{\eta }_{ij}<y_{j+1}$$
代入原方程中, 并略去高阶小量项 $R_{ij}=-\frac{h_1^2}{12}\frac{{\partial }^{4}u({\xi }_{ij},y_j)}{\partial x^4}-\frac{h_2^2}{12}\frac{{\partial }^{4}u(x_i,{\eta }_{ij})}{\partial y^4}$ 之后利用 $u_{ij}$ 代替 $U_{ij}$得到如下差分格式:

$$-({\delta }_x^2{u}_{ij}+{\delta }_y^2{u}_{ij})=f(x_i,y_j),(i,j)\in \omega$$$$u_{ij}=\phi (x_i,y_j),(i,j)\in \gamma$$

差分格式的求解

改写原方程为矩阵方程

差分格式是以 { u i j ∣ 1 ≤ i ≤ m − 1 , 1 ≤ j ≤ n − 1 } \left\lbrace u_{ij} | 1 \leq i \leq m-1, 1 \leq j \leq n-1 \right\rbrace {uij​∣1≤i≤m−1,1≤j≤n−1} 为未知量的线性方程组

改写$-({\delta }_x^2{u}_{ij}+{\delta }_y^2{u}_{ij})=f(x_i,y_j)$为:

$$-\frac{1}{h_2^2}{u}_{i,j-1}-\frac{1}{h_1^2}{u}_{i-1,j}+2(\frac{1}{h_1^2}+\frac{1}{h_2^2})u_{ij}-\frac{1}{h_1^2}{u}_{i+1,j}-\frac{1}{h_2^2}{u}_{i,j+1}=f(x_i,y_j)$$$$1\le i\le m-1,1\le j\le n-1$$

记:
u j = ( u 1 j u 2 j ⋮ u m − 1 , j ) , 0 ≤ j ≤ n \pmb{u}_j =

$$\begin{pmatrix} u_{1j} \\ u_{2j} \\ \vdots \\ u_{m-1, j} \end{pmatrix}$$ , \quad 0 \leq j \leq n uuuj​=⎝⎜⎜⎜⎛​u1j​u2j​⋮um−1,j​​⎠⎟⎟⎟⎞​,0≤j≤n
改写方程为:
$${Du\text{Du}Du}_{j-1}+{Cu\text{Cu}Cu}_j+{Du\text{Du}Du}_{j+1}={f\text{f}f}_j,1\le j\le n-1$$
$$C\text{C}C=\left(\begin{array}{cccccc}2\left(\frac{1}{h_1^2}+\frac{1}{h_2^2}\right)& -\frac{1}{h_1^2}& & & & \\ -\frac{1}{h_1^2}& 2\left(\frac{1}{h_1^2}+\frac{1}{h_2^2}\right)& -\frac{1}{h_1^2}& & & \\ & \ddots & \ddots & \ddots & & \\ & & -\frac{1}{h_1^2}& 2\left(\frac{1}{h_1^2}+\frac{1}{h_2^2}\right)& -\frac{1}{h_1^2}& \\ & & & -\frac{1}{h_1^2}& 2\left(\frac{1}{h_1^2}+\frac{1}{h_2^2}\right)& \end{array}\right)$$

D = ( − 1 h 2 2 − 1 h 2 2 ⋱ − 1 h 2 2 − 1 h 2 2 ) f j = ( f ( x 1 , y j ) + 1 h 1 2 φ ( x 0 , y j ) f ( x 2 , y j ) ⋮ f ( x m − 2 , y j ) f ( x m − 1 , y j ) + 1 h 1 2 φ ( x m , y j ) ) \pmb{D} =

$$\begin{pmatrix} -\dfrac{1}{h_2^2} & & & & \\ & -\dfrac{1}{h_2^2} & & & \\ & & \ddots & & \\ & & & -\dfrac{1}{h_2^2} & \\ & & & & -\dfrac{1}{h_2^2}\end{pmatrix}$$ \qquad \pmb{ f }_j = $$\begin{pmatrix} f(x_1, y_j) + \dfrac{1}{h_1^2} \varphi(x_0, y_j) \\ f(x_2, y_j) \\ \vdots \\ f(x_{m-2}, y_j) \\ f(x_{m-1}, y_j) + \dfrac{1}{h_1^2} \varphi(x_m, y_j)\end{pmatrix}$$ DDD=⎝⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎛​−h22​1​​−h22​1​​⋱​−h22​1​​−h22​1​​⎠⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎞​f​f​​fj​=⎝⎜⎜⎜⎜⎜⎜⎜⎜⎛​f(x1​,yj​)+h12​1​φ(x0​,yj​)f(x2​,yj​)⋮f(xm−2​,yj​)f(xm−1​,yj​)+h12​1​φ(xm​,yj​)​⎠⎟⎟⎟⎟⎟⎟⎟⎟⎞​

进一步改写, 得到:
( C D D C D ⋱ ⋱ ⋱ D C D D C ) ( u 1 u 2 ⋮ u n − 2 u n − 1 ) = ( f 1 − D u 0 f 2 ⋮ f n − 2 f n − 1 − D u n )

$$\begin{pmatrix} \pmb{C} & \pmb{D} & & & \\ \pmb{D} & \pmb{C} & \pmb{D} & & \\ & \ddots & \ddots & \ddots & \\ & & \pmb{D} & \pmb{C} & \pmb{D} \\ & & & \pmb{D} & \pmb{C} \end{pmatrix}$$ $$\begin{pmatrix} \pmb{u}_1 \\ \pmb{u}_2 \\ \vdots \\ \pmb{u}_{n-2} \\ \pmb{u}_{n-1}\end{pmatrix}$$ = $$\begin{pmatrix} \pmb{f}_1 - \pmb{Du}_0 \\ \pmb{f}_2 \\ \vdots \\ \pmb{f}_{n-2} \\ \pmb{f}_{n-1} - \pmb{Du}_n\end{pmatrix}$$ ⎝⎜⎜⎜⎜⎛​CCCDDD​DDDCCC⋱​DDD⋱DDD​⋱CCCDDD​DDDCCC​⎠⎟⎟⎟⎟⎞​⎝⎜⎜⎜⎜⎜⎛​uuu1​uuu2​⋮uuun−2​uuun−1​​⎠⎟⎟⎟⎟⎟⎞​=⎝⎜⎜⎜⎜⎜⎛​f​f​​f1​−DuDuDu0​f​f​​f2​⋮f​f​​fn−2​f​f​​fn−1​−DuDuDun​​⎠⎟⎟⎟⎟⎟⎞​

五点差分方法的简单实现

显然, 一个比较简单的方法就是直接求解上述线性方程组.在本文中直接使用求解线性方程组的方法进行暴力求解，对于矩阵规模比较大的时候，一般会使用迭代算法或是其他快速算法对上述线性方程组进行求解，在后面的文章中如果有时间会陆续进行介绍。

稀疏矩阵

上述线性方程组的矩阵为稀疏矩阵, 需要引入相关的稀疏矩阵操作,节省内存.

$Python$ 中 稀疏矩阵结构体 在 scipy.sparse 模块下.

$Python$ 中 用于求解大型稀疏矩阵方程组的函数为 scipy.sparse.linalg 模块下的 spsolve() 函数.

分别介绍这两部分

Python中的稀疏矩阵数据结构

由于 使用 spsolve() 函数时,传入的矩阵希望是 $CSR$ 或 $CSC$格式的系数矩阵, 因此仅介绍两种格式的稀疏矩阵.

以下面矩阵为例,分别介绍其 $CSR$格式存储以及 $CSC$格式存储.
A = ( 0 0 0 10 21 0 33 0 0 0 3 0 12 1 0 4 ) A =

$$\begin{pmatrix} 0 & 0 & 0 & 10 \\ 21 & 0 & 33 & 0 \\ 0 & 0 & 3 & 0 \\ 12 & 1 & 0 & 4 \end{pmatrix}$$ A=⎝⎜⎜⎛​021012​0001​03330​10004​⎠⎟⎟⎞​
$CSR(压缩行稀疏行):$
调用格式:

csr_matrix((data, indices, indptr), shape=(a, b), dtype = float)
1

csr_matrix 详解:

$ CSR $格式的稀疏矩阵对上述矩阵分别利用三个数组 data, indices, indptr 如下:

data = [10, 21, 33, 3, 12, 1, 4]
indices = [3, 0, 2, 2, 0, 1 , 3]
indptr = [0, 1, 3, 4, 7]
1
2
3

其中:
data (数据): 表示存储的非零有数值 .
indices (列索引): 按照顺序存储非零有效数值的列标号.
indptr (行偏移量): 表示indices中行元素的起始位置和终止位置, 如果indptr[i] == indptr[i+1]则表示改行元素全为0, 如果初始矩阵有 m 行, 则 len(indptr) == m+1.

示范

from scipy.sparse import csr_matrix, csc_matrix
data = [10, 21, 33, 3, 12, 1, 4]
indices = [3, 0, 2, 2, 0, 1, 3]
indptr = [0, 1, 3, 4, 7]
A = csr_matrix((data, indices, indptr), shape = (4,4))
print(A)
1
2
3
4
5
6

稀疏矩阵中提供 toarray()方法将一个稀疏矩阵转换为一个普通矩阵.

示范

print(A.toarray())
1

即可将A转换为对应的numpy矩阵.
或者可以直接通过将满矩阵 B 传入 csr_matrix 函数来获取稀疏矩阵

通过调用稀疏矩阵对象的 data, indices,indptr 属性, 查看相应的数据

示范

print('data = ', C.data)
print('indices = ', C.indices)
print('indptr = ', C.indptr)
1
2
3

输出结果:

data =  [10 21 33  3 12  1  4]
indices =  [3 0 2 2 0 1 3]
indptr =  [0 1 3 4 7]
1
2
3

$CSC(压缩列稀疏行):$

调用格式:

csc_matrix((data, indices, indptr), shape=(a, b), dtype=float)
1

$CSC$稀疏矩阵的存储方式与$CSR$格式稀疏矩阵的存储方式大致相同,将$CSR$格式中的行(列) 换为列(行)即可.

Python 中用于快速创建稀疏矩阵的函数

scipy.sparse 模块下提供了以下函数用于快速创建稀疏矩阵, 调用格式如下:

sparse.eye(20, 20, format = 'csr')   # 创建稀疏格式的单位矩阵
sparse.spdiags(array, n, row, col,format = 'csr') # 创建对角矩阵,其中n表示在主对角线上方n个单位的对角线上创建
1
2

其余更多的函数见scipy的官方文档.

Python 中用于求解稀疏矩阵方程组的函数

对于稀疏矩阵方程组的求解,需要调用 scipy.sparse.linalg模块下的spsolve()函数, 函数使用格式如下:

spsolve(spA, b)
1

其中 spA 是 $CSR$ 或者 $CSC$ 格式存储的稀疏矩阵 .

$spsolve()使用示例:$

求解线性方程组:
( 0 0 0 10 21 0 33 0 0 0 3 0 12 1 0 4 ) X = ( 1 2 3 4 )

$$\begin{pmatrix} 0 & 0 & 0 & 10 \\ 21 & 0 & 33 & 0 \\ 0 & 0 & 3 & 0 \\ 12 & 1 & 0 & 4 \end{pmatrix}$$ X = $$\begin{pmatrix} 1 \\ 2 \\ 3 \\ 4\end{pmatrix}$$ ⎝⎜⎜⎛​021012​0001​03330​10004​⎠⎟⎟⎞​X=⎝⎜⎜⎛​1234​⎠⎟⎟⎞​
其中,左端项矩阵要求用csr_matrix或者csc_matrix进行创建.

from scipy.sparse.linalg import spsolve
data = [10, 21, 33, 3, 12, 1, 4]
indices = [3, 0, 2, 2, 0, 1, 3] 
indptr = [0, 1, 3, 4, 7]
spA = csr_matrix((data, indices, indptr), shape = (4, 4))
print(spA.toarray())
b = np.arange(1, 5)
spsolve(spA, b)
1
2
3
4
5
6
7
8

结果如下:

array([-1.47619048, 21.31428571,  1.        ,  0.1       ])
1

矩阵的 Kronecker 乘积

为了快速组装左端项矩阵,可以利用$kronecker$积.
假设给定一个大小为 $m_1\times m_2$ 的矩阵 $A$ 和一个大小为 $n_1\times n_2$的矩阵, 定义 $A$ 和 $B$ 之间的 $Kronecker$ 乘积 $A\otimes B$ 如下:

A ⊗ B = ( a 11 B a 12 B ⋯ a 1 m 2 B a 21 B a 22 B ⋯ a 2 m 2 B ⋮ ⋮ ⋱ ⋮ a m 1 1 B a m 1 2 B ⋯ a m 1 m 2 B ) A \otimes B =

$$\begin{pmatrix} a_{ 11 } B & a_{ 12 }B & \cdots & a_{1m_2}B \\ a_{ 21 } B & a_{ 22 }B & \cdots & a_{2m_2}B \\ \vdots & \vdots & \ddots & \vdots \\ a_{ m_1 1 } B & a_{ m_1 2 }B & \cdots & a_{m_1 m_2}B \\ \end{pmatrix}$$ A⊗B=⎝⎜⎜⎜⎛​a11​Ba21​B⋮am1​1​B​a12​Ba22​B⋮am1​2​B​⋯⋯⋱⋯​a1m2​​Ba2m2​​B⋮am1​m2​​B​⎠⎟⎟⎟⎞​
例如 : 给定向量 $\vec{u}=\left(\begin{array}{c}1\\ 2\end{array}\right)$, $\vec{v}=\left(\begin{array}{c}3\\ 4\end{array}\right)$ 容易计算:
$$\vec{u}\otimes \vec{v}=\left(\begin{array}{c}3\\ 4\\ 6\\ 8\end{array}\right)$$
有关$Kronecker$乘积的更多描述可以查看相关数学资料.

$Numpy$ 下的 $Kronecker$积

$numpy$ 模块下的 kron函数专门用于计算 两个矩阵的 $kronecker$乘积

示范:

import numpy as np
A = np.array([ [1], [2] ], dtype = np.float64)
B = np.array([ [3], [4] ], dtype = np.float64)
np.kron(A, B)
1
2
3
4

输出结果:

[[3.]
 [4.]
 [6.]
 [8.]]
1
2
3
4

A2 = np.array([ [1, 2], [3, 1] ])
B2 = np.array([ [0, 3], [2, 1] ])
np.kron(A2, B2)
1
2
3

输出结果

array([[0, 3, 0, 6],
       [2, 1, 4, 2],
       [0, 9, 0, 3],
       [6, 3, 2, 1]])
1
2
3
4

稀疏格式下的 kron 函数

scipy.sparse模块下提供了稀疏格式的 kron 函数, 返回两个稀疏格式矩阵的 kronecker 乘积.
示范:

from scipy import sparse
A2 = sparse.csr_matrix(A2)
B2 = sparse.csr_matrix(B2)
sparse.kron(A2, B2).toarray()
1
2
3
4

输出结果:

array([[0, 3, 0, 6],
      [2, 1, 4, 2],
      [0, 9, 0, 3],
      [6, 3, 2, 1]], dtype=int32)
1
2
3
4

五点差分格式数值实验

编写程序实现五点差分格式,并利用下面的问题进行数值模拟:
$$-\left(\frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^{2}u}{\partial y^2}\right)=({\pi }^2-1)e^{x}sin(\pi y),0<x<2,0<y<1$$$$u(0,y)=sin(\pi y),u(2,y)=e^{2}sin(\pi y),0\le y\le$$$$u(x,0)=0,u(x,1)=0,0\le x\le 2$$

要求如下:

• 分别在如下网格上对问题进行求解:
  $(h_1,h_2)/(x,y)=(\frac{1}{8},\frac{1}{8});(\frac{1}{16},\frac{1}{16});(\frac{1}{32},\frac{1}{32});(\frac{1}{64},\frac{1}{64})$
• 绘制$(h_1,h_2)/(x,y)=(\frac{1}{8},\frac{1}{8})$网格上的数值解曲面以及真解曲面
• 绘制各个网格上数值解的误差曲面
• 计算上述各个数值解与真解之间的全局最大模误差, 并验证上述数值格式所得到的数值解最大模误差阶为 $O(h_1^2+h_2^2)$

已知上述方程的真解为:$u(x,y)=e^{x}sin(\pi y)$

导入需要模块

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from scipy import sparse as ss
from scipy.sparse.linalg import spsolve
1
2
3
4
5

模型参数

class PdeModel:
    def __init__(self, x, y):
        self.x = x
        self.y = y

    def space_grid(self, Nx, Ny):
        X = np.linspace(self.x[0], self.x[1], Nx + 1)
        Y = np.linspace(self.y[0], self.y[1], Ny + 1)
        h1 = (self.x[1] - self.x[0]) / Nx
        h2 = (self.y[1] - self.y[0]) / Ny
        return h1, h2, X, Y

    def source(self, X, Y):
        return (np.pi ** 2 - 1) * np.exp(X) * np.sin(np.pi * Y)

    def solution(self, X, Y):
        YY, XX = np.meshgrid(Y, X)
        return np.exp(XX) * np.sin(np.pi * YY)

    def left_solution(self, Y):
        return np.sin(np.pi * Y)

    def right_solution(self, Y):
        return np.exp(2) * np.sin(np.pi * Y)

    def upper_solution(self, X):
        return np.zeros(len(X))

    def bottom_solution(self, X):
        return np.zeros(len(X))
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30

误差函数

def fd2d_epbvp_error(solution, uh, XY):
    ee = np.abs((solution(XY) - uh))
    emax = np.max(ee)
    return emax, ee
1
2
3
4

差分格式简单实现

def fd2d_epbvp(pde, NS):
   Nx = NS[0]
   Ny = NS[1]
   h1, h2, X, Y = pde.space_grid(Nx, Ny)
   YY, XX = np.meshgrid(Y, X)
   uh = np.zeros((Nx + 1, Ny + 1))
   uh[0, :] = pde.left_solution(Y)
   uh[-1, :] = pde.right_solution(Y)
   uh[:, 0] = pde.bottom_solution(X)
   uh[:, -1] = pde.upper_solution(X)
   tmp = np.ones(Nx - 1)
   D = ss.diags([2 * tmp, -tmp, -tmp], [0, -1, 1], format = 'csc', dtype = np.float64)
   C = np.diag(-tmp / h2 ** 2)
   Id = ss.eye(Ny - 1, format = 'csc')
   DD = ss.kron(Id, D, format = 'csc') / h1 ** 2 + ss.kron(D, Id, format = 'csc') / h2 ** 2
   rhs = pde.source(XX[1: -1, 1: -1], YY[1: -1, 1: -1])
   rhs[0, :] = rhs[0, :] + uh[0, 1:-1] / h1 ** 2
   rhs[-1, :] = rhs[-1, :] + uh[-1, 1:-1] / h1 ** 2
   rhs = rhs.T.flatten()
   rhs[0: Nx - 1] = rhs[0: Nx - 1] - np.dot(C, uh[1:-1, 0])
   rhs[1 - Nx:] = rhs[1 - Nx:] - np.dot(C, uh[1:-1, -1])
   uh[1:-1, 1:-1] = spsolve(DD, rhs).reshape(Ny - 1, Nx - 1).T
   return X, Y, uh
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23

测试程序

def fd2d_epbvp_test():
    # 初始化
    x = np.array([0, 2])
    y = np.array([0, 1])
    NS = [np.array([8, 8]), np.array([16, 16]), np.array([32, 32]), np.array([64, 64])]
    X_rec = []
    Y_rec = []
    XX_rec = []
    YY_rec = []
    uh_rec = []
    ee_rec = []
    emax_rec = np.zeros(4)

    for i in range(4):
        pde = PdeModel(x, y)
        X, Y, uh = fd2d_epbvp(pde, NS[i])
        YY, XX = np.meshgrid(Y, X)
        X_rec.append(X)
        Y_rec.append(Y)
        XX_rec.append(XX)
        YY_rec.append(YY)
        uh_rec.append(uh)
        emax, ee = fd2d_epbvp_error(pde.solution, uh, X, Y)
        ee_rec.append(ee)
        emax_rec[i] = emax
    # 在最后一个网格上绘制真解图像, 以及第一个网格上的数值解图像
    plt.figure(figsize=(15, 5))
    ax1 = plt.subplot(131, projection='3d')
    ax2 = plt.subplot(132, projection='3d')
    ax3 = plt.subplot(133, projection='3d')
    ax1.set_title('numeric solution')
    ax2.set_title('exact solution')
    ax3.set_title('error')
    ax3.view_init(elev=0, azim=40)
    
    ax1.plot_surface(XX_rec[0], YY_rec[0], uh_rec[0], cmap='gist_ncar')
    ax2.plot_surface(XX_rec[-1], YY_rec[-1], pde.solution(X_rec[-1], Y_rec[-1]), cmap='gist_ncar')
    
    for i in range(4):
        ax3.plot_surface(XX_rec[i], YY_rec[i], ee_rec[i], cmap = 'gist_ncar', vmin = 0, vmax = 0.05)
    plt.show()

    print('E(2h, 2h) / E(h, h)')
    for i in range(3):
        print(emax_rec[i] / emax_rec[i + 1])
        
if __name__ == '__main__':
    fd2d_epbvp_test()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48

求解结果:

误差阶

参考资料

[1]孙志忠.偏微分方程数值解法[M].北京:科学出版社,2015:34-42.