import numpy as np
import matplotlib.pyplot as plt
from scipy.sparse import linalg,diags

def diffusion(mu, T, xmin, xmax, dx, dt):
    """ calcule et trace la solution de l'edp
         u_t - mu u_xx =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 3  schémas
         discrétisation explicite centrée
         discrétisation implicite centrée
        discrétisation de Crank-Nicolson
         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 3 schémas s=0 explicite,
         1 implicite 2 Crank-Nicolson
    """ 
    # maillage
    x=np.arange(xmin,xmax,dx)
    Nx=np.size(x)
    t=np.arange(0,T,dt)
    Nt=np.size(t)
    #initialisation pour les 3 schémas
    u=np.empty((Nx,Nt,3));
    # 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 = np.exp(-25*( x- (xmin + xmax)/2)**2);
         return y
    # donnee initiale pour les 4 schémas.
    u[:,0,0] = f(x)
    u[:,0,1] = f(x)
    u[:,0,2] = f(x)
    # CFL
    r= mu*dt/dx**2;
    ###schéma explicite
    #iterations 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)* 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)* 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)* 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, -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
    
    # iterations 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])
        
    ####schéma Crank-Nicolson####
    ####création des  matrices tridiagonales au format sparse
    A = diags([-r/2,-r/2, 1+r, -r/2,-r/2],[-Nx+1,-1,0,1,Nx-1],shape=(Nx,Nx),format="csr")
    # A[0,-1]= -r/2 et A[-1,0]= -r/2 à cause des conditions limites périodiques
    
    #
    B = diags([r/2,r/2, 1-r, r/2,r/2],[-Nx+1,-1,0,1,Nx-1],shape=(Nx,Nx), format="csr")
    #B[0,-1]= r/2 et B[-1,0]= r/2 à cause des conditions limites périodiques
    #iteration en temps
    #a chaque itération on résout A U^{n+1} = B U^n
    for n in range(Nt-1):
        u[:,n+1,2] = linalg.spsolve(A,B@u[:,n,2])

    
    # tracé des solutions à différents instants
    
    fig= plt.figure(figsize=(20, 10))
    #
    plt.subplot(131)
    plt.plot(x, u[:,0,0], label='condition initiale')
    plt.plot(x, u[:,int(T/(2*dt)),0], color='green', label=f'$t={T/2}$')
    plt.plot(x, u[:,-1,0], color='brown', label=f'$t={T}$')
    plt.xlabel('$x$')
    plt.ylabel('$u(.,t)$')
    plt.title('schéma différences finies explicite, CFL={:.2f}'.format(r))
    plt.legend();
    #
    plt.subplot(132)
    plt.plot(x, u[:,0,1], label='condition initiale')
    plt.plot(x, u[:,int(T/(2*dt)),1], color='green', label=f'$t={T/2}$')
    plt.plot(x, u[:,-1,1], color='brown', label=f'$t={T}$')
    plt.xlabel('$x$')
    plt.ylabel('$u(.,t)$')
    plt.title('schéma différences finies implicite, CFL={:.2f}'.format(r))
    plt.legend();
    #
    plt.subplot(133)
    plt.plot(x, u[:,0,2], label='condition initiale')
    plt.plot(x, u[:,int(T/(2*dt)),2], color='green', label=f'$t={T/2}$')
    plt.plot(x, u[:,-1,2], color='brown', label=f'$t={T}$')
    plt.xlabel('$x$')
    plt.ylabel('$u(.,t)$')
    plt.title('schéma différences finies Crank-Nicolson, CFL={:.2f}'.format(r))
    plt.legend();
    #
    return u

Résultats des différents schémas pour les paramètres μ = 0.5, T = 0.1, xmin = −1, xmax = 1 δx = 0.1 avec les pas de temps suivants δt = 0.005, δt = 0.01 et δt = 0.05.

u = diffusion(0.5,0.1, -1,1, 0.1, 0.005)
u = diffusion(0.5,0.1, -1,1, 0.1, 0.01)
u = diffusion(0.5,0.1, -1,1, 0.1, 0.05)

Pour le pas de temps $\delta t= 0.005, CFL=0.25$ on constate de bons résultats pour les 3 schémas. Pour $\delta t = 0.01, CFL=0.5$ le schéma explicite donne une solution irrégulière mais de valeur comparable aux deux autres schémas. Enfin pour $\delta t =0.05, CFL= 2.5$ le schéma explicite donne des résultats aberrants, tandis que le schéma implicite donne de bons résultats. Le schéma de Crank-Nicolson donne des résultats suspects près de zéro.

On change enfin μ en μ = −0.2. Que constatez-vous ?

u = diffusion(-0.2,0.1, -1,1, 0.1, 0.005)

Quand la diffusion est négative, tous les schémas donnent des résultats aberrants. Le problème mathématique est mal posé lorsque le coefficient de diffusion est négatif.