TP4 Résolution de l'équation de Burger

Question 1
Calculons la solution du problème d'évolution suivant:

u2
∂t u + u∂x ( ) = 0, x ∈ [−1, 2], t ∈ [0, T ]
2
u(−1, t) = 1, t ∈ [0, T ]
u(x, 0) = u0 (x), x ∈ [−1, 2]

où u0 est la condition initiale suivante:

⎧1 si x < 0
u0 (x) = ⎨1 − x si 0 ≤ x ≤ 1
⎩
0 si x > 1

Pour un champ donné u(x, t), et un coefficient de diffusion ν , la forme générale de

l'équation de Burger à une dimension de l'espace est le système suivant:

∂u ∂u ∂2u
+u =ν 2
∂t ∂x ∂x
Lorsque le terme de diffusion est nul, comme dans notre cas, ν = 0, on obtient une équation
de Burger sans viscosité.

∂u ∂u
+u =0 (forme advective)
∂t ∂x
Il s'agit d'un prototype de l'équation de conservation pouvant développer des discontinuités
(ondes de choc). Nous nous intéressons dans ce TP à sa forme conservatrice:

∂u 1 ∂u2
+ =0 (forme conservatrice)
∂t 2 ∂x
L'équation de Burger sans viscosité est plus généralement une équation hyperbolique non
linéaire du premier ordre.

La solution de l'équation peut être construite par la méthode des caractéristiques. Les
équations caractéristiques sont données par:
 dx du
= u et =0
dt dt
L'intégration de la deuxième équation nous permet de voir que u est constante le long de la
courbe caractéristique et que ces dernières sont des lignes droites.

u=c et x = ut + x0

où x0 est le point sur l'axe x du plan (x, t) à partir duquel la courbe caractéristique part.

Sachant que sur tout l'axe des abscisses, u est connue dès la condition initiale, et reste
invariante au fur et à mesure que l'on se déplace le long de la courbe caractéristique, on peut
écrire: u = c = u0 (x0 ) sur chaque courbe caractéristique.
La famille de trajectoires caractéristiques paramétrées par x0 est donnée par:

x = u0 (x0 )t + x0

Donc:

⎧1 ⋅ t + x 0 si x0 < 0 ⎧x0 + t si x < t

x(t) = ⎨(1 − x0) ⋅ t + x0 si 0 ≤ x0 ≤ 1 ⟹ x = ⎨(1 − t)x0 + t si t ≤ x ≤ 1
⎩ ⎩
0 ⋅ t + x0 si x0 > 1 x0 si x > 1

Après dérivation, on obtient:

⎧1 si x < t
x′ (t) = ⎨1 − x0 si t ≤ x ≤ 1
⎩
0 si x > 1

La solution de l'équation est donnée par:

u(x, t) = u0 (x0 ) = u0 (x − ut) avec x0 = x − u0 (x0 )t

Or,

⎧1 si x < t
u(x, t) = ⎨1 − x0 si 0 ≤ x ≤ 1
⎩
0 si x > 1

Pour 0 ≤x≤1,
 1−x
u(x, t) = 1 − x0 = 1 − (x − u(x, t)t) = 1 − x + u(x, t)t ⟹ u(x, t) =
1−t
Finalement,

⎧1 si x < t
u(x, t) = ⎨ 1−x si t ≤ x ≤ 1
⎩
1−t
0 si x > 1

Il s'agit d'une relation implicite qui détermine la solution de l'équation à condition que les
courbes caractéristiques ne se coupent pas.

En fonction de la condition initiale, si elles se croisent, alors il n'existe pas de solution unique,
il y a formation d'une onde de choc. Le temps de rupture ( breaking time tb ) avant qu'une
onde de choc puisse se former est donné par:

−1
tb = inf
x∈R u′0 (x)

On a :

⎧0 si x < 0
u′0 (x) = ⎨−1 si 0 ≤ x ≤ 1
⎩
0 si x > 1

Donc, tb = 1.

La solution obtenue précédemment par la méthode des caractéristiques est donc valable
pour t ≤ tb .

Afin d'établir une solution pour t > tb , il faut permettre les discontinuités de u. Cependant en
élargissant la classe des solutions, l'unicité de la solution n'est pls garantie.

Les relations de Rankine-Hugoniot permettent d'imposer des conditions de choc reliant les
sauts de u à travers la discontinuité.

Pour un système hyperbolique de la forme :

∂u ∂f (u)
+ =0
∂t ∂x
la vitesse de propagation d'un choc est donné par l'équation suivante:

Δf
s=
Δu
u2L u2
1 2 − 2R uL +uR
dans notre cas, f (u) = 2 u ainsi s= 2
uL −uR = 2
où uL et uR sont les valeurs de u à gauche et à droite de la discontinuité.

u0 (x0 ) convient comme solution, pour y < 1 on a comme valeur uL = 1 et pour y > 1 on
a uR = 0. On obtient après intégration de l'équation de discontinuité : u0 (t) = 12 (1 + t).

On en conclut que la solution discontinue, mais admissible au sens de Lax est:

Pour t <1:

⎧1 si x ≤ t
u(x, t) = ⎨ 1−x si t < x < 1
⎩
1−t
0 si x ≤ 1

et pour t ≥1:

u(x, t) = {
1 si x ≤ t+1
2
0 si x > t+1
2

Question 2
Dans un le fichier main.rs , nous définissons une fonction fluxnum qui calcule le flux
numérique pour le schéma de Rusanov, elle prend deux arguments uL et uR qui sont les
valeurs de u à gauche et à droite de la discontinuité.

Par définition du schéma de Rusanov, le flux numérique est donné par:

ρ2L + ρ2R λ
f (ρL , ρR ) = − (ρR − ρL ) ˋ λ = max(∣f ′ (ρL )∣, ∣f ′ (ρR )∣)
ou
2 2
u2
On a f (u) = 2 et f ′ (u) = u donc :
1 2 1
f (ρL , ρR ) = (ρL + ρ2R ) − max(∣ρL ∣, ∣ρR ∣)(ρR − ρL )
4 2
 fn fluxnum(ul: f64, ur: f64) -> f64 {
let abs = f64::abs;
let max = f64::max;
let lambda= max(abs(ul), abs(ur));
0.25*(ul*ul+ur*ur)-0.5*lambda*(ur-ul)
}

On definit ensuite une fonction sol_num qui calcule la solution numérique pour le schéma de
Rusanov.

fn sol_num(nx:usize, xmin:f64, xmax:f64, cfl:f64, vmax:f64, tmax:f64, mut t:f64 )-> Vec<f6

Elle prend en argument nx qui est le nombre de points de discrétisation, xmin et xmax qui
sont les bornes de l'intervalle, cfl qui est la condition CFL, vmax qui est la vitesse
maximale, tmax qui est le temps final et t qui est le temps courant.

Elle retourne un vecteur un qui contient les valeurs de la solution numérique.

Elle initialise dans un premier temps le pas de discrétisation h et le pas de temps dt en

fonction de la condition CFL.
Elle initialise également le vecteur xi qui contient les points de discrétisation ainsi que le
vecteur un et unp1 qui contiennent les valeurs de la solution numérique à l'instant n et
n+1 .

let h = (xmax -x min) / nx as f64;

let dt = cfl * h / vmax;
let xi: Vec<f64> = (0..(nx + 1)).map(|i| xmin + (i as f64) * h).collect();
let mut un: Vec<f64> = xi.iter().map(|x| sol_exacte(*x, t)).collect();
let mut unp1= un.clone();

Le calcul de la solution numérique se fait ensuite dans une boucle while qui s'arrête lorsque
le temps courant t est supérieur au temps final tmax .
 while t<tmax {
for i in 1..nx {
// schéma de Rusanov
unp1[i] += -dt/h*(fluxnum(un[i],un[i+1]) - fluxnum(un[i-1],un[i]));
}
t+=dt;
// conditions aux limites
unp1[0] = sol_exacte(xmin, t);
unp1[nx] = sol_exacte(xmax, t);
// mise à jour de un
un=unp1.clone();
}

Les conditions aux limites sont imposées grâce à la fonction sol_exacte qui retourne la
solution exacte pour un temps t et une position x .

fn sol_exacte(x: f64, t: f64) -> f64 {

if t <1.{
if x<=t{
return 1.;
}
else if x>=1.{
return 0.;
}else{
return (1.-x)/(1.-t);
}
}
if x<=(t+1.)/2.{
return 1.;
}
return 0.;
}

Nous souhaitons dans cette question comparer la solution numérique et la solution

approchée aux instants t = 1/2 , t = 1 et t = 2 avec des grilles à N = 50 , N = 100 ,
N = 1000 , N = 10000 points.

Pour ce faire dans le main.rs , nous initialisons les paramètres de notre simulation.
 let nx_tab: Vec<usize> = vec![50, 100, 1000, 10000];

let xmin = -1.;

let xmax= 2.;

let cfl = 0.8;

let vmax = 1.;

let tmax = 0.5;

let t = 0.;

Pour chaque valeur de nx dans nx_tab , on calcule la solution numérique, le pas de

discrétisation h et les vecteurs xi respectifs.

let un_tab: Vec<Vec<f64>> = nx_tab.iter().map(|nx| sol_num(*nx, xmin, xmax, cfl, vmax,

let h_tab: Vec<f64> = nx_tab.iter().map(|nx| (xmax-xmin)/(*nx as f64)).collect();

let xi_tab: Vec<Vec<f64>> = nx_tab.iter().zip(h_tab.iter()).map(|(nx, h)| (0..(*nx + 1

Afin de pouvoir tracer plusieurs graphes sur la même fenêtre, nous utilisons la bibliothèque
gnuplot .

// Création de la fenêtre
let mut fig = Figure::new();
// Création des axes
let mut ax = fig.axes2d();
// Titre de la fenêtre
ax.set_title("Solution numérique et exacte pour différentes valeurs de nx", &[]);
let color = ["light-blue", "dark-green", "blue", "pink"];
for i in 0..nx_tab.len() {
// Tracé des graphes pour chaque valeur de nx
ax.lines(&xi_tab[i], &un_tab[i], &[Caption(&format!("nx={}", nx_tab[i])), Color(co
}
// Tracé de la solution exacte
ax.lines(&xi_tab[3], &uex, &[Caption("Solution exacte"), Color("red")]);
fig.show();

Pour t=1/2 , on obtient les graphes suivants:

Pour t=1 , on obtient les graphes suivants:
Pour t=2 , on obtient les graphes suivants:
Le schéma de Rusanov est donc convergent.
On observe que plus le nombre de points de discrétisation est grand, plus la solution
approchée est proche de la solution exacte.

Question 3
Le code est similaire à celui du schéma de Rusanov, la seule différence réside dans la
fonction sol_num . On modifie l'expression de unp1 dans la boucle while comme suit:
 while t<tmax {
for i in 1..nx {
// schéma differences finies
unp1[i] += -dt/h*un[i-1]*(un[i]-un[i-1]);
}
t+=dt;
unp1[0] = sol_exacte(xmin, t);
unp1[nx] = sol_exacte(xmax, t);
un=unp1.clone();
}
return un;

Pour t=1/2 , on obtient les graphes suivants:

Pour t=1 , on obtient les graphes suivants:

Pour t=2 , on obtient les graphes suivants:
Tout comme le schéma de Rusanov, plus le nombre de points de discrétisation est grand,
plus le schéma aux différences finies est précis. Toutefois la convergence diminue avec le
temps, lorsque t=2s la solution approchée est décalée par rapport à la solution exacte. Ainsi,
on préférera le schéma aux différences finies lorsque t est petit les solutions approchées sont
plus précises.

Question 4
Nous souhaitons dans cette question résoudre le problème d'évolution suivant:

⎧ ∂t ρ + ∂x f (ρ) = 0 ∀x ∈]0, L[ et t ∈ [0, T ]

⎨ ρ(x, 0) = ρ0 (x)
⎩
∀x ∈]0, L[
ρ(0, t) = 0 ∀t ∈]0, T [

où ρ ∈ [0, ρ0 ] est la densité de véhicules et ρ0 (x) est la densité maximale de véhicules sur la
route. ρ = ρ0 correspond à une congestion totale de la route et ρ = 0 correspond à une
route vide.
 = 1 km. La densité et la
Afin de simuler ce trafic, on considère une route de longueur L
vitesse maximales sont respectivement ρ0 = 200 véhicules/km et v0 = 130 km/h.

La fonction f représentant le flux de véhicules est donnée par:

f (ρ) = ρv0 (1 − )
ρ
ρ0

On définit les fonctions suivantes, elle correspondent respectivement à la fonction f et à sa

dérivée:

fn f(rho : f64) -> f64 {

return rho * V_MAX * (1. - rho / RHO_MAX);
}

fn f_prime(rho : f64) -> f64 {

return V_MAX * (1. - 2. * rho / RHO_MAX);
}

fn flux_num_traffic(u1: f64, u2: f64) -> f64 {

let abs = f64::abs;
let max = f64::max;
let lambda = max(abs(f_prime(u1)), abs(f_prime(u2)));
return 0.5 * (f(u1) + f(u2)) - 0.5 * lambda * (u2 - u1);
}

On initialise dans le main également les éléments suivants en plus des variables déjà
définies:
 // le temps maximal de simulation
let tmaxmax = 5./60.;
// le pas de temps maximal auquel on observe l'évolution de la densité
let dtmax = 0.2/60.;

// vecteurs contenant les valeurs de la densité initiale

let mut un = vec![0.; nx as usize + 1];
let mut unp1 = vec![0.; nx as usize + 1];

// vecteur contenant l'état initial du bouchon

let mut uinit = vec![0.; nx as usize + 1];
for i in 0..nx as usize + 1 {
un[i] = rho_zero(xi[i]);
uinit[i] = rho_zero(xi[i]);
}

Nous pouvons créer une boucle while qui va itérer sur le pas de temps maximal tmaxmax . A
chaque itération, on va résoudre le problème d'évolution avec le schéma de Rusanov et on va
tracer la densité de véhicules en fonction de la position sur la route. On va également calculer
la vitesse maximale sur la route et on va l'afficher à chaque itération.
 while tmax < tmaxmax {
while t < tmax {
// calcul de la vitesse maximale
let mut vmax = 0.;
for i in 0..nx as usize + 1 {
let v= f_prime(un[i]);
if v > vmax {
vmax = v;
}
}
let dt = dx / vmax;
for i in 1..nx as usize {
// schéma de Rusanov
unp1[i] = un[i] - dt / dx * (flux_num_traffic(un[i], un[i+1]) - flux_num_t
}
unp1[0] = 0.;
unp1[nx] = 0.;
for i in 0..nx as usize + 1 {
un[i] = unp1[i];
}
t += dt;
}
plot1d(&xi, &uinit, &un);
tmax += dtmax;
}

Pour tmax = 0 secondes :

Pour tmax = 24 secondes :
`

Pour tmax = 72 secondes :

Pour tmax = 96 secondes :
Pour tmax = 108 secondes :
Pour tmax = 120 secondes :
C'est autour de tmax=108 secondes que le trafic automobile redevient fluide.