# -*- coding: utf-8 -*-

import numpy as np
import numpy.linalg as npl
import matplotlib.pyplot as plt

##################################################################################################
# Résolution approchée de -u'' + alpha u = f dans ]0, 1[ avec conditions de Neumann homogènes    #
# par différences finies.                                                                        #
##################################################################################################

def laplacien_nh(J):
	d = 2.*np.ones(J)
	s = np.ones(J-1)
	A = np.diag(d) - np.diag(s,-1) - np.diag(s,1)
	A[0,0] = 1. # Approximation à l'ordre 1 de Neumann homogène. 
	A[J-1,J-1] = 1. # Approximation à l'ordre 1 de Neumann homogène. 
	#A[0,1] = -2. # Approximation à l'ordre 2 de Neumann homogène. 
	#A[J-1,J-2] = -2.# Approximation à l'ordre 2 de Neumann homogène. 
	return(A)

def f(x):
	#return(1 + 0.*x)
	return 8.*x*(1. - x) - 2.*x*x - 2*(1. - x)*(1. - x) + x*x*(1. - x)*(1. - x)

alpha = 1.
J = 100 # Nombre de points à l'intérieur de [0,1].
dx = 1/(J-1)
X = dx*np.arange(J)

Mat = 1/dx/dx*laplacien_nh(J) + alpha*np.eye(J)
print(Mat)
Mat[0,0] = 1./dx/dx
Mat[J-1,J-1] = 1./dx/dx
print(Mat)
F = f(X)
U = npl.solve(Mat,F)

plt.plot(X,U,label = 'Solution approchée')
plt.legend()
plt.grid('on')
plt.show()

print("Erreur L² : ", npl.norm(U - X*X*(1. - X)*(1. - X),2)/np.sqrt(J))
print("Erreur L^infty : ", npl.norm(U - X*X*(1. - X)*(1. - X),np.inf))


#############################