TP: Résolution Numérique Des Edps Par DF: - Équation de La Chaleur
TP: Résolution Numérique Des Edps Par DF: - Équation de La Chaleur
TP: Résolution Numérique Des Edps Par DF: - Équation de La Chaleur
Considérons le problème de diffusion de la chaleur dans une barre homogène, de longueur L, de coefficient
de conduction λ, de coefficient calorifique Cp , sans production d’énergie interne. La température u(x, t)
au point x et à l’instant t vérifie l’équation suivante :
∂u ∂ 2u
(x, t) = a 2 (x, t), 0 ⩽ x ⩽ L, 0 ⩽ t ⩽ T
(P) u(x, ∂t ∂x
0) = u0 (x), 0 ⩽ x ⩽ L, condition initiale
u(0, t) = b0 (t), u(L, t) = bL (t), 0 ⩽ t ⩽ T, conditions aux bords.
λ
avec a = la diffusivité thermique.
ρCp
Nous allons chercher une solution approchée définie sur un maillage en espace et en temps. Pour
abdellatif.elghazi@uir.ac.ma
L
un entier N strictement positif donné, on pose h = = ∆x, c’est le pas d’espace et on choisit
N
T
k = ∆t = le pas de temps. Nous noterons ui j la solution approchée aux points xi = ih, i = 0, · · · , N
M
au temps tj = jk, j = 0, · · · , M.
Les conditions initiales et aux limites s’écrivent : t
ui j+1 − ui j = ruj+1
i−1 − 2rui
j+1
+ ruj+1
i+1 , i = 1, · · · , N − 1.
4. Montrer que la résolution du schéma d’Euler implicite est équivalente à la résolution du système
linéaire :
AUj+1 = Uj + f
où Uj = u1 j , u2 j , · · · , ujM−1 . Expliciter A et f.
1|??
22 novembre 2023
5. Montrer que la matrice A est inversible.
6. Application
La fonction suivante :
|??
22 novembre 2023
32 U (2: N ,j +1) = A \ U (2: N ,j ) ;
33 end ;
34 % Solution Exacte
35 n = length ( t ) ;
36 for j =1: n
37 uu (: , j ) = uexac (x , t ( j ) ) ;
38 end
39 e = norm ( U (: , n ) - uu (: , n ) , inf ) ;
40 end
1 clear all ;
2 clc ;
3 I =[0 ,1];
4 J =[0 ,0 .4 ];
5 N =10;
6 dt =0 .1 ;
7 u0 = @ ( x ) sin ( pi * x ) ; % condition initiale
8 uexac = @ (x , t ) sin ( pi * x ) . * exp ( - pi * pi * t ) ; % Solution exacte
9 [U ,x ,t , uu , e ]= chaleur (N , dt ,I ,J , u0 , uexac ) ;
.@..
10 fprintf ( 'l erreur obtenue pour % d est % f \ n ' ,N , e ) ;
11 % On trace la solution exacte et la solution approchee a l ' instant ...
finale
12 n = length ( t ) ;
13 plot (x , uu (: , n ) ',x , U (: , n ) ', '+ ')
14 legend ( ' Solution Exacte ' , ' Solution Approché ')
15 % % Code à complter
.
(c) Compléter le code avec les instructions pour afficher l’erreur pour d’autres valeurs de dt =
0.001 et N = 20.
(d) On prends N = 100. Compléter le code avec les instructions pour tracer la courbe de l’erreur
en fonction de dt.
(e) Pour dt = 0.01 et N = 50, le scripte suivant trace en 3D les courbes de la solution exacte
et la solution approchée à un instant t.
|??
22 novembre 2023
6 xlabel ( 'x ')
7 ylabel ( 't ')
8 title ( ' solution approchee ')
9 axis ([0 1 0 0 .4 0 1])
10 % Matrice avec les valeurs exactes
11 subplot (2 ,1 ,2)
12 surf ( xx , tt ,uu ')
13 xlabel ( 'x ')
14 ylabel ( 't ')
15 title ( ' solution exacte ')
16 axis ([0 1 0 0 .4 0 1])
.
Commenter les lignes de code et expliquer la valeur de sortie de chaque ligne.
(f) Pour étudier l’ordre du schéma on va calculer l’erreur commise entre la solution exacte et
la solution approchée quand le pas dt ∈ [0.01, 0.005, 0.0025, 0.0013, 0.0005, 00025] avec
N = 200.
On considère une membrane rectangulaire fixée sur les bords (On néglige les forces élastiques dues à la
déformation de la membrane) qui se déforme sous l’effet d’une force transversale f(x, y) par élément de
surface et on cherche la réflexion verticale u. Le problème modèle est une équation de poisson donnée
par :
−∆u(x, y) = f(x, y), ∀(x, y) ∈ Ω = [a, b] × [c, d]
(P 2 )
u(a, y) = u(b, y) = u(x, c) = u(x, d) = 0, Conditions aux limites
∂ 2u ∂ 2u
avec ∆u(x, y) = (x, y) + (x, y).
∂x2 ∂y2
Introduisons une partition régulière du segment [a, b] en Nx sous intervalles [xi , xi+1 ] de longueur
b−a
hx = avec x0 = a, xNx = b et xi = a + ihx pour i = 0, · · · , Nx
N
De même nous introduisons une partition régulière du segment [c, d] en Ny sous intervalles [yj , yj+1 ]
d−c
de longueur hy = avec y0 = c, yNy = d et yi = c + jhy pour j = 0, · · · , Ny .
Ny
4|??
22 novembre 2023
Schéma à 5 points pour le laplacien
Les noeuds de la grille sont les points Mi,j de coordonnées (xi , yj ). Ainsi, les inconnus sont les
valeurs uij qui approchent la solution exacte u(xi , yj ) aux noeuds intérieurs pour 1 ⩽ i ⩽ Nx − 1, et
1 ⩽ j ⩽ Ny − 1.
Il y a (Nx − 1) × (Ny − 1) inconnus.
1. Proposer un schéma basé sur une formule de différence finies centrées pour résoudre le problème(P 2 ).
2. Pour simplifier on prends hx = hy = h. Montrer que la résolution du schéma numérique est
équivalent à la résolution de
AU = f
Expliciter U, A et f.
3. On considère le problème :
−∆u(x, y) = 0, ∀(x, y) ∈ Ω = [0, 20] × [0, 10]
(P 3 )
u(0, y) = u(x, 0) = u(x, 10) = 0, et u(20, y) = 100 Conditions aux limites
5 5 5
On prends hx = hy = h ∈ {5, , , }
2 4 8
4. En variant h, résoudre le système (P 3 ). Donner à chaque fois le nombre d’inconnues.
5
5. pour h = 5 et h = construire manuellement la matrice du système. Dessiner une grille avec les
2
points inconnus.
6. Application : Le scripte suivant permet de résoudre l’equation de Laplace en 2D.
1 clc ; clear
2 h =0 .25 ;
3 a =0;
4 b =20;
5 c =0; d =10;
6 nx = floor (( b - a ) / h ) ;
7 ny = floor (( d - c ) / h ) ;
8 n =( nx -1) *( ny -1) ;
9 % %% ( remplissage des elements de la matrice A )
10 A = zeros ( n ) ;
|??
22 novembre 2023
11 for i =1:( n -1)
12 A (i , i ) = -4;
13 A ( i +1 , i ) =1; A (i , i +1) =1;
14 if ( mod (i ,( nx -1) ) ==0)
15 A ( i +1 , i ) =0;
16 A (i , i +1) =0;
17 end
18 end
19 for i =1: n - nx +1
20 A ( nx -1+ i , i ) =1;
21 A (i , nx -1+ i ) =1;
22 end
23 A (n , n ) = -4;
24 % %% ( remplissage des ´el´ements de la matrice B )
25 for i =1: n
26 B ( i ) =0;
27 if ( mod (i , nx -1) ==0)
28 B ( i ) = -100;
29 end
30 end
31 % %% ( resolution du systeme Av = B et transformation du vecteur ~v en la ...
matrice U ) .
32 V = A \B ';
33 k =1;
34 for j =1: ny -1
35 for i =1: nx -1
36 u (j , i ) = V ( k ) ;
37 k = k +1;
38 end
39 end
40 % %% ( d´ecallage des ´elements pour ins´erer les conditions aux limites )
41 for j = ny : -1:2
42 for i = nx : -1:2
43 u (j , i ) = u (( j -1) ,i -1) ;
44 end
45 end
46 for i =1: nx
47 for j =1: ny
48 u (1 , i ) =0;
49 u (j ,1) =0;
50 u ( ny +1 , i ) =0;
51 u (j , nx +1) =100;
52 end
53 end
54 u (1 , nx +1) =0;
55 % %% les vecteurs x et y
56 x =0: h : b ; y =0: h : d ;
57 % %% ( affichage la courbe en tenant compte des vecteurs x ,y et de la ...
matrice U )
58 figure (2)
59 mesh (x ,y , u )
6|??
22 novembre 2023
Commenter les lignes et la sortie du code.
7. On suppose que pour h = 0.1 le scripte donne une bonne approximation de la solution exacte.
Donner le code Matlab pour donner l’erreur commise entre les solutions pour h = 0.1 et h = 0.25
et tracer dans la même figure les courbes des deux solutions.
Exercice 3
∆u(x, y) + 4π 2 (x2 + y2 )u(x, y) = f(x, y), ∀(x, y) ∈ Ω = [0, 1] × [0, 1]
(P 3 ) u(0, y) = sin(πy2 ); u(1, y) = sin(π(y2 + 1))
u(x, 0) = sin(πx2 ) et u(x, 1) = sin(π(x2 + 1)), Conditions aux limites
abdellatif.elghazi@uir.ac.ma
AU = f
Expliciter U, A et f.
3. Résoudre numériquement le problème et tracer en 3D les solutions.
⊛
7|??
22 novembre 2023