# -*- 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 non 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)) 

alpha = 1.
a = 1. # Valeur à gauche. 
b = 2. # Valeur à droite. 
J = 100 # 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)
F[0] = F[0] + a/dx/dx
F[J-1] = F[J-1] + b/dx/dx
print(F)
U = npl.solve(Mat,F)

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

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

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