import numpy as np
import matplotlib.pyplot as plt

def transport(c, T, xmin, xmax, dx, dt):
    """ calcule et trace la solution de l'edp
         u_t + c u_x =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 4  schémas explicites
         discrétisation  décentrée amont, centrée, Lax-Friedrichs
         et Lax-Wendroff
         on affiche la CFL  = c*dt/dx
        le script retourne la matrice u(j,n,s) = u(x_j, t_n,s)
         pour les 4 schémas s=1 décentré,
         2 centré 3 Lax-Friedrichs 4 Lax-Wendroff
    """
    # maillage
    x = np.arange(xmin, xmax + dx, dx)
    t = np.arange(0, T + dt, dt)
    # initialisation pour les 4 schémas
    u = np.empty((len(x), len(t), 4))
    # donnee initiale pour les 4 schémas
    # u(x,t,i) représente la solution du schéma i.
    # on pourrait economiser de la mémoire en ne
    # conservant que la valeur finale u(:,T,i) et
    # ecrasant les valeurs aux temps intermédiaires

    u[:, 0, 0] = f(x)
    u[:, 0, 1] = f(x)
    u[:, 0, 2] = f(x)
    u[:, 0, 3] = f(x)
    # CFL
    r = c * dt / dx
    # boucle en temps
    for n in range(len(t) - 1):
        # on utilise le SLICING pour éviter les boucles en espace
        # v[1:-1] = v[1],v[2], ...v[N-1]
        
        # schema decentre amont
        u[1:, n + 1, 0] = r * u[:-1, n, 0] + (1 - r) * u[1:, n, 0]
        # condition au bord x = xmin constant
        u[0, n + 1, 0] = u[0, n, 0]

        # schema centre
        u[1:-1, n + 1, 1] = (r / 2) * u[:-2, n, 1] + \
            u[1:-1, n, 1] - (r / 2) * u[2:, n, 1]
        # conditions au bord constantes
        u[0, n + 1, 1] = u[0, n, 1]
        u[-1, n + 1, 1] = u[-1, n, 1]
        # schema Lax-Friedrichs
        u[1:-1, n + 1, 2] = (1 + r) / 2 * u[:-2, n, 2] + \
            (1 - r) / 2 * u[2:, n, 2]
        # conditions au bord constantes
        u[0, n + 1, 2] = u[0, n, 2]
        u[-1, n + 1, 2] = u[-1, n, 2]
        # schema Lax-Wendroff
        u[1:-1, n + 1, 3] = (r ** 2 + r) / 2 * u[:-2, n, 3] + \
            (1 - r ** 2) * u[1:-1, n, 3] + (r ** 2 - r) / 2 * u[2:, n, 3]
        # conditions au bord constantes
        u[0, n + 1, 3] = u[0, n, 3]
        u[-1, n + 1, 3] = u[-1, n, 3]

    # tracé des courbes pour chaque schéma
    plt.close('all')
    fig = plt.figure(figsize=(10, 8))
    fig.suptitle(
        'Solution de l\'équation de transport pour différents schémas')

    # subplot 1
    plt.subplot(221)
    plt.plot(x, f(x), '-b', label='t=0')
    plt.plot(x, u[:, -1, 0], '-k', label='t=T')
    plt.plot(x, f(x - c * T), '--k', label='sol. exacte t=T')
    plt.xlabel('x')
    plt.ylabel('u(.,t)')
    plt.title('Schéma décentré amont, CFL={:.2f}'.format(r))
    plt.legend(loc='lower left')

    # subplot 2
    plt.subplot(222)
    plt.plot(x, f(x), '-b', label='t=0')
    plt.plot(x, u[:, -1, 1], '-k', label='t=T')
    plt.plot(x, f(x - c * T), '--k', label='sol. exacte t=T')
    plt.xlabel('x')
    plt.ylabel('u(.,t)')
    plt.title('Schéma centré, CFL={:.2f}'.format(r))
    plt.legend(loc='lower left')

    # subplot 3
    plt.subplot(223)
    plt.plot(x, f(x), '-b', label='t=0')
    plt.plot(x, u[:, -1, 2], '-k', label='t=T')
    plt.plot(x, f(x - c * T), '--k', label='sol. exacte t=T')
    plt.xlabel('x')
    plt.ylabel('u(.,t)')
    plt.title('Schéma Lax-Friedrichs, CFL={:.2f}'.format(r))
    plt.legend(loc='lower left')

    # subplot 4
    plt.subplot(224)
    plt.plot(x, f(x), '-b', label='t=0')
    plt.plot(x, u[:, -1, 3], '-k', label='t=T')
    plt.plot(x, f(x - c * T), '--k', label='sol. exacte t=T')
    plt.xlabel('x')
    plt.ylabel('u(.,t)')
    plt.title('Schéma Lax-Wendroff, CFL={:.2f}'.format(r))
    plt.legend(loc='lower left')

    plt.tight_layout()
    plt.show()

# donnée initiale
def f(x):
    # decommenter la fonction choisie.
    # fonction Heaviside lissée
    # y = np.zeros_like(x)
    # y[x >= 0] = 4 * x[x >= 0] ** 2 * (3 - 4 * x[x >= 0]) * (x[x >= 0] <= 0.5) + 1 * (x[x >= 0] > 0.5)
    # fonction 1_(x<0)
    y = np.zeros_like(x)
    y[x < 0] = 1
    # fonction creneau 1_(0,1/2)
    # y = np.zeros_like(x)
    # y[(x >= 0) & (x < 0.5)] = 1
    # cloche de Gauss.
    # y = np.exp(-25 * (x - (xmin + xmax) / 2) ** 2)
    return y

# Paramètres
c = 0.5
T = 0.75
xmin = -1
xmax = 1
dx = 0.01
dt = 0.01

transport(c, T, xmin, xmax, dx, dt)

Le schéma centré explose. Les schémas décentrés amont et de Lax-Friedrichs donnent des résultats similaires. Le schéma de Lax-Wendroff semble plus précis mais il présente des petites oscillations au voisinage de la discontinuité. Le schéma n'est pas monotone.

c = 0.5
T = 0.75
xmin = -1
xmax = 1
dx = 0.01
dt = 0.02
transport(c, T, xmin, xmax, dx, dt)

Les schémas donnent des bons résultats sauf le schéma centré qui explose.

c = 0.5
T = 1
xmin = -1
xmax = 1
dx = 0.01
dt = 0.023
transport(c, T, xmin, xmax, dx, dt)

Tous les schémas donnent de mauvais résultats. C'est normal car le nombre $ CFL: = \frac{c \delta t}{\delta x } >1$.

c = -0.5
T = 0.75
xmin = -1
xmax = 1
dx = 0.01
dt = 0.01
transport(c, T, xmin, xmax, dx, dt)

Lorsque $c<0$ le schéma décentré amont devient de fait décentré aval donc il explose encore plus violemment que le schéma centré. Les schémas de Lax-Friedrichs et Lax-Wendroff donnent de bons résultats, mais il y toujours des oscillations au voisinage de la discontinuité dans le schéma de Lax-Friedrichs.

PS Rappels sur le slicing

v=[3,4,6,7,11]

v[-1]

11

v[2:3]

[6]

v[1:-1]

[4, 6, 7]

len(v)

5

v[1:len(v)]

[4, 6, 7, 11]