I wrote a code on Python 2.7.5 to solve numerically the following differential equation.
S is chosen that way in order to have $f= \sin(\pi x)$ as the exact solution, that can be compared with the constructed solution.
The differential operator is discretized using a second-order central finite difference method.
When running the code, the order of magnitude of the evolution of the discretization error
Err is correct (2nd order) for the first batch of points of
N (10,20,40,80,100,800,2000,5000) (goes down to 10e-11).
But I found that for very small
dx meaning, very big
Err doesn't stagnate at the machine precision error (as expected) but goes up for several magnitude orders, and even worse, the discretization error
Err deteriorates (up to 10e-7).
Is it a problem in the precision/round managment of Python for very small numbers, or is it normal to have even worse Error for that amount of points ?
I tried to write the problem in different ways (Matrix, Loops, Matrix inversion technices) but the same phenomenon happened.
import numpy as np import matplotlib.pyplot as plt import scipy.sparse as sp from scipy.sparse.linalg.dsolve import spsolve Err= N=np.array([10,20,40,80,100,800,2000,5000,10000,50000,100000,600000,900000]) #Loop on the number of points used for the discretization of the 1D domain for Nx in N: dx=(1.0)/Nx# fixed distance between two consecutive points in the domain x = np.linspace(dx,1.0-dx,Nx-1)# points of the 1D domain #contruction of the diff operator discretization Matrix data = [np.ones(Nx-1), -2*np.ones(Nx-1), np.ones(Nx-1)]# Diagonal terms offsets = np.array([-1, 0, 1])# Their positions LAP = sp.dia_matrix((data, offsets), shape=(Nx-1, Nx-1))#Sparse matrix S = np.pi**2*np.sin(np.pi*x)*dx**2# Second member f = spsolve(LAP,-S)#Resolution of the matrix system f_ana=np.sin(np.pi*x)#The exact solution Err.append((np.absolute((f_ana-f))).max()) plt.figure() plt.plot(N,Err,N,N**float(-2),'--') plt.legend(['Err','Theoretical asymptotic behaviour']) plt.xlabel('dx') plt.ylabel('Err') plt.xscale('log') plt.yscale('log') plt.show()