%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt

plt.rc('figure',figsize=(14,9))

from IPython.html.widgets import interactive
from IPython.display import display

Résolution de l'équation de transport :

$$\frac{\partial u}{\partial t} - c \frac{\partial u}{\partial x} = 0\quad\hbox{dans }\mathbb{R} \times (0,T)\; .$$

avec la donnée initiale $g$

def g(x):
    y=np.where(x<0.0,1.0,-1.0)
    return y

Mise en place des schémas :

T=0.50
c=0.5
N=1000
J=1000
x=np.linspace(-1.0,1.0,J+1)

dx=2.0/J
dt=T/N
ll=c*dt/dx

print ll

u0=g(x)

plt.plot(x,u0,color='green')
plt.xlim([-1.0,1.0])
plt.ylim([-2.0,2.0])

b=np.matrix(np.zeros((J-1,1)))
b[0,0]=+ll*0.5*g(-1.0)
b[J-2,0]=-ll*0.5*g(1.0)

0.125

La donnée initiale : dans ce type d'équations, des données initiales discontinues sont classiques

Schéma explicite :

def schExpB(u0,ll,TY,J,N):
    u=np.copy(u0)
    v=np.copy(u0)

    if TY==1:
        for n in range(1,N):
            for j in range(1,J-1):
                v[j]=u[j]-ll*(u[j+1]-u[j])
            u=np.copy(v)    
    else:
        for n in range(1,N):
            for j in range(1,J-1):
                v[j]=u[j]-ll*(u[j]-u[j-1])
            u=np.copy(v)
        
    return u

def visu(N):
    u=schExpB(u0,ll,1,J,N)
    x=np.linspace(-1.0,1.0,J+1)
    plt.plot(x,u,color='red',)
    z=g(x-0.5*T)
    plt.plot(x,z,color='green')
    plt.xlim([-1.0,1.0])
    plt.ylim([-10.0,10.0])  

wdgt=interactive(visu,N=(0,50,5))
display(wdgt)

Schéma instable : la solution calculée est en rouge : apparition des premières oscillations qui vont "exploser", c'est-à-dire devenir de plus en plus grandes. Ici nous sommes dans les premiers pas de temps...

u=schExpB(u0,ll,0,J,N)
x=np.linspace(-1.0,1.0,J+1)
plt.plot(x,u,color='red',)
z=g(x-0.5*T)
plt.plot(x,z,color='green')
plt.xlim([-1.0,1.0])
plt.ylim([-2.0,2.0])

(-2.0, 2.0)

Schéma upwind "classique"

Un schéma vedette : Lax-Friedrichs

def schLF(u0,ll,TY,J,N):
    u=np.copy(u0)
    v=np.copy(u0)
    for n in range(1,N):
        for j in range(1,J-1):
#            v[j]=u[j]-0.5*ll*(u[j+1]-u[j-1])+0.5*ll*ll*(u[j+1]+u[j-1]-2.0*u[j])
            v[j]=(u[j+1]+u[j-1])*0.5-0.5*ll*(u[j+1]-u[j-1])
        u=np.copy(v)
    return u

u=schLF(u0,ll,0,J,N)
x=np.linspace(-1.0,1.0,J+1)
plt.plot(x,u,color='red',)
z=g(x-0.5*T)
plt.plot(x,z,color='green')
plt.xlim([-1.0,1.0])
plt.ylim([-2.0,2.0])

(-2.0, 2.0)

Schéma dissipatif : $|g|<1$ et donc on perd au niveau de l'amplitude de la solution...

Schéma implicite :

A=np.eye(J-1)
for i in range(J-2):
    A[i,i+1]=ll*0.5
    A[i+1,i]=-ll*0.5

def schImpB(u0,ll,A,J,N,b):
    v=np.matrix(np.zeros((J-1,1)))
    u=np.matrix(np.zeros((J-1,1)))
    for i in range(J-1):
        u[i,0]=u0[i+1]
    for n in range(1,N+1):
        v=np.linalg.solve(A,u+b)
        u=np.copy(v)
    return u

u=schImpB(u0,ll,A,J,N,b)
uu=np.copy(u0)
for i in range(1,J):
    uu[i]=u[i-1,0]
x=np.linspace(-1.0,1.0,J+1)
plt.plot(x,uu,color='red',)
z=g(x-0.5*T)
plt.plot(x,z,color='green')
plt.xlim([-1.0,1.0])
plt.ylim([-1.0,2.0])

(-1.0, 2.0)

La stabilité n'empêche pas quelques oscillations d'apparaître...