# -*- 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 Dirichlet homogènes  #
# par différences finies.                                                                        #
##################################################################################################

def laplacien_dh(J):
	d = 2.*np.ones(J)
	s = np.ones(J-1)
	A = np.diag(d) - np.diag(s,-1) - np.diag(s,1)
	return(A)

def f(x):
	return(1 + 0.*x)
	#return(4*np.pi*np.pi*np.sin(2*np.pi*x) + np.sin(2*np.pi*x)*np.sin(2*np.pi*x)) #XXX

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

Mat = 1/dx/dx*laplacien_dh(J) + alpha*np.eye(J)
F = f(X)
U = npl.solve(Mat,F)

print(U)

X = np.concatenate([[0],X,[1]]) # Ici on ajoute les points sur le bord pour l'illustration, 
#U = np.concatenate([[0],U,[0]]) # avec les conditions de Dirichlet homogènes. 
U = np.hstack((0,U,0))

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

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