In [3]:
u = diffrea(0.1,0.5, 2, -np.pi,np.pi, 0.1, 0.01)
import numpy as np
import matplotlib.pyplot as plt
from scipy.sparse import linalg,diags
def diffrea(a, c, T, xmin, xmax, dx, dt):
""" calcule et trace la solution de l'edp
u_t - c u_xx - a u =0
avec la donnée initiale u(., t=0) =f(x)
pour f(x) de différentes formes.
sur xmin<x<xmax, 0<t<T
avec la méthode des différences finies
on code 2 schémas
discrétisation explicite centrée
discrétisation implicite centrée
on affiche la CFL parabolique = c*dt/dx^2
le script retourne la matrice u(j,n,s) = u(x_j, t_n,s)
pour les 2 schémas s=1 explicite, 2 implicite
"""
# maillage
x=np.arange(xmin,xmax+dx,dx)
Nx=np.size(x)
t=np.arange(0,T+dt,dt)
Nt=np.size(t)
#initialisation pour les 2 schémas
u=np.empty((Nx,Nt,2));
# condition initiale
def f(x):
#decommenter la fonction choisie.
# fonction Heaviside lissée
# y = 0.*(x<0)+4*x.^2.*(3-4*x).*(0<=x).*(x<=0.5)+ 1.*(x>0.5);
#fonction 1_(x<0)
# y = 1.*(x<0);
#fonction creneau 1_(-1/2,1/2)
# y = 1.*(-0.5<=x).*(x<1/2);
# cloche de Gauss.
# y = exp(-25*( x- (xmin + xmax)/2).^2);
# sinus
y=np.sin(x)
return y;
# donnee initiale pour les 2 schémas.
u[:,0,0] = f(x)
u[:,0,1] = f(x)
# CFL
r= c*dt/dx**2
s= a*dt
#boucle en temps
for n in range(Nt-1):
# on utilise le slicing plus efficace que les boucles
u[1:-1,n+1,0] = r*u[:-2,n,0] + (1-2*r+s)* u[1:-1,n,0] + r*u[2:,n,0]
#condition au bord x=xmin periodique
u[0,n+1,0] = r*u[-1,n,0] +(1-2*r +s )* u[0,n,0] + r*u[1,n,0]
#condition au bord x=xmax periodique
u[-1,n+1,0] = r*u[-2,n,0] +(1-2*r +s)* u[-1,n,0] + r*u[0,n,0]
####schema implicite###
#création de la matrice tridiagonale
#on utilise le format sparse (creux)
A = diags([-r,-r, 1+2*r-s, -r,-r],[-Nx+1,-1, 0, 1,Nx-1], shape=(Nx,Nx),
format="csr")
# A[0,-1]= -r et A[-1,0]= -r à cause des conditions limites périodiques
#
# boucle en temps
# a chaque itération on résout A U^{n+1} = U^n
for n in range(Nt-1):
u[:,n+1,1] = linalg.spsolve(A,u[:,n,1])
# tracé des solutions à différents instants
fig= plt.figure(figsize=(20, 10))
plt.subplot(121)
plt.plot(x, u[:,0,0], label='condition initiale')
plt.plot(x, u[:,-1,0], color='brown', label=f'$t={T}$')
plt.plot(x,np.exp((a-c)*T)*np.sin(x) , '*', label='Exact solution at ' f'$t={T}$')
plt.xlabel('$x$')
plt.ylabel('$u(.,t)$')
plt.title('schéma différences finies explicite')
plt.legend();
plt.subplot(122)
plt.plot(x, u[:,0,1], label='condition initiale')
plt.plot(x, u[:,-1,1], color='brown', label=f'$t={T}$')
plt.plot(x,np.exp((a-c)*T)*np.sin(x) , '*', label='Exact solution at ' f'$t={T}$')
plt.xlabel('$x$')
plt.ylabel('$u(.,t)$')
plt.title('schéma différences finies implicite')
plt.legend();
return u
u = diffrea(0.1,0.5, 2, -np.pi,np.pi, 0.1, 0.01)
CFL = 0.5*0.01/(0.1)**2
CFL
0.4999999999999999
u = diffrea(0.1,0.5, 2, -np.pi,np.pi, 0.1, 0.0102)
CFL = 0.5*0.0102/(0.1)**2
CFL
0.5099999999999999