Introduction à Scilab pour le Calcul Scientiﬁque

Alexandre Ern
(last modiﬁcation date: October 10, 2017)

Version pdf de ce document
Version sans bandeaux

1 Introduction
2 Premier script Scilab : programmer un schéma
3 Deuxième script Scilab : variations sur la condition initiale et le schéma numérique
4 Troisième script Scilab : animation graphique
5 Quatrième script Scilab : un schéma avec inversion matricielle

1 Introduction

On considère le problème modèle

( { ∂u- ∂u- ∂t (x,t) + a ∂x(x, t) = 0, x ∈ ℝ,t ≥ 0, ( u(x,0) = u0(x), x ∈ ℝ.

(1)

L’inconnue est la fonction u qui dépend de la coordonnée spatiale x et du temps t. Les données du problème sont la fonction u₀ appelée condition initiale et le paramètre réel a appelé vitesse d’advection. On vériﬁe aisément que la solution de (1) est

u(x,t) = u0(x − at).

(2)

Une interprétation physique de l’équation (1) est la suivante : on considère un ﬂeuve rectiligne représenté par la droite réelle et qui s’écoule à la vitesse a. La quantité u(x,t) représente la concentration au point x et à l’instant t d’un produit polluant qui a été déversé accidentellement dans la ﬂeuve à t = 0 selon le proﬁl u₀. Le polluant est transporté par l’écoulement et à un temps t ﬁxé, la fonction xu(x,t) modélise la répartition de polluant le long du ﬂeuve.

Pour simpliﬁer, nous allons supposer que la donnée initiale u₀ est périodique (en espace) de période L si bien que nous pouvons nous restreindre à l’intervalle x ∈ [0,L] et chercher la solution u satisfaisant

( |{ ∂∂ut(x,t) + a∂∂ux(x,t) = 0, x ∈ [0,L ],t ≥ 0, u(x,0) = u (x), x ∈ [0,L ], |( 0 u(L,t) = u(0,t), t ≥ 0.

(3)

La dernière équation est une condition limite de périodicité.

Aﬁn d’approcher numériquement la solution u, on se donne

un entier N permettant de déﬁnir un pas d’espace δx = et un maillage uniforme de l’intervalle [0,L] constitué des N points x_i = (i − 1)δx, 1 ≤ i ≤ N;
un pas de temps δt permettant de construire la suite d’instants discrets tⁿ = nδt pour n ≥ 0.

Notre objectif est d’évaluer des quantités U_iⁿ, 1 ≤ i ≤ N et n ≥ 0, telles que

U ni ≃ u (xi,tn).

L’idée est d’approcher les dérivées partielles dans (3) à l’aide de développements de Taylor et d’écrire par exemple

{ n+1 n ∂u n U i ≃ Ui + δt∂t(xi,t), U n ≃ U n+ δx∂u(x ,tn). i+1 i ∂x i

On obtient

Uin+1−--Uni- Uin−-U-ni−-1 δt + a δx = 0.

(4)

En réarrangeant (4), il vient

n+1 n n n Ui = Ui − ν (Ui − U i− 1),

(5)

où nous avons posé

aδt ν = δx-.

(6)

Ce paramètre, sans dimension physique et connu sous le nom de nombre de Courant, jouera un rôle déterminant par la suite.

L’équation (5) déﬁnit un schéma numérique qu’on appelera schéma upwind (la terminologie, inspirée de l’anglais, rappelle le fait que la dérivée en espace a été décalée dans la direction d’où provient le vent (ou le courant)). Notons que (5) est un schéma explicite : étant données les valeurs de U_iⁿ pour 0 ≤ i ≤ N, (5) permet de calculer explicitement les valeurs de U_iⁿ⁺¹ pour 0 ≤ i ≤ N. La simulation numérique s’arrête lorsque le temps courant tⁿ dépasse un temps physique de simulation choisi à l’avance.

Il nous reste à approcher la condition initiale et la condition limite en écrivant

{ U 0= u (x), 0 ≤ i ≤ N, in 0n i U1 = UN , n ≥ 0.

2 Premier script Scilab : programmer un schéma

On écrira un script Scilab qu’on appelera script1.sce et qui comprendra 4 parties :

initialisation des paramètres physiques : vitesse de propagation, longueur du domaine de calcul, temps maximum de simulation et condition initiale. Pour cette dernière, on prendra une fonction en créneau
${ 1 si 1 < x < 1.5, u0(x) = 0 sinon.$ (7)
initialisation des paramètres numériques; on choisira N = 201 et on en déduira δx puis on choisira le nombre de Courant ν = 0.8 et on en déduira le pas de temps δt.
boucle en temps; initialiser (U_i⁰)_0≤i≤N puis évaluer (U_iⁿ⁺¹)_0≤i≤N en fonction de (U_iⁿ)_0≤i≤N tant que tⁿ ≤ T.
visualisation graphique des résultats; aﬃcher sur une même fenêtre la solution exacte et la solution approchée.

Conseils Scilab.

initialiser le script par la commande clear; qui nettoye l’espace mémoire occupé avant le lancement du script.
paramètres physiques :
  a=0.1;// vitesse d'advection
  Tmax=10;// temps maximum de simulation
  L=5;// longueur du domaine

Pour la condition initiale, on utilisera la fonction bool2s dont on consultera l’aide en ligne (help → Scilab Programming → bool2s). La condition initiale sera programmée sous forme de fonction Scilab, ce qui devrait donner quelque chose comme

  function v=condinit(x)
    v=bool2s((x > 1. ) & (x < 1.5));
  endfunction

On notera que v est une variable locale à la fonction condinit.
paramètres numériques :
  N=201;// nb de points de discrétisation en espace
  dx=L/(N-1);// pas d'espace
  nu=0.8;// nombre de Courant-Friedrichs-Levy
  dt=nu*dx/a;// pas de temps
initialisation du maillage et de la donnée initiale :
x=linspace(0,L,N)';
u=condinit(x);

(voir l’aide en ligne; la transposition permet de manipuler des vecteurs colonne ce qui est plus pratique pour la sortie graphique).
boucle en temps : à chaque temps discret tⁿ, on eﬀectue le calcul (5) de façon vectorielle en tenant compte de la condition de périodicité
  tps=dt:dt:Tmax;
  for t=tps do
    uold=u;
    u(2:N)=uold(2:N)-nu*(uold(2:N)-uold(1:N-1));
    u(1)=u(N);
  end
  t=tps($);// $ pointe sur le dernier élément du vecteur
sortie graphique : voir l’aide en ligne Help → Graphic Library. On pourra par exemple utiliser les commandes suivantes :
  xselect()
  xbasc()
  xtitle('temps='+string(t)+', dt='+string(dt)+', NU='+string(nu));
  plot2d([x,x],[u,uexact],style = [5,2],strf = "111",leg = "sol. numerique@sol. exacte", ...
         rect = [0,-1.5,L,1.5]);
pour exécuter le script, taper dans la fenêtre Scilab la commande exec script1.sce.
aﬁn de mesurer le temps d’exécution du programme, on insérera la commande
timer();

avant le démarrage de la boucle en temps et la commande

timer()

en dernière ligne du script Scilab. Le temps d’exécution écoulé entre les 2 appels de la fonction timer sera aﬃché dans la fenêtre Scilab.

Question 1 Augmenter progressivement le paramètre ν et observer les résultats. Quelle est la valeur critique?

3 Deuxième script Scilab : variations sur la condition initiale et le schéma numérique

Outre la condition initiale (7), on souhaite également pouvoir considérer la condition initiale

u (x) = sin (8πx∕L ). 0

(8)

De plus, on souhaite avoir le choix entre 2 schémas numériques : le schéma upwind (5) et le schéma de Lax-Wendroﬀ

n+1 n ν- n n ν2- n n n U i = U i − 2 (Ui+1 − Ui−1) + 2 (Ui+1 − 2U i + Ui−1),

(9)

dont la justiﬁcation sera étudiée en cours de Calcul Scientiﬁque.

Ecrire un script Scilab qu’on appelera script2.sce et qui oﬀre ces 2 options.

Conseils Scilab.

on introduira deux nouvelles variables qu’on pourra par exemple appeler choixCI et choixSCH et qui prendront la valeur 1 ou 2 en fonction de l’option retenue. On utilisera la commande select de Scilab dont la syntaxe est
  function v=condinit(x)
    select choixCI
      case 1 then
       v=bool2s((x > 1. ) & (x < 1.5));
      case 2 then
       v=sin(8*%pi*x/L);
    end
  endfunction
à chaque temps discret tⁿ, on eﬀectue le calcul (9) de façon vectorielle en tenant compte de la condition de périodicité (il y a maintenant deux cas particuliers à prendre en compte, car le schéma peut “déborder” sur la gauche ou sur la droite) :
  u(2:N-1)=uold(2:N-1)-nu/2*(uold(3:N)-uold(1:N-2))+ ...
           nu^2/2*(uold(3:N)-2*uold(2:N-1)+uold(1:N-2));
  u(N)=uold(N)-nu/2*(uold(2)-uold(N-1))+nu^2/2*(uold(2)-2*uold(N)+uold(N-1));
  u(1)=u(N);

Question 2 Repartir de ν = 0.8 et essayer les 4 possibilités du couple condition initiale / schéma numérique. Quelles conclusions tirez-vous ? Augmenter progressivement ν pour le schéma de Lax-Wendroﬀ et observer.

4 Troisième script Scilab : animation graphique

L’objectif de cette nouvelle étape est de pouvoir visualer la solution numérique au ﬁl des itérations au lieu de visualiser uniquement sa valeur ﬁnale.

Partant du script script2.sce, on réalisera un nouveau script qu’on appelera script3.sce.

Conseils Scilab.

avant la boucle en temps, il convient pour réaliser une animation d’initialiser la fenêtre graphique de la façon suivante :
  driver("X11")
  xselect()
  xbasc()
  xset("pixmap",1)

La commande driver("X11") permet d’accélérer l’aﬃchage en ne stockant pas en mémoire toutes les données relatives au tracé de la courbe mais seulement les plus importantes (perte de la possibilité du zoom par exemple). La commande xset("pixmap",1) évitera le clignotement de la fenêtre lors de la boucle en temps.
au sein de la boucle en temps, la fenêtre graphique est rafraichie de la façon suivante :
  xset("wwpc")
  uexact=a_ecrire;
  xtitle(a_ecrire);
  plot2d(a_ecrire);
  xset("wshow")

La commande xset("wwpc") eﬀace le contenu courant de la fenêtre graphique en mémoire. La commande xset("wshow") aﬃche dans la fenêtre graphique le contenu courant de la mémoire.
à l’issue de la boucle en temps, on revient à la conﬁguration initiale de la fenêtre graphique
xset("pixmap",0)
driver("Rec")

Question 3 Reprendre les essais précédents puis augmenter progressivement le paramètre ν et observer (pour le schéma de Lax-Wendroﬀ, prendre ν = 1.2, ν = 1.5 et ν = 1.8). Quelles conclusions en tirer?

5 Quatrième script Scilab : un schéma avec inversion matricielle

Nous allons modiﬁer le schéma de Lax-Wendroﬀ (9) en “implicitant le terme de diﬀusion” de la façon suivante

n+1 n ν- n n ν2- n+1 n+1 n+1 Ui = Ui − 2(U i+1 − U i−1) + 2 (U i+1 − 2U i + Ui−1 ).

(10)

En introduisant les vecteurs Wⁿ = (U_iⁿ − ν
2 (U_i+1ⁿ − U_i−1ⁿ))_{0≤i≤N−1} et Uⁿ⁺¹ = (U_iⁿ⁺¹)_{0≤i≤N−1} ainsi que la matrice A d’ordre N − 1 donnée par

( 1 + ν2 − ν2- 0 ... 0 − ν2 ) | 2 2 2 | | − ν2- ... ... 0 | || .. .. .. .. || A = || 0 . . . . || | ... ... ... ... 0 | || . . 2 || ( 0 .. .. − ν2 ) − ν2 0 ... 0 − ν2 1 + ν2 2 2

le schéma (10) s’écrit sous la forme

AU n+1 = W n.

La condition de périodicité est utilisée directement dans le schéma numérique aﬁn d’éliminer U_Nⁿ. La taille des vecteurs est donc N − 1 et non plus N comme dans les sections précédentes. Ce choix a été fait pour simpliﬁer la structure de la matrice A.

On constate que le schéma (10) n’est plus explicite mais implicite : si Uⁿ est connu, l’évaluation de Uⁿ⁺¹ nécessite l’inversion d’un système linéaire. On s’attend donc à ce que le coût d’une itération dans (10) soit plus élevé que dans (9). Cependant, comme nous le verrons dans les expériences numériques, le schéma implicite (10) n’est pas limité par la valeur du nombre de Courant ν ce qui permet de considérer des pas de temps plus grands et donc de réaliser moins d’itérations en temps. Il n’est donc pas évident de déterminer a priori qui du schéma explicite ou implicite sera le plus performant.

Aﬁn de réaliser les expériences numériques, on partira du script script3.sce dont on enlevera les parties relatives au schéma numérique. Le nouveau script sera appelé script4.sce.

Conseils Scilab.

initialisation de la matrice A : on utilisera la fonction diag de Scilab :
  aa=nu^2/2;
  A=diag((1+2*aa)*ones(1,dimA))-diag(aa*ones(1,dimA-1),-1)-diag(aa*ones(1,dimA-1),1);
  A(1,dimA)=-aa;
  A(dimA,1)=-aa;
initialisation de Wn :
  Wn(2:dimA-1)=u(2:dimA-1)-nu/2*(u(3:dimA)-u(1:dimA-2));
  Wn(1)=u(1)-nu/2*(u(2)-u(dimA));
  Wn(dimA)=u(dimA)-nu/2*(u(1)-u(dimA-1));
deux options sont possibles pour la boucle en temps : soit résoudre le système linéaire à chaque pas de temps
  u=A\Wn;//  (dans la boucle en temps)

soit inverser la matrice A avant les itérations en temps puis utiliser son inverse

  B=inv(A);
  tps=dt:dt:Tmax;
  for t=tps do
    // a_ecrire
    u=B*Wn;
  end

Question 4 Comparer les 2 options pour l’inversion du système linéaire. Réaliser des expériences numériques en faisant varier le nombre de Courant. Discuter des performances relatives des schémas explicite et implicite.

Introduction à Scilab pour le Calcul Scientiﬁque

Contents

1 Introduction

2 Premier script Scilab : programmer un schéma

3 Deuxième script Scilab : variations sur la condition initiale et le schéma numérique

4 Troisième script Scilab : animation graphique

5 Quatrième script Scilab : un schéma avec inversion matricielle