Discretization Error amplification instead of stagnation to machine precision
I wrote a code on Python 2.7.5 to solve numerically the following differential equation.
* $\frac{\partial^2f}{\partial x^2}=-S$
* $S=\pi^{2}\sin(\pi x)$
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 Nx (10000,50000,100000,600000,900000), 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() 509