Résolution d’une EDP de type parabolique

2
∂u ∂ u
=α 2
Soit à résoudre l’équation de la chaleur ∂t ∂x avec α ∈ R fixe. Considérons le
système:

∂u k ∂2 u
{ = ¿ ¿
∂t cρ ∂x 2 ¿¿
¿
3
Prendre k=0.13 , c=0.11 ,g ρ=7.8 g/cm
Δx=0.25 , Δt sera donnée par les
,
k . Δt 1
r= ≺
2 2
conditions de stabilité définies par l’inéquation: cρΔx . La solution analytique est
donnée par:
∞
1 π (2 n+1)(x−1) −0 . 3738( 2 n+1)2 t
uexact ( x , t )=800 x ∑ 2 2
x cos xe
n=0 π (2 n+1) 2
1 1
r≻ r⊲
a) Utiliser la méthode directe (méthode explicite) en supposant 2 puis 2
Que constatez-vous?
b) Utiliser maintenant la méthode de Crank- Nicolson et faites varier le r. Que constatez-
vous?
c) Comparaisons des méthodes:

 Comparer, en fonction de r, graphiquement

uexact (0.25,t ) et
uapp ( 0.25,t ) puis
uexact (0.75,t ) et
uapp ( 0.75,t ) t
 Même question pour puis
uexact ( x ,099) et uapp ( x ,0.99) puis
uexact ( x ,1.98) et
uapp ( x ,1.98) .
d) Conclure.

1. Discrétisation de l’espace et du temps:

a) Discrétisation de l’espace:

Nous devons résoudre notre équation aux dérivées partielles dans l’espace compris

entre 0 et 2 avec un pas Δx . Pour passer donc trouver le point

xi à partir du
x
point i−1 , nous devons ajouter à ce dernier le terme Δx , de sorte que nous
avons:
 x i=x i−1 + Δx

⇒
{
x i−1=x i−2 + Δx
:
x 2 =xi + Δx
x1=x 0 + Δx
En sommant membres `a membres les équations ainsi obtenues et après
simplifications, nous trouvons:
x i=x 0 +(i−1+1) Δx=x0 +iΔx or
x 0 n’est autre que le point a de

l’intervalle [ a,b ] , donc, finalement:

x i=a+iΔx 1≤i≤n. x
b−a 2
n x= =
Avec Δx Δx

b) Discrétisation du temps:
Notre équation, c’est à dire le processus de résolution doit se refaire dans un intervalle

de temps égal à Δt jusqu’ à arriver à

T max . Pour passer de
tj à
t j+1 nous

devons à chaque fois ajouter Δt à

t j . Par un raisonnement par récurrence
comme précédemment, nous tirons:
t j = jΔt car et correspond à l’instant t=0
T max
t j = jΔt 1≤ j≤nt nt =
, et avec Δt

2. Programmation de la solution analytique:

∞
1 π (2 n+1)(x−1) −0 . 3738( 2 n+1)2 t
uexact ( x , t )=800 x ∑ 2 2
x cos xe
Soit: n=0 π (2 n+1) 2
Puisque nous sommes en analyse numérique et non en analyse fondamentale nous devons
savoir que nous ne travaillons qu’avec des points d’un maillage et non dans tout le domaine

comme on aurait cru le faire. Pour évaluer la valeur numérique de exact pour tout point u
appartenant dans la grille, nous devons faire une boucle qui, pour chaque point (xi; tj), calcule
uexact ( x i ,t j
. Posons ∞=100 acceptable en analyse numerique, nous avons donc le
programme suivant qui calcule les éléments de la matrice v représentant exactement les

valeurs de
u( x i ,t j ) dans le quadrillage.

***********************************
clc; clear;
k=0.13;c=0.11; p=7.8;dx=0.25;
r=1/4
dt=dx*dx*c*p*r/k;
Tmax=100*dt;%%(Nous avons fait Tmax un multiple de dt pour qu’en le divisant par dt, nt
soit toujours un entier, car sinon le programme ne compilera pas)%
cla=0; clb=0;
a=0;b=2;
nx=(b-a)/dx;
nt=Tmax/dt;
x=0:dx:b; t=0:dt:Tmax;
% *********** la solution analytique ********************
v=zeros(nx+1,nt+1);
n=0;
while(n · 100)
for i=1:nx+1
for j=1:nt+1

2 2 2
u(i,j)=v(i,j)+ 800∗1/(pi ∗(é∗n+1) *cos(pi∗(2∗n+1)∗(x(i)−1)/2)∗¿¿exp(−3.738∗(2n+1) ∗t(j))
v(i,j)=u(i,j);
end
end
n=n+1;
end
mesh(t,x,v)
*****************************
1
r= .
Voici la courbe que nous avons obtenue pour 4
3. Discrétisation de l’équation différentielle:

a) Rappel de la notation spatiale et temporelle:

Si la fonction u prend comme variables le temps

tj et de l’esapce

x i , y i . .. . , par commodité du langage, on notera

u( x i ,t j ) par
uJi et en

dimension 2D de l’espace et 1D en temps,

u( x i , y j ,t k ) sera notée par
uki, j .
 b) Discrétisation des dérivées partielles:
2
∂y
D’après les approximations, nous pouvons approximer ∂ x2 par
j
∂2
u ui+1−2 uij +ui+1
j
∂u
(x , t )≈
2 i j 2
∂x ( Δx) puisque ∂t
est une dérivée simple c’est à dire
une dérivée première, nous utilisons ici l’une des trois formules vues en projet.
j+1 j
∂u u i −ui
( xi ,t j )≈¿ ¿
Prenons la formule avale (droite): ∂t ( Δt )
Nous avons remplacé
h x par Δx et ht par Δt .

c) Equation aux dérivées partielles discrétisée:

Remplaçons 3.1 et 3.2 dans l’EDP donnée en énoncé. On aura:

j +1 j j j j
u i −ui n ui +1−2 ui +ui −1 Δt k
≈ x 2
⇒uij+1−u ij = 2
j
x x (u i+1 −2u ij + ui−1
j
)
( Δt ) cρ ( Δx) ( Δx ) cρ

r=
Δt k
x
Posons : ( Δx)2 cρ
l’équation devient: i i+1u j+1 =ru j +(1−2 r)u j +ru j
i i−1
Ceci est l’équation discrétisée de l’équation de la chaleur avec la méthode des différences
finies. On l’appelle équation de la méthode directe.
Remarque 4. :
a) r regroupe les constantes physiques du problème et celles de la discrétisation de

l’espace a,b et du temps

[ ] [ 0,T max ] . Donc les pas de discrétisations influencent
l’équation discrétisée.

b) Pour j=0 , t=0 et i i+1 u1=ru0 +(1−2 r)u0 +ru0

i i−1 l’´equation nécessite la
connaissance de la condition initiale pour démarrer le processus. Cette condition

n’est autre que

f (x , 0)=f ( x i ) donnée par l’énoncée du problème.

d) Résolution numérique du système:

 Quand on fixe j et on varie i de
1 à n x−1 , on obtient le système linéaire

de
n x−1 équations à
n x−1 inconnues suivant:

u 1j+1=ru 2j +(1−2 r )u1j +ru 0j

{
Avec:
u2j+1 =ru 3j +(1−2 r )u2j +ru 1j

j+ 1
x
j

u0j
:
:
j j
u n −1 =ru nx +(1−2 r )u2 +run x−2

la condition aux limites en

j+1
⇔ U =MxU +N
{
j

0≤ j≤nt −1
a , on la notera dans la programmation
j j
u0j =cla et
un x
la condition aux limites en b qui sera notée
un x =clb
.
Nous avons donc:

rU 0j

)( )
1−2r r 0 ....... ... 0 0

(
r 1−2r r ... : :
0 r 1−2r r : : N= :
: : : : : 0 :
0
: : ; r 1−2r r j
0 : ; 0 r 1−2r rU n
x

et
La matrice M est une matrice tri diagonale ou trigonale et la résolution est d’autant plus
simple à programmer avec une des méthodes itératives ou directes vues en analyse numérique
II. Mais, encore une fois, nous allons nous suffire de la résolution directe par la division
matricielle. C’est à dire si on a à résoudre le système Au=B alors on écrira en Matlab
x= A \ B et nous aurons automatiquement la réponse.
Attention: B doit être un vecteur colonne de même nombre d’éléments que le nombre des
lignes de A et le vecteur ⃗u sera donné en un vecteur colonne. Il ne nous est pas
nouveau de donner la transposée d’une matrice pour avoir tout ce que nous voulons.
1 1 1
r= r≺ r≻
4 pour le cas ou 2 et r=0. 625 pour le cas ou 2 ,
Nous avons pris
on a le programme suivant qui nous calcule les valeurs du vecteur ⃗u , remplacé par le
vecteur
⃗h dans le programme pour pouvoir lier le programme avec le premier calculant les
éléments de la matrice v pour la solution exacte, à chaque instant j et qui les stocke
dans une matrice w , en commanant par les valeurs initiales qui ne sont d’autres que les
conditions initiales. Pour avoir la courbe du deuxième cas, il suffit de remplacer r par
0.625.
clc; clear;
k=0.13;c=0.11;
p=7.8;dx=0.25;
r=1/4; dt=dx*dx*c*p*r/k;
Tmax=100*dt;
a=0;b=2;
cla=0;clb=0;
nx=(b-a)/dx;
nt=Tmax/dt;
x=0:dx:b; t=0:dt:Tmax;
for i=1:nx-1
N(i)=0;
end
N(1)=r*cla;
N(nx-1)=r*clb;
for i=1:nx-2
M(i,i)=1-2*r;
M(i,i+1)=r;
M(i+1,i)=r;
end
M(nx-1,nx-1)=1-2*r;
for i=1:nx+1
if x(i) < 1
Ci(i)=100*x(i);
else
Ci(i)=100*(2-x(i));
end
end
for i=1:nx-1
h(i)=Ci(i+1);
end
j=1;
h=h’;
while(j < nt + 2)
for i=1:nx-1
w(i,j)=h(i);
end
h=M*h+N’;
j=j+1;
end
for i=nx:-1:2
for j=nt+1:-1:1
w(i,j)=w(i-1,j);
end
end
for j=1:nt+1
w(1,j)=0;
w(nx+1,j)=0;
end
mesh(t,x,w);

Voici les courbes obtenues après compilation pour les deux cas:

Constatations:
1 1
r≺ r≻
Nous avons remarqué, d’après les deux graphes correspondants à 2 et 2 , que
1
r≺
la méthode présente deux régions de stabilité. Si 2
la méthode est stable et converge
vers la solution analytique. Nous avons une courbe semblable à celle de la solution
analytiquement établie. Plus j s’accroit, plus les valeurs du vecteur u deviennent
nulles. Ce qui est en entier accord avec la solution exacte qui est une série des fonctions
convergente, alors son terme général doit tendre forcement vers 0.
1
r≻
Mais si 2 , la méthode présente des failles et la courbe ne converge même pas.
L’allure obtenue dans ce cas est totalement en désaccord avec celle de la solution ´exacte.
Pour avoir un bon résultat dans cette méthode, il faut donc tenir compte de la valeur de r
car elle influencera le résultat.

4. Méthode de Crank - Nicolson ou méthode implicite:

En mathématiques, en analyse numérique, la méthode de Crank-Nicolson est un algorithme

simple permettant de résoudre des systèmes d’équations aux dérivées partielles. Cette
méthode utilise les différences finies pour approcher une solution du problème : elle est
numériquement stable, et quadratique pour le temps. On peut facilement la généraliser à des
problèmes à deux ou trois dimensions.
Cette méthode, publiée en 1947, est le résultat des travaux de la mathématicienne britannique
Phyllis Nicolson(1917 ¡ 1968) et du physicien John Crank (1916 ¡ 2006). Ils l’utilisèrent
dans la résolution de l’équation de la chaleur. Mais, peu après la méthode est devenue usuelle
dans plusieurs problèmes numériques.
 Methode:

Pour assurer une résolution numérique stable et ce dans toute la région de l’´etude, Crank et
Nicolson ont propose de prendre la moyenne des dérivées partielles spatiales de u entre les
t
instants j et
t j+1 au lieu de considérer seulement la dérivée spatiale à l’instant
tj ,
c’est `a dire:
α ∂2 u ∂2 u
∂u
( x ,t )=
[
( x , t )+
∂t i j 2 ∂ x 2 i j ∂ x 2 i j+1
(x ,t )
]
 Résolution numérique de l’EDP :

Nous remplaçons les dérivées partielles par leurs approximations déjà établies pour avoir
l’équation numérique suivante:

j j j+1 j+11 j+1 j+1 j

1 ui+1−2 ui ui+1 −2 ui + ui−1 cρ ui −ui
(
2 ( Δx )2
+
( Δx)2
= x
k Δt )
Δt k j
⇒2 (u ij+1−u ij )= (u −2u ij +u i−1
2 cρ i+1
j j+1
+ui+1 −2 uij+1 +ui+11
j+1
)
( Δx)
Δt k
r=
( Δx)2 cρ , on trouve l’équation:
Posons
j+1
−rui−1 +(2+2 r )uij+1−ru i+1
j+1 j
=ru i−1 +(2−2 r)uij +rui+11
j

Ceci est le schéma numérique de Crank - Nicolson.

Lorsque nous fixons j et faisons varier i , nous obtenons le système linéaire suivant:

−ru 0j+1 +(2+2 r )u1j+1 −ru 2j+1=ru0j +(2−2r )u 1j +ru 2j

{ −ru1j+1 +(2+2 r )u2j+1 −ru 3j+1 =ru1j +(2−2r )u 2j +ru 3j

−ru 2j+1 +(2+2r )u 3j+1 −ru 4j+1=ru 2j +(2−2 r )u3j +ru 4j
:
:
−run x − 2 +(2+2 r )un x−1−ru nx =runj x−1 +(2−2 r )unj x +runj x +1
j+1 j+1 j+ 1

⇒ M 1 U j+1 + N 1= M 2 U j +N 2

Avec :
 2+2 r −r 0 . . .. .. . .. 0 −rU 0j+1

( )( )
−r 2+2 r −r : : 0
0 −r 2+2r −r : : :
M 1= , N 1=
: : : : : 0 :
: : −r 2+2 r −r 0
0 .. . .. . 0 −r 2+2 r −rU nj+1
x

⇒U j+1=M −1 j
1 [ M 2 U +( N 2 −N 1 ) ]

Ce système matriciel est linéaire d’ordre

n x−1 .

 Programmation de la solution :

Voici le programme qui sert à donner la courbe de la solution numérique de la méthode de

Crank - Nicolson. Pour avoir la solution (la courbe) pour le deuxième cas, il suffit de
remplacer r par 0.25 :

clc; clear;
k=0.13;c=0.11;
p=7.8;dx=0.25;
r=0.625;
dt=dx*dx*c*p*r/k;
Tmax=100*dt;
a=0;b=2;
cla=0;clb=0;
nx=(b-a)/dx; nt=Tmax/dt;
x=0:dx:b; t=0:dt:Tmax;
for i=1:nx-1
N1(i)=0;
N2(i)=0;
end
N2(1)=-r*cla;
N2(nx-1)=-r*clb;
for i=1:nx-2
M1(i,i)=2+2*r;
M1(i,i+1)=-r;
M1(i+1,i)=-r;
M2(i,i)=2-2*r;
M2(i,i+1)=r;
M2(i+1,i)=r;
end
M1(nx-1,nx-1)=2+2*r;
M2(nx-1,nx-1)=2-2*r;
for i=1:nx+1
if x(i)<1
Ci(i)=100*x(i);
else
Ci(i)=100*(2-x(i));
end
end
for i=1:nx-1
h(i)=Ci(i+1);
end
j=1;
h=h’;
I=eye(nx-1)
while(j<nt+2)
for i=1:nx-1
w(i,j)=h(i);
end
h=(M1nI)*(M2*h+(N2’-N1’));
j=j+1;
end
for i=nx:-1:2
for j=nt+1:-1:1
w(i,j)=w(i-1,j);
end
end
for j=1:nt+1
w(1,j)=0;
w(nx+1,j)=0;
end
mesh(t,x,w);

Voici les courbes représentatives de w dans les deux situations:

Constatations
Nous remarquons ici que la méthode de Crank - Nicolson est inconditionnellement stable
 1
r≻
( c’est `a dire quel que soit e ∈ R ). Toutefois, quand 2
la méthode converge plus
vite vers 0 que si c’est le cas contraire. Il faut aussi savoir qu’on paie cette stabilité par la

résolution, à chaque instant

tj , d’un système linéaire qui nécessite une inversion de la
matrice [M1] tri diagonale, soit:
j+1 j j+1 −1 j
M1U +N 1 =M 2 U ° N 2 ⇒ U =M 1 [ M 2 U +( N 2 −N 1 )]

2−2r r 0 .. .. . . .. . 0 rU 0j

( ) ()
r 2−2r r : : 0
0 r 2−2 r r : :
M 2= , N 2= :
: : : : : 0 :
: : r 2−2 r r 0
0 . .. ... 0 r 2−2 r rU j nx

a. Comparaison des méthodes.

En fonction de r , nous voulons comparer graphiquement la méthode numérique explicite

(directe) avec la solution analytique pour voir comment la solution approchée se rapproche à
la solution exacte.

 Comparaison de uexact (0.25,t ) et uapp ( 0.25,t ) , uexact (0.75,t ) et

uapp ( 0.75,t )
Nous combinons le programme de la solution analytique avec celui de la solution numerique
directe pour pouvoir faire la comparaison, puis on ajoute à la fin du nouveau programme, la
partie suivante:
Remarquons que x=0 . 25 i=2 et x=0 .75 pour i=4 .

Sa1=v(2,:);
Sn1=w(2,:);
plot(t,Sa1,’r’,t,Sn1,’g’);
Sa2=v(4,:);
Sn2=w(4,:)
figure(2)
plot(t,Sa2,’r’,t,Sn2,’g’);

Voici donc les courbes comparatives de ces deux méthodes aux points 0.25 et 0.75:
Figure 2.6:
 Comparaison de uexact ( x ,0.99) et uapp ( x ,0.99) , uexact ( x ,1.98) et
uapp ( x ,1.98)

1
r= t j =0.99 t j =1.98
Remarquons qu’avec 4 , on n’a pas un pointj tel que ni

mais on voit que pour j=10 ;

t j =0.9281≈0.99 et pour j=20 ;
t j =0.9504≈1.98 . On refait ensuite la même chose mais on ajoute à présent la partie
suivante à la fin de la combinaison:

Sa1=v(:,10);
Sn1=w(:,10);
plot(t,Sa1,’r’,t,Sn1,’g’);
Sa2=v(:,20);
Sn2=w(:,20)
figure(2) plot(t,Sa2,’r’,t,Sn2,’g’);

Voici donc les courbes comparatives des solutions pour les deux méthodes à
t j =0.9281
et
t j =1.9594 .
b. Etudes théoriques de la méthode explicite :

 Convergence numérique :

La notion de Convergence signifie que la solution numérique (donnée par le schéma

numérique) approche la solution exacte (analytique) quand x →0 et Δt →0 , ce qui
signifie que la maille devient très petite.

Considérons
uij la solution numérique approchée et
U ij la solution ´exacte de
l’équation :
∂U k ∂2 U
=
∂ t cρ ∂ x 2 :
k Δt
j j j r= x
au point
x=x i et
t=t
j notons e i =U i −u i pour cp ( Δx )2 ,
uij+1 =r(ui+1
j j
+ui−1 )+(1−2r )u ij , e ij =U ij−u ij ⇒ uij =U ij −eij que l’on remplace
dans le schéma numérique précédente :
U ij+1 −e ij+1 =r(U i+1
j j
−ei+1 j
+U i−1 j
−e i−1 )+2r(U ij −e ij )
⇒ eij+1 =[ r( ei+1
j j
+e i−1 )+(1−2 r )e ij ] +U ij+1 −r(U i+1
j j
+U i−1 )−(1−2 r )U ij .
Avec la méthode de Taylor, nous trouvons:
∂U
U ij+1 =U ij−ei+1
j
+ Δ tx ( x , ξ ), t j ≤ξ≤t j+1
∂t i
2
∂U ( Δx ) ∂2 U
U ij+1 =U ij + Δx. ( xi , t j )+ . 2 (η1 ,t j ) x i ≤η1 ≤x i+1
∂x 2 ∂x
2
j ∂U ( Δx) ∂2 U
U i−1 =U ij + Δx . ( x , t )+ . 2 (η2 , t j ) x i−1 ≤η2 ≤x i
∂x i j 2 ∂x
que l’on remplace dans l’équation précédemment trouvée pour avoir:

∂U ∂U
e ij+1 =[ r(e i+1
j j
+ei−1 )+(1−2r )eij ]+U ij + Δt . (x i , ξ )−r (U ij +Δx . (x ,t )
∂t ∂x i j
( Δx)2 ∂2 U j ∂U ( Δx)2 ∂2 U
+ . 2 (η1 , t j ))−r (U i − Δx. ( x i ,t j )+ . 2 (η2 ,t j ))−(1−2r )U ij
2 ∂x ∂x 2 ∂x
Apres simplification des quantités en couleurs, l’équation devient:

∂U ( Δx)2 ∂ 2 U
e ij+1 =[ r (e i+1
j j
+ei−1 )+(1−2r )eij ]+Δt . ( x i , ξ )−r . 2 (η 1 ,t j )
∂t 2 ∂x
( Δx)2 ∂2 U
+ . 2 (η2 ,t j ))
2 ∂x

Dans un maillage très fin, les points de la discrétisation sont très proches les uns des autres.

Nous avons donc

x i≈η1≈x i+1≈η2≈ x i ainsi que
t j≈ξ≈t j+1 . Pour ce, nous

allons remplacer
η1 et
η2 par
xi et ξ par
t j puis r par
k Δt
cρ ( Δx )2 pour avoir l’équation sous la forme:
∂U k Δt ( Δx )2 ∂2 U
e ij+1 =[ r (e i+1
j j
+ei−1 )+(1−2r )eij ]+Δt . ( x i , t j )− .( ( x ,t )
∂t cρ ( Δx)2 2 ∂ x2 i j
( Δx)2 ∂ 2 U
+ . 2 ( x i ,t j ))
2 ∂x
2
Apres simplification de la quantité ( Δx ) , factorisons Δt :
2
∂U k ∂U
e ij+1 =[ r(e i+1
j j
+ei−1 )+(1−2r )eij ]+Δt .[ ( x i , t j )− ( x , t )]
∂t cρ ∂ x 2 i j
Le deuxième crochet, qui n’est autre que l’équation de la chaleur initialement posée dans
l’énoncée, est nul. Ce qui nous donc l’égalité finale:
e ij+1 =rei+1
j j
+rei−1 +(1−2 r )e ij
⇒|eij+1|=|rei+1
j j
+re i−1 +(1−2r)eij|

⇒ 1−2r|=1−2r,|
Si n< 1/2 , on a : 1-2r > 0  alors :
 j+1 j j j
|e i |≤r|ei+1|+r|e i−1|+(21−2r)|ei |

E j=max 1≤i≤n −1 |e i |
j
x

Posons :

⇒|e j+1|≤rE j + (1−2 r) E j ⇒|e j+1|≤E j ∀ 1≤ j≤n t −1

Par réccurence :

|e j+1|≤ E j ≤....E1≤E 0
E0 =U 0 −u 0 à t =0 ,
U 0 −u0 =0 , E0 =0,
Or, d’après les conditions initiales , donc
1
r≺
|e j+1|≤ 0 ⇒ e j =0 , ∀ j 2
et
1
r≺ .
Donc, la méthode explicite converge vers zéro quand 2

 Stabilité numérique :

On dit qu’un schéma numérique est stable si et seulement si au cours des calculs, l’erreur
commise d’une itération à l’autre (d’une maille à l’autre) n’infecte pas le calcul suivant.
Dans le cas contraire, le calcul peut exploser et on n’aura pas une bonne convergence. soit
l’équation de la chaleur:

2
∂U k ∂ U
−
∂ t cρ ∂ x 2
La numérisation de cette équation nous donne la relation matricielle:

j+1 j 1 0
U = AU ⇒U = AU
U0
etant donné par les conditions initiales.
 1 0 2 1 20 3 2 30
{U=AU¿{ =AU=AU¿{ =AU=AU¿{: ¿ U
0
par une certaine erreur et on donne Ū
0

Supposons qu’on introduit

1 0 2 1 20 3 2 30
{Ū=AŪ¿{ =AŪ=AŪ¿{ =AŪ=AŪ¿{: ¿
|e 0=Ū 0−U 0 ⇒ e j =Ū j−U j =A j U⃗ 0−A j U 0 =A j ( Ū 0−U 0 ) 
|⇒e j= A j e 0
Avec :
A admet N valeurs propres distinctes, donc ses vecteurs propres associés forment une
basse

A{ X1=λ1X¿{AX2=λ2X¿{AX3=λ3X {:¿ ¿
λi : est la valeur propre associée au vecteur propre
Xi : . Il existe alors

c i , c 2 , c 3 .. . .. c n tels que:
e 0 =c i X 1 +c 2 X 2 +. .. .+c n X n
N N N N
⇒ e = A ( ∑ ci X i )=∑ A c i X i =∑ c i A X i =∑ c i λij X i
j j j j

i=1 i=1 i=1 i=1

N
e j =∑ c i λij X i
i=1
N j
j
⇒|e |=|∑ c i λij X i|≤ ∑|c i|x|λ ij|x|X i|
i=1 i=1
L’erreur est maitrisable si
∀ |λi|≺1 . Ceci est vérifié si et seulement si:

max1≤i≤N |λi|<1
Pour la matrice A de cet exemple, les valeurs propres sont données.
iπ
λi =1−4 r sin2 , i=1 ,. . .. . , N
2(N +1)
Alors :
iπ k Δt
|λi|<1 ⇒|1−4 r sin 2 |<1 avec r= x
2(N +1) cρ Δx
iπ
⇒−1<1−4 r sin2 <1
2( N +1)

Pour :
2 iπ 1 sin2 iπ −1
−1<1−4 r sin ⇒ r< . [ ] , ∀ i=1,2 ,. .. . , N
2( N +1) 2 2( N +1)
1 sin 2 iπ −1
⇒r< . mini=1 ,. . .. , N [ ]
2 2( N +1 )
2 sin2 x sin 2 x −1
Or : ∀ x , sin x ≤1⇒ ≤1 , ∀ a>1 ⇒[ ] ≥1 , donc
a a
sin2 iπ −1 1
min i=1,.. .. , N [ ] ¿1 ⇒r<
2( N +1) 2
1
r<
On tire encore une fois que la méthode directe (explicite) est stable si 2 . Ce qui est en
parfait commun accord avec l’expérience réalisée.

Conclusion
Il est intéressant de remarquer que l’équation de la chaleur, introduite initialement pour d
´écrire la conduction thermique, apparait également dans d’autres branches de la physique
théorique. Elle permet par exemple de d´écrire :
 Le phénomène de diffusion ;

 Certains aspects probabilistes du mouvement brownien

Les méthodes de résolution de cette équation sont multiples et diversifiées. Toutefois, il faut
faire attention à la stabilité et la convergence de la méthode choisie car, simple et facile à
manipuler qu’elle soit, elle pourrait causer des instabilités inattendues et si on n’a pas la
solution analytique pour faire la comparaison, des confusions pourraient s’engendrer. Il serait
sage de payer le cout de résolution et avoir une méthode stable et convergente
inconditionnellement que de se contenter d’une résolution aléatoire ( valable pour quelques
valeurs particulières) vue qu’elle soit maniable facilement.