import numpy as np
import matplotlib.pyplot as plt
Nx,L=100,1
def T0(Te,Tm,x):
    return Te+(Tm-Te)*np.sin(np.pi*x/L)
x=np.linspace(0,L,Nx+1)
plt.plot(x,T0(50,350,x))
plt.xlabel('x(m)')
plt.ylabel('T(x,0) en degre Celsius')
plt.grid()
plt.show()
rho,c,lamba,Te,Tm,Nx=2150,1000,1.65,50,400,100
D=lamba/rho/c
tau=L**2/D
dx=L/Nx
Nt=100000
dt=tau/Nt
r=D*dt/dx**2 #$ on doit vérifier que $r<1/2$ : c'est le critère de convergence 
tabT=np.zeros((Nx+1,Nt+1)) #$ tabT est un tableau de $Nx+1$ lignes et $Nt+1$ colonnes ne comportant que des zéros 
tabT[:,0]=T0(50,350,x) #$ désigne la première colonne du tableau 
tabT[0,:]=50 #$ désigne la première ligne du tableau 
tabT[Nx,:]=50 #$ désigne la dernière ligne du tableau 
for i in range(1,Nx):
    tabT[i,1]=tabT[i,0]+r*(tabT[i-1,0]-2*tabT[i,0]+tabT[i+1,0])
for j in range(1,Nt):
        for i in range(1,Nx):
                tabT[i,j+1]=tabT[i,j]+r*(tabT[i-1,j]-2*tabT[i,j]+tabT[i+1,j])
plt.plot(x,tabT[:,0])
for j in range(5):
    plt.plot(x,tabT[:,j*3000])
    plt.xlabel('x(m)')
    plt.ylabel('T(x,0)  en degre Celsius')
    plt.grid()
plt.show()   
plt.plot(np.linspace(0,tau,Nt+1),tabT[25,:])
plt.plot(np.linspace(0,tau,Nt+1),tabT[50,:])
plt.grid()
plt.show()