# -*- coding: utf-8 -*-
"""
Exercice 13.2f : résolution d'une équation différentielle par la
méthode des éléments finis. 
Created on Mon Sep 11 14:10:15 2017

@author: ribot/grivet
"""
import numpy as np
import numpy.linalg as lg
import matplotlib.pyplot as plt
# Paramètres
alfa = 10; epsilon = 1
# DiscrĂŠtisation en espace
xmin = 0; xmax = 1
npt = 21; nint = npt-1; nx = npt-2
h = (xmax-xmin)/nint
xx = np.linspace(xmin,xmax,npt)
#xx=xx'
xxint = xx[1:nx+1]
# Matrice de rigiditĂŠ
K = (2*np.eye(nx,nx)-np.diag(np.ones(nx-1),1)
-np.diag(np.ones(nx-1),-1))/h
#Matrice de masse
M = (4*np.eye(nx,nx)+np.diag(np.ones(nx-1),1)
+np.diag(np.ones(nx-1),-1))*h/6
# Second membre
F = xxint
# RÊsolution - système 1D
u = lg.solve(epsilon*K+alfa*M,F)
uu=np.hstack((0,u,0))
plt.plot(xx,uu)
plt.title("solution de $-u''(x)+10u(x) = x$, avec $u(0) = u(1) = 0$")