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.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')