# -*- coding: utf-8 -*-
"""
Created on Mon Feb 02 17:56:19 2015

Listing 13.7 : équation des ondes, résolution par 
différences finies

@author: Ribot/Grivet

"""
import numpy as np

#  Vitesse
c = 1
#   Discrétisation en espace
xmin = 0.0; xmax = 4.0; npt = 81; nx = npt-2; h = (xmax-xmin)/(npt -1)
xx = np.linspace(xmin ,xmax, npt )
xx = xx.transpose()
xxint = xx[1:nx+1]
#    Matrice de discrétisation en espace 1D
K = 2*np.eye(nx,nx) - np.diag(np.ones(nx-1),1) - np.diag(np.ones(nx-1),-1)
K1D = K/h**2
meth = input('choix de la methode (chaine exp ou imp): ')
#   Discrétisation en temps
#   Cas explicite
if meth == 'exp':
    	Tfin = 1; dt = 0.9*h/c ; ntps = int(Tfin/dt)
#   Cas Implicite
elif meth == 'imp':
    	Tfin = 1; dt = 0.1; ntps = int(Tfin/dt)
u0=np.hstack((np.zeros_like(xxint[xxint <=1.5]), xxint[(1.5< xxint)&(xxint <=2)]-1.5,\
   2.5- xxint[(2<xxint)&(xxint <=2.5)], np.zeros_like(xxint[xxint>2.5])))
u1 = np.copy(u0)
for t in range(ntps):
    #  Euler Explicite
    if meth == 'exp':
        u=2*u1-c**2*dt**2*np.dot(K1D,u1)-u0
	#   Euler Implicite
    elif meth == 'imp':
        u= np.linalg.solve(np.eye(nx ,nx)+dt**2*c**2*K1D,2*u1-u0)
        u0 = u1 ; u1 = u;
uu = np.hstack((0,u,0)) 
matplotlib.pyplot.plot(xx,uu)