Global truncation error behavior at fixed time step
I am trying to solve the following diffusion equation problem:
$\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 $ is choosen adequately in order to have $f_{ana} = 1+\sin(\pi x)\sin(\pi t)$ as the analytical solution of the equation to be compared with the numerical one.
The equation was discretized using a finite difference/ Central method in space (2nd order in space) and Euler implicit in time (1st order in time).
The global error behavior should follow asymptotically $\dfrac{A}{dx^{2}}+\dfrac{B}{dt}$ with $A$ and $B$ that do not depend on $dt$ and $dx$.
Now after implementing the resolution algorithm on Python, and constructing the error plot variation in terms of $dx$ (2nd norm and $\infty $ norm ) for a fixed $dt$, I was expecting a behavior following $C+\dfrac{A}{dx^{2}}$, with a "smooth" decrease in the convergence rate as C would become more predominant for smaller $dx$ till reaching the constant. What I got instead (and what I don't understand) is a small increase in the convergence rate, before recovering to the constant value ("reversed bell" behavior).
Is it normal to have an "acceleration" of the convergence rate in this case ? Is it a problem related to rounding effect of the very small space error ?
Here is the ready-to-use Python code :
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. 935