我正在尝试解决以下扩散方程问题：

$\frac{\partial f}{\partial t}=\frac{\partial (D\frac{\partial f}{\partial x})}{\partial x}+S$

$D=1+x^{2}+\sin(x)$

$f(x,0)=1 , f(0,t)=f(1,t)=1$

$S$被充分选择，以便拥有$f_{ana} = 1+\sin(\pi x)\sin(\pi t)$作为与数值方程进行比较的方程的解析解。

使用空间中的有限差分/中心方法（空间中的二阶）和时间隐含的欧拉（时间中的一阶）对该方程进行离散化。

全局错误行为应渐进地遵循$\dfrac{A}{dx^{2}}+\dfrac{B}{dt}$和$A$和$B$不依赖于$dt$和$dx$.

现在在 Python 上实现了解析算法，并根据以下方面构建了误差图变化$dx$ （第二规范和$\infty$norm ) 对于一个固定的$dt$，我期待以下行为$C+\dfrac{A}{dx^{2}}$，收敛速度“平稳”下降，因为 C 在较小的情况下将变得更加占主导地位$dx$直到达到常数。相反，我得到的（以及我不明白的）是收敛速度的小幅增加，然后恢复到恒定值（“倒钟”行为）。

在这种情况下，收敛速度“加速”是否正常？是否与非常小的空间误差的舍入效果有关？

这是现成的 Python 代码：

import numpy as np
import matplotlib.pyplot as plt
import scipy as sci
from scipy.sparse.linalg.dsolve import spsolve
import scipy.sparse as sp
import sympy as sym

#Symbolic calculation
x, t= sym.symbols('x t')
Dxx_1=1+x**2+sym.sin(x)

fi_1=1+sym.sin(sym.pi*x)*sym.sin(sym.pi*t) #Analytical solution
Soi_1=-sym.diff(Dxx_1*sym.diff(fi_1,x),x)+sym.diff(fi_1,t)#Source

Dxx_2=sym.lambdify((x),Dxx_1,"numpy")

#Conversion from symbolic to function
Soi=sym.lambdify((x,t),Soi_1,"numpy")
fi=sym.lambdify((x,t),fi_1,"numpy")    


#Lists that will contain the trucation errors for plot
Err00=[]
Err2=[]

#list of the number of points in the domain 
N=np.array([10,20,40,60,80,100,120,140,160,180,200,320,640])

#Loop on Nx
for Nx in N:

    #Space domain definition
    lx=1.0
    ox=0.0
    dx=(lx-ox)/(Nx)

    x1=np.linspace(ox+dx,lx-dx,Nx-1)#Points of the domain where f is unknown
    x2=np.linspace(ox,lx,Nx+1)#All points of the domain 

    #Diffusion coefficient
    Dxx=Dxx_2(x2)

    #Initialisation
    f=fi(x1,0)

    #Time grid
    t1=0.5
    dt=0.0001
    Nt=int(t1/dt)

    #Implicit matrix definition - Sparse method
    diag0= -(2*Dxx[1:-1]/(dx**2))
    diag1= ((Dxx[2:-1]-Dxx[:-3])/(4*dx**2))+((Dxx[1:-2])/(dx**2))
    diag_1= (-((Dxx[3:]-Dxx[1:-2])/(4*dx**2))+((Dxx[2:-1])/(dx**2)))
    data = [-dt*np.append(diag_1,[0]),1-dt*diag0,-dt*np.append([0],diag1)]     # Diagonal terms
    offsets = np.array([-1, 0, 1])                     # Their positions
    J = sp.dia_matrix((data, offsets), shape=(Nx-1,Nx-1))

    #Iteration and resolution
    for k in range(1,Nt+1,1):
        #Boundary condition vector
        S1=np.zeros(Nx-1)
        S1[0]=fi(0.0,k*dt)*(-(Dxx[2]-Dxx[0])/(4*dx**2)+(Dxx[1])/(dx**2))
        S1[-1]=fi(1.0,k*dt)*((Dxx[-1]-Dxx[-3])/(4*dx**2)+(Dxx[-2])/(dx**2))
        #2nd member vector
        S=Soi(x1,k*dt)
        f=spsolve(J,f+dt*S+dt*S1)
        sol_ana=fi(x1,k*dt)
        print('iteration '+str(k)+' - t = '+str(k*dt)+' - Implicit - Nx ='+str(Nx))


    print('Nx = '+str(Nx))
    Err2.append(100*np.linalg.norm(f-sol_ana)/np.linalg.norm(sol_ana))
    print('Err2 = ' + str(Err2[-1]))
    Err00.append(100*(abs(f-sol_ana).max()/abs(sol_ana).max()))
    print('Err00 = '+ str(Err00[-1]))


plt.figure()
plt.plot(N,Err00,'o',N,Err2,'o',N,N[0]**2*Err00[0]*N**float(-2),N,N[0]*Err00[0]*N**float(-1))
plt.title(' Err = f(Nx) - 1D - DF - Dirichlet  - Impl - dt='+str(dt)+' - t= '+str(t1))

plt.xscale('log')
plt.yscale('log')
plt.legend(['$Err_{\infty}$','$Err_{2}$','$\dfrac{1}{Nx^{2}}$','$\dfrac{1}{Nx}$'])
plt.xlabel('Nx')
plt.ylabel('Err')
plt.annotate('$f_{ana}='+sym.latex(fi_1)+'$', xy=(0.4, 0.14), xycoords='axes fraction')
plt.annotate('$D_{xx}='+sym.latex(Dxx_1)+'$', xy=(0.4, 0.1), xycoords='axes fraction')