一维空间抛物型方程差分法求解-稳定性分析

稳定性和后验误差

问题提出
u t + a u x = b u x x + c u + f , x ∈ Ω , t > 0 , u_t + au_x = bu_{xx} + cu + f,x \in \Omega,t > 0, ut​+aux​=buxx​+cu+f,x∈Ω,t>0,
u ( x , 0 ) = u 0 ( x ) , x ∈ Ω , u(x,0) = u_0(x),x \in \Omega, u(x,0)=u0​(x),x∈Ω,
u ( α , t ) = u α ( t ) , t > 0 , u(\alpha,t) = u_{\alpha}(t),t > 0, u(α,t)=uα​(t),t>0,
u ( β , t ) = u β ( t ) , t > 0. u(\beta,t) = u_{\beta}(t),t > 0. u(β,t)=uβ​(t),t>0.

其中 a , b , c a,b,c a,b,c都是常数， u α ( t ) , u β ( t ) , u 0 ( x ) u_{\alpha}(t),u_{\beta}(t),u_0(x) uα​(t),uβ​(t),u0​(x)在 Ω = [ α , β ] \Omega = [\alpha,\beta] Ω=[α,β]的边界满足相容性条件。在我们接下来的报告中，我们只考虑 Ω = [ 0 , π ] , t ∈ [ 0 , 2 ] \Omega=[0,\pi],t\in[0,2] Ω=[0,π],t∈[0,2]这部分时空区域得数值解。针对上面的抛物方程，我们取 a = − 1 , b = 1 , c = 1 , α = 0 , β = π a = -1,b = 1,c = 1,\alpha = 0,\beta = \pi a=−1,b=1,c=1,α=0,β=π。那么我们求解的方程就是：
u t − u x = u x x + u + f , x ∈ Ω , t > 0 , u_t -u_x = u_{xx} + u + f,x \in \Omega,t > 0, ut​−ux​=uxx​+u+f,x∈Ω,t>0,
测试函数
我们先给出3个测试函数来进行稳定性分析和差分法的修正比较：
u 1 ( x , t ) = l n ( 10 ( x + t ) 2 + ( x − t ) 2 + 0.5 ) u_{1}(x,t) = ln(10(x + t)^{2} + (x - t)^{2} + 0.5) u1​(x,t)=ln(10(x+t)2+(x−t)2+0.5)
u 2 ( x , t ) = ( x 3 − x ) ∗ c o s h ( 2 t ) , u 3 ( x , t ) = x 2 − t 2 x 2 + t 2 + 0.1 u_{2}(x,t) = (x^{3} - x)*cosh(2t),\\ u_{3}(x,t) = \frac{x^{2} - t^{2}}{x^{2} + t^{2} + 0.1} u2​(x,t)=(x3−x)∗cosh(2t),u3​(x,t)=x2+t2+0.1x2−t2​
针对空间区域 Ω = [ 0 , π ] \Omega = [0,\pi] Ω=[0,π]以及时间区域 t i m e = [ 0 , 2 ] time= [0,2] time=[0,2]，我们在空间方向和时间方向分别取适当的分化长度 h , τ h,\tau h,τ，使得 J = π h + 1 , M = 2 τ + 1 J = \frac{\pi}{h} + 1,M = \frac{2}{\tau} + 1 J=hπ​+1,M=τ2​+1都为正整数。
引入一些记号：
e r r o r = ∑ ( j = 0 , m = 0 ) ( j = J − 1 , m = M − 1 ) ( u i , j − u e x a c t i , j ) 2 s u m ( j = 0 , m = 0 ) ( j = J − 1 , m = M − 1 ) ( u e x a c t i , j ) 2 error = \frac{\sum_{(j=0,m=0)}^{(j = J-1,m = M-1)} (u_{i,j} - uexact_{i,j})^{2} }{\\sum_{(j=0,m=0)}^{(j = J-1,m = M-1)} (uexact_{i,j})^{2}} error=sum(j=0,m=0)(j=J−1,m=M−1)​(uexacti,j​)2∑(j=0,m=0)(j=J−1,m=M−1)​(ui,j​−uexacti,j​)2​
L e r r o r ∞ = m a x ∣ u i , j − u e x a c t i , j ∣ L_{error}^{\infty} = max |u_{i,j} - uexact_{i,j}| Lerror∞​=max∣ui,j​−uexacti,j​∣
这两个是我们误差评估的公式。

显格式求解

U j m + 1 − U j m τ + a U j + 1 m − U j − 1 m 2 h = b U j + 1 m − 2 U j m + U j − 1 m h 2 + c U j m + f j m + 1 \frac{U_{j}^{m+1}-U_{j}^{m}}{\tau}+a \frac{U_{j+1}^{m}-U_{j-1}^{m}}{2 h}=b \frac{U_{j+1}^{m}-2 U_{j}^{m}+U_{j-1}^{m}}{h^{2}}+c U_{j}^{m} + f_{j}^{m+1} τUjm+1​−Ujm​​+a2hUj+1m​−Uj−1m​​=bh2Uj+1m​−2Ujm​+Uj−1m​​+cUjm​+fjm+1​
⇒ \Rightarrow ⇒
L E ( U j m + 1 ) = U j m + 1 − U j m τ + a U j + 1 m − U j − 1 m 2 h − b ( U j + 1 m − 2 U j m + U j − 1 m h 2 ) − c U j m = f j m + 1 L_{E}(U_{j}^{m+1}) = \frac{U_{j}^{m+1}-U_{j}^{m}}{\tau}+a \frac{U_{j+1}^{m}-U_{j-1}^{m}}{2 h}-b (\frac{U_{j+1}^{m}-2 U_{j}^{m}+U_{j-1}^{m}}{h^{2}})-c U_{j}^{m} = f_{j}^{m+1} LE​(Ujm+1​)=τUjm+1​−Ujm​​+a2hUj+1m​−Uj−1m​​−b(h2Uj+1m​−2Ujm​+Uj−1m​​)−cUjm​=fjm+1​
⇒ \Rightarrow ⇒
误差 e j m = U j m − u j m e_{j}^{m} = U_{j}^{m} - u_{j}^{m} ejm​=Ujm​−ujm​，其中 u j m u_{j}^{m} ujm​表示抛物方程的精确解 u u u在 ( x j , t m ) (x_j,t_m) (xj​,tm​)的取值。差分算子作用在 e j m e_{j}^{m} ejm​上，得到：
L E ( e j m + 1 ) = e j m + 1 − e j m τ + a e j + 1 m − e j − 1 m 2 h − b ( e j + 1 m − 2 e j m + e j − 1 m h 2 ) − c e j m = T j m L_{E}(e_{j}^{m+1}) = \frac{e_{j}^{m+1}-e_{j}^{m}}{\tau}+a \frac{e_{j+1}^{m}-e_{j-1}^{m}}{2 h}-b (\frac{e_{j+1}^{m}-2 e_{j}^{m}+e_{j-1}^{m}}{h^{2}})-c e_{j}^{m} = T_{j}^{m} LE​(ejm+1​)=τejm+1​−ejm​​+a2hej+1m​−ej−1m​​−b(h2ej+1m​−2ejm​+ej−1m​​)−cejm​=Tjm​
截断误差：
T j m = − u t t ( x j , t m ′ ) ∗ τ 2 − a ∂ 2 u x 2 ( x j ′ , t m ) ∗ h 2 6 + b ∂ 4 u x 4 ( x j ′ ′ , t m ) ∗ h 2 12 = O ( τ + h 2 ) T_{j}^{m} = -u_{tt}(x_j,t_{m}^{'})*\frac{\tau}{2} - a \frac{\partial^2 u}{x^2}(x_{j}^{'},t_m)*\frac{h^2}{6} + b \frac{\partial^4 u}{x^4}(x_{j}^{''},t_m)*\frac{h^2}{12}=O(\tau+h^2) Tjm​=−utt​(xj​,tm′​)∗2τ​−ax2∂2u​(xj′​,tm​)∗6h2​+bx4∂4u​(xj′′​,tm​)∗12h2​=O(τ+h2)
稳定性分析
U j m + 1 = ( μ + 1 2 ν ) U j − 1 m + ( 1 + c τ − 2 μ ) U j m + ( μ − 1 2 ν ) U j + 1 m U_{j}^{m + 1} = (\mu + \frac{1}{2}\nu)U_{j - 1}^{m} + (1 + c \tau - 2\mu)U_{j}^{m} + (\mu - \frac{1}{2}\nu)U_{j+1}^{m} Ujm+1​=(μ+21​ν)Uj−1m​+(1+cτ−2μ)Ujm​+(μ−21​ν)Uj+1m​
其中 μ = b τ h 2 , ν = a τ h \mu = \frac{b \tau}{h^2},\nu = \frac{a \tau}{h} μ=h2bτ​,ν=haτ​，很明显，根据差分格式，我们很容易得到，当 ∣ ν ∣ ≤ 2 μ ≤ 1 |\nu| \leq 2 \mu \leq 1 ∣ν∣≤2μ≤1时，显式差分格式满足最大值原理，因此是 L ∞ L^{\infty} L∞稳定的。在最大值原理满足的条件下，我们根据截断误差的分析公式进一步推导，我们有下面的不等式估计：
e j m + 1 = ( μ + ν 2 ) e j − 1 m + ( 1 + c τ − 2 μ ) e j m + ( μ − ν 2 ) e j + 1 m + τ T j m e_{j}^{m+1} = (\mu + \frac{\nu}{2})e_{j - 1}^{m} + (1 + c \tau - 2 \mu)e_{j}^{m} + (\mu - \frac{\nu}{2})e_{j + 1}^{m} + \tau T_{j}^{m} ejm+1​=(μ+2ν​)ej−1m​+(1+cτ−2μ)ejm​+(μ−2ν​)ej+1m​+τTjm​
⇒ \Rightarrow ⇒
∣ e j m + 1 ∣ ≤ ( μ + ν 2 ) ∣ e m ∣ + ( 1 + c τ − 2 μ ) ∣ e m ∣ + ( μ − ν 2 ) ∣ e m ∣ + τ ∣ T m ∣ |e_{j}^{m+1}| \leq (\mu + \frac{\nu}{2})|e^{m}| + (1 + c \tau - 2 \mu)|e^{m}| + (\mu - \frac{\nu}{2})|e^{m}| + \tau |T^{m}| ∣ejm+1​∣≤(μ+2ν​)∣em∣+(1+cτ−2μ)∣em∣+(μ−2ν​)∣em∣+τ∣Tm∣

⇒ \Rightarrow ⇒

∣ e m + 1 ∣ ≤ ( 1 + c τ ) ∣ e m ∣ + τ ∣ T m ∣ |e^{m+1}| \leq (1 + c \tau)|e^{m}| + \tau |T^{m}| ∣em+1∣≤(1+cτ)∣em∣+τ∣Tm∣

⇒ \Rightarrow ⇒

∣ e m ∣ ≤ ( 1 + c τ ) m ∣ e 0 ∣ + τ ∑ j = 0 m ∣ T j ∣ |e^{m}| \leq (1 + c \tau)^{m}|e^{0}| + \tau \sum_{j = 0}^{m}|T^{j}| ∣em∣≤(1+cτ)m∣e0∣+τ∑j=0m​∣Tj∣

⇒ \Rightarrow ⇒

∣ e m ∣ ≤ e x p ( c t m a x ) ∣ e 0 ∣ + t m a x ∣ T ∞ , t m a x ∣ |e^{m}| \leq exp(ct_{max})|e^{0}| + t_{max}|T_{\infty,t_{max}}| ∣em∣≤exp(ctmax​)∣e0∣+tmax​∣T∞,tmax​​∣
根据范数的等价性，显然这种格式当 ∣ ν ∣ ≤ 2 μ ≤ 1 |\nu| \leq 2 \mu \leq 1 ∣ν∣≤2μ≤1时也是 L 2 L^{2} L2稳定的，当边界上误差为0时，我们进一步可以得到显格式误差为 O ( τ + h 2 ) O(\tau + h^2) O(τ+h2)。
U j m + 1 = ( μ + 1 2 ν ) U j − 1 m + ( 1 + c τ − 2 μ ) U j m + ( μ − 1 2 ν ) U j + 1 m + τ f j m + 1 , A U m + 1 = B U m + r U_{j}^{m + 1} = (\mu + \frac{1}{2}\nu)U_{j - 1}^{m} + (1 + c \tau - 2\mu)U_{j}^{m} + (\mu - \frac{1}{2}\nu)U_{j+1}^{m} + \tau f_{j}^{m + 1}, AU^{m+1} = BU^{m} + r Ujm+1​=(μ+21​ν)Uj−1m​+(1+cτ−2μ)Ujm​+(μ−21​ν)Uj+1m​+τfjm+1​,AUm+1=BUm+r
， U m = ( U 1 m , U 2 m , … , U J − 2 m ) T U^{m} = (U_{1}^{m},U_{2}^{m},\ldots,U_{J - 2}^{m})^{T} Um=(U1m​,U2m​,…,UJ−2m​)T。A是 ( J − 2 ) × ( J − 2 ) (J - 2)\times(J -2) (J−2)×(J−2)的单位矩阵，B是三对角矩阵， r ∈ R J − 1 r \in R^{J-1} r∈RJ−1，\
B = ( 1 + c τ − 2 μ μ − 1 2 ν μ + 1 2 ν 1 + c τ − 2 μ μ − 1 2 ν ⋱ ⋱ ⋱ μ + 1 2 ν 1 + c τ − 2 μ μ − 1 2 ν μ + 1 2 ν 1 + c τ − 2 μ ) B=\left(\begin{array}{ccccc} 1 + c \tau - 2\mu & \mu - \frac{1}{2}\nu & & & \\ \mu + \frac{1}{2}\nu & 1 + c \tau - 2\mu & \mu - \frac{1}{2}\nu & & \\ & \ddots & \ddots & \ddots & \\ & & \mu + \frac{1}{2}\nu & 1 + c \tau - 2\mu & \mu - \frac{1}{2}\nu \\ & & & \mu + \frac{1}{2}\nu & 1 + c \tau - 2\mu \end{array}\right) B=⎝⎜⎜⎜⎜⎛​1+cτ−2μμ+21​ν​μ−21​ν1+cτ−2μ⋱​μ−21​ν⋱μ+21​ν​⋱1+cτ−2μμ+21​ν​μ−21​ν1+cτ−2μ​⎠⎟⎟⎟⎟⎞​ \
r = ( τ f 1 m + 1 + ( μ + 1 2 ν ) U 0 m , τ f 2 m + 1 , … , τ f J − 2 m + 1 + ( μ − 1 2 ν ) U J − 1 m ) T r = (\tau f_{1}^{m+1} + (\mu + \frac{1}{2}\nu)U_{0}^{m},\tau f_{2}^{m+1},\ldots,\tau f_{J-2}^{m+1} + (\mu - \frac{1}{2}\nu)U_{J-1}^{m})^{T} r=(τf1m+1​+(μ+21​ν)U0m​,τf2m+1​,…,τfJ−2m+1​+(μ−21​ν)UJ−1m​)T
不过在具体求解过程中，没必要使用矩阵。
针对我们设计的参数，方程\ref{pao}，我们给定在x方向的分化长度 δ x = h = 0.1 \delta x = h = 0.1 δx=h=0.1，那么根据a,b,c的选取，我们必须要取空间方向的分化 τ ≤ 0.005 ∗ π \tau \leq 0.005*\pi τ≤0.005∗π，这里我们取得就是 τ = 0.005 \tau = 0.005 τ=0.005，否则在迭代过程中会出现不稳定。事实上，在我们设计代码的过程中，我们可以自己做测试，就会发现当 τ > 0.005 ∗ π \tau > 0.005*\pi τ>0.005∗π时，误差打印出来就是inf。代码如下：

# -*- coding: utf-8 -*-
"""
Created on Sat Nov  7 19:17:15 2020


"""

import numpy as np
import matplotlib.pyplot as plt
import torch
import time
import torch.nn as nn
st = time.time()
def UU(X, order,prob):#X表示(x,t)
    if prob==1:
        temp = 10*(X[:,0]+X[:,1])**2 + (X[:,0]-X[:,1])**2 + 0.5
        if order[0]==0 and order[1]==0:
            return torch.log(temp)
        if order[0]==1 and order[1]==0:#对x求偏导
            return temp**(-1) * (20*(X[:,0]+X[:,1]) + 2*(X[:,0]-X[:,1]))
        if order[0]==0 and order[1]==1:#对t求偏导
            return temp**(-1) * (20*(X[:,0]+X[:,1]) - 2*(X[:,0]-X[:,1]))
        if order[0]==2 and order[1]==0:
            return - temp**(-2) * (20*(X[:,0]+X[:,1])+2*(X[:,0]-X[:,1])) ** 2 \
                   + temp**(-1) * (22)
        if order[0]==1 and order[1]==1:
            return - temp**(-2) * (20*(X[:,0]+X[:,1])+2*(X[:,0]-X[:,1])) \
                   * (20*(X[:,0]+X[:,1])-2*(X[:,0]-X[:,1])) \
                   + temp**(-1) * (18)
        if order[0]==0 and order[1]==2:
            return - temp**(-2) * (20*(X[:,0]+X[:,1])-2*(X[:,0]-X[:,1])) ** 2 \
                   + temp**(-1) * (22)
    if prob==2:
        if order[0]==0 and order[1]==0:
            return (X[:,0]*X[:,0]*X[:,0]-X[:,0]) * \
                   0.5*(torch.exp(2*X[:,1])+torch.exp(-2*X[:,1]))
        if order[0]==1 and order[1]==0:
            return (3*X[:,0]*X[:,0]-1) * \
                   0.5*(torch.exp(2*X[:,1])+torch.exp(-2*X[:,1]))
        if order[0]==0 and order[1]==1:
            return (X[:,0]*X[:,0]*X[:,0]-X[:,0]) * \
                   (torch.exp(2*X[:,1])-torch.exp(-2*X[:,1]))
        if order[0]==2 and order[1]==0:
            return (6*X[:,0]) * \
                   0.5*(torch.exp(2*X[:,1])+torch.exp(-2*X[:,1]))
        if order[0]==1 and order[1]==1:
            return (3*X[:,0]*X[:,0]-1) * \
                   (torch.exp(2*X[:,1])-torch.exp(-2*X[:,1]))
        if order[0]==0 and order[1]==2:
            return (X[:,0]*X[:,0]*X[:,0]-X[:,0]) * \
                   2*(torch.exp(2*X[:,1])+torch.exp(-2*X[:,1]))
    if prob==3:
        temp1 = X[:,0]*X[:,0] - X[:,1]*X[:,1]
        temp2 = X[:,0]*X[:,0] + X[:,1]*X[:,1] + 0.1
        if order[0]==0 and order[1]==0:
            return temp1 * temp2**(-1)
        if order[0]==1 and order[1]==0:
            return (2*X[:,0]) * temp2**(-1) + \
                   temp1 * (-1)*temp2**(-2) * (2*X[:,0])
        if order[0]==0 and order[1]==1:
            return (-2*X[:,1]) * temp2**(-1) + \
                   temp1 * (-1)*temp2**(-2) * (2*X[:,1])
        if order[0]==2 and order[1]==0:
            return (2) * temp2**(-1) + \
                   2 * (2*X[:,0]) * (-1)*temp2**(-2) * (2*X[:,0]) + \
                   temp1 * (2)*temp2**(-3) * (2*X[:,0])**2 + \
                   temp1 * (-1)*temp2**(-2) * (2)
        if order[0]==1 and order[1]==1:
            return (2*X[:,0]) * (-1)*temp2**(-2) * (2*X[:,1]) + \
                   (-2*X[:,1]) * (-1)*temp2**(-2) * (2*X[:,0]) + \
                   temp1 * (2)*temp2**(-3) * (2*X[:,0]) * (2*X[:,1])
        if order[0]==0 and order[1]==2:
            return (-2) * temp2**(-1) + \
                   2 * (-2*X[:,1]) * (-1)*temp2**(-2) * (2*X[:,1]) + \
                   temp1 * (2)*temp2**(-3) * (2*X[:,1])**2 + \
                   temp1 * (-1)*temp2**(-2) * (2)
def FF(prob,X,a,b,c):
    return UU(X,[0,1],prob) + a*UU(X,[1,0],prob) - b*UU(X,[2,0],prob) - c*UU(X,[0,0],prob)
class FD():
    def __init__(self,bound,hx,prob):
        self.prob = prob
        self.dim = 2
        self.hx = hx
        self.nx = [int((bound[0,1] - bound[0,0])/self.hx[0]) + 1,int((bound[1,1] - bound[1,0])/self.hx[1]) + 1]
        self.size = self.nx[0]*self.nx[1]
        self.X = torch.zeros(self.size,self.dim)
        for m in range(self.nx[1]):
            for j in range(self.nx[0]):
                self.X[m*self.nx[0] + j,0] = bound[0,0] + j*self.hx[0]
                self.X[m*self.nx[0] + j,1] = bound[1,0] + m*self.hx[1]
        self.u_acc = UU(self.X,[0,0],prob).view(-1,1)
        self.u_init()
    def u_init(self):
        self.u = torch.zeros(self.size,1)
        for m in range(self.nx[1]):
            for j in range(self.nx[0]):
                X = self.X[m*self.nx[0] + j:m*self.nx[0] + j + 1,:]
                if m == 0 or j == 0 or j == self.nx[0] - 1:
                    self.u[m*self.nx[0] + j] = UU(X,[0,0],self.prob)
    def solve(self,a,b,c):
        mu = b*self.hx[1]/(self.hx[0]**2)
        nu = 0.5*a*self.hx[1]/self.hx[0]
        for m in range(self.nx[1] - 1):
            for j in range(1,self.nx[0] - 1):
                X = self.X[m*self.nx[0] + j:m*self.nx[0] + j + 1,:]
                ind = m*self.nx[0] + j
                self.u[(m + 1)*self.nx[0] + j] = (mu + nu)*self.u[ind - 1] + (1 + c*self.hx[1] - 2*mu)*self.u[ind] + (mu - nu)*self.u[ind + 1]\
                + self.hx[1]*FF(self.prob,self.X[(m + 1)*self.nx[0] + j:(m + 1)*self.nx[0] + j + 1,:],a,b,c)
        return self.u
def error(u_pred, u_acc):
    fenzi = ((u_pred - u_acc)**2).sum()
    fenmu = (u_acc**2).sum()
    return (fenzi/fenmu)**(0.5)     
              
bound = torch.tensor([[0,np.pi],[0,2]]).float()
hx = [0.1*np.pi,0.005]
prob = 1
fd = FD(bound,hx,prob)
a = -1
b = 1
c = 1
u_acc = fd.u_acc
u_pred = fd.solve(a,b,c)
ela = time.time() - st
print("the prob is %d, error is %.4f,run time is %.3f"%(prob,error(u_pred,u_acc),ela))
print(max(abs(u_acc - u_pred)))
#--------------------------------------------------
#下面是做后验误差分析过程
def host_error(hx,h):
    fd = FD(bound,hx,prob)
    ffd = FD(bound,h,prob)
    J = int(hx[0]/h[0]);M = int(hx[1]/h[1])
    
    fd_u = fd.solve(a,b,c)
    ffd_u = ffd.solve(a,b,c)
    u = torch.zeros(fd.size,1)
    
    for m in range(fd.nx[1]):
        for j in range(fd.nx[0]):
            u[m*fd.nx[0] + j] = ffd_u[M*m*ffd.nx[0] + J*j]
    return max(abs(u - fd_u))
#------------------------------------
h = 0.1*np.pi
t = 0.0025
ex = [h,t]
print("ex:%.3f"%(t/(h**2)))
ed = FD(bound,ex,prob)

fx = [0.5*h,t]
fd = FD(bound,fx,prob)
print("fx:%.3f"%(t/((0.5*h)**2)))

gx = [0.25*h,t]
gd = FD(bound,gx,prob)
print("gx:%.3f"%(t/((0.25*h)**2)))

err = host_error(ex,fx)
er = host_error(fx,gx)
print(err,er)
print(err/er)

代码最后一部分是在计算后验误差，这个在没有解析精确解情况下特别重要，后文会提到。

隐格式求解

U j m + 1 − U j m τ + a U j + 1 m − U j − 1 m 2 h = b U j + 1 m + 1 − 2 U j m + 1 + U j − 1 m + 1 h 2 + c U j m + 1 + f j m + 1 \frac{U_{j}^{m+1}-U_{j}^{m}}{\tau}+a \frac{U_{j+1}^{m}-U_{j-1}^{m}}{2 h}=b \frac{U_{j+1}^{m+1}-2 U_{j}^{m+1}+U_{j-1}^{m+1}}{h^{2}}+c U_{j}^{m+1} + f_{j}^{m+1} τUjm+1​−Ujm​​+a2hUj+1m​−Uj−1m​​=bh2Uj+1m+1​−2Ujm+1​+Uj−1m+1​​+cUjm+1​+fjm+1​
L I ( U j m + 1 ) = U j m + 1 − U j m τ + a U j + 1 m − U j − 1 m 2 h − b ( U j + 1 m + 1 − 2 U j m + 1 + U j − 1 m + 1 h 2 ) − c U j m + 1 = f j m + 1 L_{I}(U_{j}^{m+1}) = \frac{U_{j}^{m+1}-U_{j}^{m}}{\tau}+a \frac{U_{j+1}^{m}-U_{j-1}^{m}}{2 h}-b (\frac{U_{j+1}^{m+1}-2 U_{j}^{m+1}+U_{j-1}^{m+1}}{h^{2}})-c U_{j}^{m+1} = f_{j}^{m+1} LI​(Ujm+1​)=τUjm+1​−Ujm​​+a2hUj+1m​−Uj−1m​​−b(h2Uj+1m+1​−2Ujm+1​+Uj−1m+1​​)−cUjm+1​=fjm+1​
隐格式的截断误差仍然是 O ( τ + h 2 ) O(\tau + h^2) O(τ+h2)
稳定性分析
( 1 − c τ + 2 μ ) U j m + 1 − μ ( U j − 1 m + 1 + U j + 1 m + 1 ) = ν 2 U j − 1 m + U j m − ν 2 U j + 1 m (1 - c \tau + 2 \mu)U_{j}^{m + 1} - \mu (U_{j - 1}^{m + 1} + U_{j + 1}^{m + 1}) = \frac{\nu}{2} U_{j - 1}^{m} + U_{j}^{m} - \frac{\nu}{2} U_{j + 1}^{m} (1−cτ+2μ)Ujm+1​−μ(Uj−1m+1​+Uj+1m+1​)=2ν​Uj−1m​+Ujm​−2ν​Uj+1m​
其中 μ = b τ h 2 , ν = a τ h \mu = \frac{b \tau}{h^2},\nu = \frac{a \tau}{h} μ=h2bτ​,ν=haτ​,
令 U j m = λ k m e x p ( i k π δ x ) U_{j}^{m} = \lambda_{k}^{m}exp(ik \pi \delta x) Ujm​=λkm​exp(ikπδx) 以及令 f = 0 f = 0 f=0，
( 1 − c τ + 2 μ ( 1 − c o s ( 2 β ) ) λ k m + 1 = ( 1 − i ν s i n ( 2 β ) ) λ k m (1 - c \tau + 2\mu (1 - cos(2\beta))\lambda_{k}^{m+1} = (1 - i\nu sin(2\beta))\lambda_{k}^{m} (1−cτ+2μ(1−cos(2β))λkm+1​=(1−iνsin(2β))λkm​
⇒ \Rightarrow ⇒

λ k = 1 − i ν s i n ( 2 β ) 1 − c τ + 2 μ ( 1 − c o s ( 2 β ) ) = 1 − i 2 ν s i n ( β ) c o s ( β ) 1 − c τ + 4 μ s i n 2 ( β ) \lambda_{k} = \frac{1 - i\nu sin(2\beta)}{1 - c \tau + 2\mu (1 - cos(2\beta))}= \frac{1 - i 2\nu sin(\beta)cos(\beta)}{1 - c \tau + 4\mu sin^2(\beta)} λk​=1−cτ+2μ(1−cos(2β))1−iνsin(2β)​=1−cτ+4μsin2(β)1−i2νsin(β)cos(β)​

其中 β = k π δ x 2 \beta = \frac{k\pi \delta x}{2} β=2kπδx​，
⇒ \Rightarrow ⇒

( λ k 2 − 1 ) 1 − c τ + 2 μ ( 1 − c o s ( 2 β ) ) 2 = ( 1 − i ν s i n ( 2 β ) ) 2 − ( 1 − c τ + 2 μ ( 1 − c o s ( 2 β ) ) ) 2 (\lambda_{k}^2 - 1)1 - c \tau + 2\mu (1 - cos(2\beta))^2 = (1 - i\nu sin(2\beta))^2 - (1 - c \tau + 2\mu (1 - cos(2\beta)))^2 (λk2​−1)1−cτ+2μ(1−cos(2β))2=(1−iνsin(2β))2−(1−cτ+2μ(1−cos(2β)))2

省略 c τ c\tau cτ，有
⇒ \Rightarrow ⇒
( λ k 2 − 1 ) 1 − c τ + 2 μ ( 1 − c o s ( 2 β ) ) 2 = 4 s i n 2 ( β ) ( 4 μ 2 s i n 2 ( β ) + 2 μ − ν 2 c o s 2 ( β ) ) (\lambda_{k}^2 - 1)1 - c \tau + 2\mu (1 - cos(2\beta))^2 = 4sin^2(\beta)(4\mu^2 sin^2(\beta) + 2\mu - \nu^2 cos^2(\beta)) (λk2​−1)1−cτ+2μ(1−cos(2β))2=4sin2(β)(4μ2sin2(β)+2μ−ν2cos2(β))

⇒ \Rightarrow ⇒

ν 2 ≥ 4 μ 2 s i n 2 ( β ) + 2 μ c o s 2 ( β ) \nu^2 \geq \frac{4\mu^2 sin^2(\beta) + 2 \mu}{cos^2(\beta)} ν2≥cos2(β)4μ2sin2(β)+2μ​
隐格式无条件稳定。
隐格式求解
( 1 − c τ + 2 μ ) U j m + 1 − μ ( U j − 1 m + 1 + U j + 1 m + 1 ) = ν 2 U j − 1 m + U j m − ν 2 U j + 1 m + τ f j m + 1 (1 - c \tau + 2 \mu)U_{j}^{m + 1} - \mu (U_{j - 1}^{m + 1} + U_{j + 1}^{m + 1}) = \frac{\nu}{2} U_{j - 1}^{m} + U_{j}^{m} - \frac{\nu}{2} U_{j + 1}^{m} + \tau f_{j}^{m + 1} (1−cτ+2μ)Ujm+1​−μ(Uj−1m+1​+Uj+1m+1​)=2ν​Uj−1m​+Ujm​−2ν​Uj+1m​+τfjm+1​
A = ( 1 − c τ + 2 μ − μ − μ 1 − c τ + 2 μ − μ ⋱ ⋱ ⋱ − μ 1 − c τ + 2 μ − μ − μ 1 − c τ + 2 μ ) A=\left(\begin{array}{ccccc} 1 - c \tau + 2 \mu & -\mu& & & \\ -\mu &1 - c \tau + 2 \mu & -\mu & & \\ & \ddots & \ddots & \ddots & \\ & & -\mu & 1 - c \tau + 2 \mu & -\mu\\ & & & -\mu & 1 - c \tau + 2\mu \end{array}\right) A=⎝⎜⎜⎜⎜⎛​1−cτ+2μ−μ​−μ1−cτ+2μ⋱​−μ⋱−μ​⋱1−cτ+2μ−μ​−μ1−cτ+2μ​⎠⎟⎟⎟⎟⎞​ \
B = ( 1 − ν 2 ν 2 1 − ν 2 ⋱ ⋱ ⋱ ν 2 1 − ν 2 ν 2 1 ) B=\left(\begin{array}{ccccc} 1 & -\frac{\nu}{2} & & & \\ \frac{\nu}{2}& 1 & -\frac{\nu}{2} & & \\ & \ddots & \ddots & \ddots & \\ & & \frac{\nu}{2} & 1 & -\frac{\nu}{2} \\ & & & \frac{\nu}{2} & 1 \end{array}\right) B=⎝⎜⎜⎜⎜⎛​12ν​​−2ν​1⋱​−2ν​⋱2ν​​⋱12ν​​−2ν​1​⎠⎟⎟⎟⎟⎞​ \
r = ( τ f 1 m + 1 + ν 2 U 0 m , τ f 2 m + 1 , … , τ f J − 2 m + 1 − ν 2 U J − 1 m ) T r = (\tau f_{1}^{m+1} + \frac{\nu}{2} U_{0}^{m},\tau f_{2}^{m + 1},\ldots,\tau f_{J-2}^{m + 1} - \frac{\nu}{2} U_{J-1}^{m} )^{T} r=(τf1m+1​+2ν​U0m​,τf2m+1​,…,τfJ−2m+1​−2ν​UJ−1m​)T
代码如下：

# -*- coding: utf-8 -*-
"""
Created on Sun Nov  8 09:38:27 2020


"""

import numpy as np
import matplotlib.pyplot as plt
import torch
import time
import torch.nn as nn
import torch.nn.functional as F
st = time.time()
def UU(X, order,prob):#X表示(x,t)
    if prob==1:
        temp = 10*(X[:,0]+X[:,1])**2 + (X[:,0]-X[:,1])**2 + 0.5
        if order[0]==0 and order[1]==0:
            return torch.log(temp)
        if order[0]==1 and order[1]==0:#对x求偏导
            return temp**(-1) * (20*(X[:,0]+X[:,1]) + 2*(X[:,0]-X[:,1]))
        if order[0]==0 and order[1]==1:#对t求偏导
            return temp**(-1) * (20*(X[:,0]+X[:,1]) - 2*(X[:,0]-X[:,1]))
        if order[0]==2 and order[1]==0:
            return - temp**(-2) * (20*(X[:,0]+X[:,1])+2*(X[:,0]-X[:,1])) ** 2 \
                   + temp**(-1) * (22)
        if order[0]==1 and order[1]==1:
            return - temp**(-2) * (20*(X[:,0]+X[:,1])+2*(X[:,0]-X[:,1])) \
                   * (20*(X[:,0]+X[:,1])-2*(X[:,0]-X[:,1])) \
                   + temp**(-1) * (18)
        if order[0]==0 and order[1]==2:
            return - temp**(-2) * (20*(X[:,0]+X[:,1])-2*(X[:,0]-X[:,1])) ** 2 \
                   + temp**(-1) * (22)
    if prob==2:
        if order[0]==0 and order[1]==0:
            return (X[:,0]*X[:,0]*X[:,0]-X[:,0]) * \
                   0.5*(torch.exp(2*X[:,1])+torch.exp(-2*X[:,1]))
        if order[0]==1 and order[1]==0:
            return (3*X[:,0]*X[:,0]-1) * \
                   0.5*(torch.exp(2*X[:,1])+torch.exp(-2*X[:,1]))
        if order[0]==0 and order[1]==1:
            return (X[:,0]*X[:,0]*X[:,0]-X[:,0]) * \
                   (torch.exp(2*X[:,1])-torch.exp(-2*X[:,1]))
        if order[0]==2 and order[1]==0:
            return (6*X[:,0]) * \
                   0.5*(torch.exp(2*X[:,1])+torch.exp(-2*X[:,1]))
        if order[0]==1 and order[1]==1:
            return (3*X[:,0]*X[:,0]-1) * \
                   (torch.exp(2*X[:,1])-torch.exp(-2*X[:,1]))
        if order[0]==0 and order[1]==2:
            return (X[:,0]*X[:,0]*X[:,0]-X[:,0]) * \
                   2*(torch.exp(2*X[:,1])+torch.exp(-2*X[:,1]))
    if prob==3:
        temp1 = X[:,0]*X[:,0] - X[:,1]*X[:,1]
        temp2 = X[:,0]*X[:,0] + X[:,1]*X[:,1] + 0.1
        if order[0]==0 and order[1]==0:
            return temp1 * temp2**(-1)
        if order[0]==1 and order[1]==0:
            return (2*X[:,0]) * temp2**(-1) + \
                   temp1 * (-1)*temp2**(-2) * (2*X[:,0])
        if order[0]==0 and order[1]==1:
            return (-2*X[:,1]) * temp2**(-1) + \
                   temp1 * (-1)*temp2**(-2) * (2*X[:,1])
        if order[0]==2 and order[1]==0:
            return (2) * temp2**(-1) + \
                   2 * (2*X[:,0]) * (-1)*temp2**(-2) * (2*X[:,0]) + \
                   temp1 * (2)*temp2**(-3) * (2*X[:,0])**2 + \
                   temp1 * (-1)*temp2**(-2) * (2)
        if order[0]==1 and order[1]==1:
            return (2*X[:,0]) * (-1)*temp2**(-2) * (2*X[:,1]) + \
                   (-2*X[:,1]) * (-1)*temp2**(-2) * (2*X[:,0]) + \
                   temp1 * (2)*temp2**(-3) * (2*X[:,0]) * (2*X[:,1])
        if order[0]==0 and order[1]==2:
            return (-2) * temp2**(-1) + \
                   2 * (-2*X[:,1]) * (-1)*temp2**(-2) * (2*X[:,1]) + \
                   temp1 * (2)*temp2**(-3) * (2*X[:,1])**2 + \
                   temp1 * (-1)*temp2**(-2) * (2)
def FF(prob,X,a,b,c):
    return UU(X,[0,1],prob) + a*UU(X,[1,0],prob) - b*UU(X,[2,0],prob) - c*UU(X,[0,0],prob)
class FD():
    def __init__(self,bound,hx,prob):
        self.prob = prob
        self.dim = 2
        self.hx = hx
        self.nx = [int((bound[0,1] - bound[0,0])/self.hx[0]) + 1,int((bound[1,1] - bound[1,0])/self.hx[1]) + 1]
        self.size = self.nx[0]*self.nx[1]
        self.X = torch.zeros(self.size,self.dim)
        for m in range(self.nx[1]):
            for j in range(self.nx[0]):
                self.X[m*self.nx[0] + j,0] = bound[0,0] + j*self.hx[0]
                self.X[m*self.nx[0] + j,1] = bound[1,0] + m*self.hx[1]
        self.u_acc = UU(self.X,[0,0],prob).view(-1,1)
        
    def matrix(self,a,b,c):
        mu = b*self.hx[1]/(self.hx[0]**2)
        nu = 0.5*a*self.hx[1]/self.hx[0]
        A = torch.zeros(self.nx[0],self.nx[0])
        B = torch.zeros(self.nx[0],self.nx[0])
        for i in range(1,self.nx[0] - 1):
            A[i,i] = 1 - c*self.hx[1] + 2*mu
            A[i,i + 1] = - mu
            A[i,i - 1] = - mu
            B[i,i + 1] = - nu
            B[i,i - 1] = nu
            B[i,i] = 1
        A[0,0] = 1;A[-1,-1] = 1
        return A,B
    def right(self,a,b,c,m):
        r = torch.zeros(self.nx[0],1)
        for j in range(1,self.nx[0] - 1):
            X = self.X[(m + 1)*self.nx[0] + j:(m + 1)*self.nx[0] + j + 1,:]
            r[j] = self.hx[1]*FF(self.prob,X,a,b,c)
        r[0] = UU(self.X[(m + 1)*self.nx[0]:(m + 1)*self.nx[0] + 1,:],[0,0],self.prob)
        r[-1] = UU(self.X[(m + 1)*self.nx[0] + self.nx[0] - 1:(m + 1)*self.nx[0] + self.nx[0],:],[0,0],self.prob)
        return r
    def solve(self,a,b,c):
        A,B = self.matrix(a,b,c)
        u = torch.zeros(self.size,1)
        
        u[:self.nx[0]] = UU(self.X[:self.nx[0],:],[0,0],self.prob).view(-1,1)
        #print(u[:self.nx[0]].shape, UU(self.X[:self.nx[0],:],[0,0],self.prob).shape)
        for m in range(self.nx[1] - 1):
            rig = [e-mail][m*self.nx[0]:(m + 1)*self.nx[0]] + self.right(a,b,c,m)
            u[(m + 1)*self.nx[0]:(m + 2)*self.nx[0]],lu = torch.solve(rig,A)
        return u
def error(u_pred, u_acc):
    fenzi = ((u_pred - u_acc)**2).sum()
    fenmu = (u_acc**2).sum()
    return (fenzi/fenmu)**(0.5)     
              
bound = torch.tensor([[0,np.pi],[0,2]]).float()
hx = [0.1*np.pi,0.05]
prob = 1
fd = FD(bound,hx,prob)
a = -1
b = 1
c = 1
u_acc = fd.u_acc
u_pred = fd.solve(a,b,c)
ela = time.time() - st
print("the prob is %d, error is %.4f,run time is %.3f"%(prob,error(u_pred,u_acc),ela))
print(max(abs(u_acc - u_pred)))
#--------------------------------------------------
#下面是做后验误差分析过程
def host_error(hx,h):
    fd = FD(bound,hx,prob)
    ffd = FD(bound,h,prob)
    J = int(hx[0]/h[0]);M = int(hx[1]/h[1])
    
    fd_u = fd.solve(a,b,c)
    ffd_u = ffd.solve(a,b,c)
    u = torch.zeros(fd.size,1)
    
    for m in range(fd.nx[1]):
        for j in range(fd.nx[0]):
            u[m*fd.nx[0] + j] = ffd_u[M*m*ffd.nx[0] + J*j]
    return max(abs(u - fd_u))
#------------------------------------
h = 0.1*np.pi
t = 0.05
ex = [h,t]
print("ex:%.3f"%(t/(h**2)))
ed = FD(bound,ex,prob)

fx = [0.5*h,t]
fd = FD(bound,fx,prob)
print("fx:%.3f"%(t/((0.5*h)**2)))

gx = [0.25*h,t]
gd = FD(bound,gx,prob)
print("gx:%.3f"%(t/((0.25*h)**2)))

err = host_error(ex,fx)
er = host_error(fx,gx)
print(err,er)
print(err/er)

无解析解情况

此前我们一直都是求解那种能给出精确解的方程，并且通过差分解和精确解比较来测试我们的差分法效果。如果都有精确解了，还要差分法干什么,这里我们重点引入后验误差估计的方法来评估差分法的优劣。
不管是哪种格式，都有下面的公式成立：
U h , τ , ( j ) m = u j m + ψ j m τ + φ j m h 2 + O ( τ 2 + h 4 ) U_{h,\tau,(j)}^{m} = u_{j}^{m} + \psi_{j}^{m} \tau + \varphi_{j}^{m} h^2 + O(\tau^2 + h^4) Uh,τ,(j)m​=ujm​+ψjm​τ+φjm​h2+O(τ2+h4)
则如果我们取分化长度为 ( τ 2 , h ) , ( τ 4 , h ) , ( τ , h 2 ) , ( τ , h 4 ) , ( τ 2 , h 2 ) (\frac{\tau}{2},h),(\frac{\tau}{4},h),(\tau,\frac{h}{2}),(\tau,\frac{h}{4}),(\frac{\tau}{2},\frac{h}{2}) (2τ​,h),(4τ​,h),(τ,2h​),(τ,4h​),(2τ​,2h​)，我们会得到形式一样的公式，具体参考下面：
⇒ \Rightarrow ⇒
U h 2 , τ , ( j ) m = u j m + ψ j m τ + φ j m ( h 2 ) 2 + O ( τ 2 + h 4 ) U_{\frac{h}{2},\tau,(j)}^{m} = u_{j}^{m} + \psi_{j}^{m} \tau + \varphi_{j}^{m} (\frac{h}{2})^2 + O(\tau^2 + h^4) U2h​,τ,(j)m​=ujm​+ψjm​τ+φjm​(2h​)2+O(τ2+h4)
⇒ \Rightarrow ⇒
U h 4 , τ , ( j ) m = u j m + ψ j m τ + φ j m ( h 4 ) 2 + O ( τ 2 + h 4 ) U_{\frac{h}{4},\tau,(j)}^{m} = u_{j}^{m} + \psi_{j}^{m} \tau + \varphi_{j}^{m} (\frac{h}{4})^2 + O(\tau^2 + h^4) U4h​,τ,(j)m​=ujm​+ψjm​τ+φjm​(4h​)2+O(τ2+h4)
我们做一个记号：用 U h 2 , m U_{\frac{h}{2},m} U2h​,m​表示分化长度为 h 2 , τ \frac{h}{2},\tau 2h​,τ的网格剖分上的数值解， u j m u_{j}^{m} ujm​表示当分化长度为 h , τ h,\tau h,τ在对应网格点上的精确解， u j u_j uj​表示分化长度为 h , τ h,\tau h,τ的整个网格上的精确解（相当于是精确的网格函数）。
U h , m − u j = 4 3 ( U h , m − U h 2 , m ) + O ( τ + h 4 ) U_{h,m} - u_{j} = \frac{4}{3}(U_{h,m}-U_{\frac{h}{2},m}) + O(\tau + h^4) Uh,m​−uj​=34​(Uh,m​−U2h​,m​)+O(τ+h4)
U h , m − u j = 4 3 ( U h , m − U h 2 , m ) + O ( τ ) + C h 4 , U h 2 , m − u j = 4 3 ( U h 2 , m − U h 4 , m ) + O ( τ ) + C ( h 2 ) 2 , ∥ U h , m − U h 2 , m ∥ ∞ ∥ U h 2 , m − U h 4 , m ∥ ∞ ≈ 4. U_{h,m} - u_{j} = \frac{4}{3}(U_{h,m}-U_{\frac{h}{2},m}) + O(\tau) + Ch^4,\\ U_{\frac{h}{2},m} - u_{j} = \frac{4}{3}(U_{\frac{h}{2},m}-U_{\frac{h}{4},m}) + O(\tau) + C(\frac{h}{2})^2,\\ \frac{\parallel U_{h,m} - U_{\frac{h}{2},m}\parallel_{\infty}}{\parallel U_{\frac{h}{2},m} - U_{\frac{h}{4},m}\parallel_{\infty}}\approx 4. Uh,m​−uj​=34​(Uh,m​−U2h​,m​)+O(τ)+Ch4,U2h​,m​−uj​=34​(U2h​,m​−U4h​,m​)+O(τ)+C(2h​)2,∥U2h​,m​−U4h​,m​∥∞​∥Uh,m​−U2h​,m​∥∞​​≈4.
根据后验误差估计，我们可以判断这种网格剖分的合理性。尤其是当微分方程的精确解未知的时候，根据后验误差估计公式的第一个式子，我们可以根据加细剖分的方式来得到差分格式近似解与未知精确解的相似程度。在代码运行的过程中，计算的第三个式子来确定网格剖分的合理性。\par
这里进行一个小说明，由于显格式稳定有条件，当我们取 δ x = h = 0.1 \delta x = h = 0.1 δx=h=0.1，那么根据a,b,c的选取，我们必须要取空间方向的分化 τ ≤ 0.005 \tau \leq 0.005 τ≤0.005，在显格式具体做法中，我们取得就是 τ = 0.005 \tau = 0.005 τ=0.005。因此显格式难以进行后验误差验证，在改变空间步长 h h h的同时必须要改变时间步长 τ \tau τ，因此在后面的微分方程求解过程中，要特别注意事先就取一个相对小的空间步长，不然无法比较后验误差。
u t − u x = u x x + u + e − 0.5 x s i n ( 5 x ) , x ∈ Ω , t ∈ ( 0 , 2 ] , u_t - u_x = u_{xx} + u + e^{-0.5x}sin(5x),x \in \Omega,t \in (0,2], ut​−ux​=uxx​+u+e−0.5xsin(5x),x∈Ω,t∈(0,2],
u ( x , 0 ) = 0 , x ∈ Ω , u ( 0 , t ) = 0 , t > 0 , u ( π , t ) = 0 , t > 0. u(x,0) = 0,x \in \Omega, u(0,t) = 0,t > 0, u(\pi,t) = 0,t > 0. u(x,0)=0,x∈Ω,u(0,t)=0,t>0,u(π,t)=0,t>0.
针对显格式
我们取最原始的分化长度 h = 0.1 ∗ π , τ = 0.0025 4 h = 0.1*\pi,\tau = \frac{0.0025}{4} h=0.1∗π,τ=40.0025​，可以证明，对于下面做后验误差的4种剖分组合，格式都满足最大值定理。
( h , τ ) (h,\tau) (h,τ) & ( h 2 , τ ) (\frac{h}{2},\tau) (2h​,τ) & ( h 4 , τ ) (\frac{h}{4},\tau) (4h​,τ) & ( h 8 , τ ) (\frac{h}{8},\tau) (8h​,τ)
E h , h 2 , h 4 = ∥ U h , m − U h 2 , m ∥ ∞ ∥ U h 2 , m − U h 4 , m ∥ ∞ E_{h,\frac{h}{2},\frac{h}{4}} = \frac{\parallel U_{h,m} - U_{\frac{h}{2},m}\parallel_{\infty}}{\parallel U_{\frac{h}{2},m} - U_{\frac{h}{4},m}\parallel_{\infty}} Eh,2h​,4h​​=∥U2h​,m​−U4h​,m​∥∞​∥Uh,m​−U2h​,m​∥∞​​
\
E h 2 , h 4 , h 8 = ∥ U h 2 , m − U h 4 , m ∥ ∞ ∥ U h 4 , m − U h 8 , m ∥ ∞ E_{\frac{h}{2},\frac{h}{4},\frac{h}{8}} = \frac{\parallel U_{\frac{h}{2},m} - U_{\frac{h}{4},m}\parallel_{\infty}}{\parallel U_{\frac{h}{4},m} - U_{\frac{h}{8},m}\parallel_{\infty}} E2h​,4h​,8h​​=∥U4h​,m​−U8h​,m​∥∞​∥U2h​,m​−U4h​,m​∥∞​​



非常奇怪的是，对于隐格式，针对不同的剖分，我们生成了两个近似解图片，但是图像差异很大，不过针对同一种剖分，隐格式和显格式的图像几乎没有差距。
基于后验误差，没有解析精确解的代码如下：

# -*- coding: utf-8 -*-
"""
Created on Thu Nov 12 19:41:06 2020


"""

# -*- coding: utf-8 -*-
"""
Created on Sun Nov  8 09:38:27 2020


"""

import numpy as np
import matplotlib.pyplot as plt
import torch
import time
import torch.nn as nn
import torch.nn.functional as F
st = time.time()

def FF(X):
    return torch.exp(-0.5*X[:,0])*torch.sin(5*X[:,0])
class FD():
    def __init__(self,bound,hx,prob):#这里的prob是无用参数，仅仅是为了偷懒，所以把有试验函数的代码进行了简单修改
        self.prob = prob
        self.dim = 2
        self.hx = hx
        self.nx = [int((bound[0,1] - bound[0,0])/self.hx[0]) + 1,int((bound[1,1] - bound[1,0])/self.hx[1]) + 1]
        self.size = self.nx[0]*self.nx[1]
        self.X = torch.zeros(self.size,self.dim)
        for m in range(self.nx[1]):
            for j in range(self.nx[0]):
                self.X[m*self.nx[0] + j,0] = bound[0,0] + j*self.hx[0]
                self.X[m*self.nx[0] + j,1] = bound[1,0] + m*self.hx[1]
        self.u_init()
    def u_init(self):
        self.u = torch.zeros(self.size,1)
        for m in range(self.nx[1]):
            for j in range(self.nx[0]):
                X = self.X[m*self.nx[0] + j:m*self.nx[0] + j + 1,:]
                if m == 0 or j == 0 or j == self.nx[0] - 1:
                    self.u[m*self.nx[0] + j] = 0
    def solve(self,a,b,c):
        mu = b*self.hx[1]/(self.hx[0]**2)
        nu = 0.5*a*self.hx[1]/self.hx[0]
        for m in range(self.nx[1] - 1):
            for j in range(1,self.nx[0] - 1):
                X = self.X[m*self.nx[0] + j:m*self.nx[0] + j + 1,:]
                ind = m*self.nx[0] + j
                self.u[(m + 1)*self.nx[0] + j] = (mu + nu)*self.u[ind - 1] + (1 + c*self.hx[1] - 2*mu)*self.u[ind] + (mu - nu)*self.u[ind + 1]\
                + self.hx[1]*FF(self.X[(m + 1)*self.nx[0] + j:(m + 1)*self.nx[0] + j + 1,:])
        return self.u
     
              
bound = torch.tensor([[0,np.pi],[0,2]]).float()

prob = 1

a = -1
b = 1
c = 1

#--------------------------------------------------
#下面是做后验误差分析过程
def host_error(hx,h):
    fd = FD(bound,hx,prob)
    ffd = FD(bound,h,prob)
    J = int(hx[0]/h[0]);M = int(hx[1]/h[1])
    
    fd_u = fd.solve(a,b,c)
    ffd_u = ffd.solve(a,b,c)
    u = torch.zeros(fd.size,1)
    
    for m in range(fd.nx[1]):
        for j in range(fd.nx[0]):
            u[m*fd.nx[0] + j] = ffd_u[M*m*ffd.nx[0] + J*j]
    return max(abs(u - fd_u))
#------------------------------------
h = 0.1*np.pi
t = 0.0025/4
ex = [h,t]
print("ex:%.3f"%(t/(h**2)))
ed = FD(bound,ex,prob)

fx = [0.5*h,t]
fd = FD(bound,fx,prob)
print("fx:%.3f"%(t/((0.5*h)**2)))
x_train = np.linspace(bound[0,0],bound[0,1],fd.nx[0])
y_train = np.linspace(bound[1,0],bound[1,1],fd.nx[1])
x,y = np.meshgrid(x_train,y_train)
u_pred = fd.solve(a,b,c)
plt.contourf(x,y,u_pred.reshape(fd.nx[0],fd.nx[1]).t(),40,cmap = 'Blues')
plt.colorbar()
plt.xlabel('x')
plt.ylabel('t')
plt.title('the FD solution')
plt.savefig('xian.jpg')

gx = [0.25*h,t]
gd = FD(bound,gx,prob)
print("gx:%.3f"%(t/((0.25*h)**2)))

hx = [0.125*h,t]
hd = FD(bound,hx,prob)
print("gx:%.3f"%(t/((0.125*h)**2)))

err = host_error(ex,fx)
er = host_error(fx,gx)
e = host_error(gx,hx)
print(err,er)
print(err/er)
print(er/e)

然后是隐格式

# -*- coding: utf-8 -*-
"""
Created on Thu Nov 12 19:41:06 2020


"""

# -*- coding: utf-8 -*-
"""

"""

import numpy as np
import matplotlib.pyplot as plt
import torch
import time
import torch.nn as nn
import torch.nn.functional as F
st = time.time()

def FF(X):
    return torch.exp(-0.5*X[:,0])*torch.sin(5*X[:,0])
class FD():
    def __init__(self,bound,hx,prob):#这里的prob是无用参数，仅仅是为了偷懒，所以把有试验函数的代码进行了简单修改
        self.prob = prob
        self.dim = 2
        self.hx = hx
        self.nx = [int((bound[0,1] - bound[0,0])/self.hx[0]) + 1,int((bound[1,1] - bound[1,0])/self.hx[1]) + 1]
        self.size = self.nx[0]*self.nx[1]
        self.X = torch.zeros(self.size,self.dim)
        for m in range(self.nx[1]):
            for j in range(self.nx[0]):
                self.X[m*self.nx[0] + j,0] = bound[0,0] + j*self.hx[0]
                self.X[m*self.nx[0] + j,1] = bound[1,0] + m*self.hx[1]
        
        
    def matrix(self,a,b,c):
        mu = b*self.hx[1]/(self.hx[0]**2)
        nu = 0.5*a*self.hx[1]/self.hx[0]
        A = torch.zeros(self.nx[0],self.nx[0])
        B = torch.zeros(self.nx[0],self.nx[0])
        for i in range(1,self.nx[0] - 1):
            A[i,i] = 1 - c*self.hx[1] + 2*mu
            A[i,i + 1] = - mu
            A[i,i - 1] = - mu
            B[i,i + 1] = - nu
            B[i,i - 1] = nu
            B[i,i] = 1
        A[0,0] = 1;A[-1,-1] = 1
        return A,B
    def right(self,a,b,c,m):
        r = torch.zeros(self.nx[0],1)
        for j in range(1,self.nx[0] - 1):
            X = self.X[(m + 1)*self.nx[0] + j:(m + 1)*self.nx[0] + j + 1,:]
            r[j] = self.hx[1]*FF(X)
        r[0] = 0
        r[-1] = 0
        return r
    def solve(self,a,b,c):
        A,B = self.matrix(a,b,c)
        u = torch.zeros(self.size,1)
        
        u[:self.nx[0]] = 0
        #print(u[:self.nx[0]].shape, UU(self.X[:self.nx[0],:],[0,0],self.prob).shape)
        for m in range(self.nx[1] - 1):
            rig = [e-mail][m*self.nx[0]:(m + 1)*self.nx[0]] + self.right(a,b,c,m)
            u[(m + 1)*self.nx[0]:(m + 2)*self.nx[0]],lu = torch.solve(rig,A)
        return u
     
              
bound = torch.tensor([[0,np.pi],[0,2]]).float()
prob = 1
a = -1
b = 1
c = 1

#--------------------------------------------------
#下面是做后验误差分析过程
def host_error(hx,h):
    fd = FD(bound,hx,prob)
    ffd = FD(bound,h,prob)
    J = int(hx[0]/h[0]);M = int(hx[1]/h[1])
    
    fd_u = fd.solve(a,b,c)
    ffd_u = ffd.solve(a,b,c)
    u = torch.zeros(fd.size,1)
    
    for m in range(fd.nx[1]):
        for j in range(fd.nx[0]):
            u[m*fd.nx[0] + j] = ffd_u[M*m*ffd.nx[0] + J*j]
    return max(abs(u - fd_u))
#------------------------------------
h = 0.1*np.pi
t = t = 0.05
ex = [h,t]
print("ex:%.3f"%(t/(h**2)))
ed = FD(bound,ex,prob)

fx = [0.5*h,t]
fd = FD(bound,fx,prob)
print("fx:%.3f"%(t/((0.5*h)**2)))
x_train = np.linspace(bound[0,0],bound[0,1],fd.nx[0])
y_train = np.linspace(bound[1,0],bound[1,1],fd.nx[1])
x,y = np.meshgrid(x_train,y_train)
u_pred = fd.solve(a,b,c)
plt.contourf(x,y,u_pred.reshape(fd.nx[0],fd.nx[1]).t(),40,cmap = 'Blues')
plt.colorbar()
plt.xlabel('x')
plt.ylabel('t')
plt.title('the FD solution')
plt.savefig('yin.jpg')
gx = [0.25*h,t]
gd = FD(bound,gx,prob)
print("gx:%.3f"%(t/((0.25*h)**2)))

hx = [0.125*h,t]
hd = FD(bound,hx,prob)
print("gx:%.3f"%(t/((0.125*h)**2)))

err = host_error(ex,fx)
er = host_error(fx,gx)
e = host_error(gx,hx)
print(err,er)
print(err/er)
print(er/e)

重点说明，之前我放了一个二维空间抛物型方程，里面列了一些差分格式求解，当时说显格式差分法不可取，原因是稳定性条件不满足，需要重新设计剖分长度。