1
$\begingroup$

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

enter image description here

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.