Examains

Examains Master
Master Informatique
Informatique (I):
(I):
Méthodes
Méthodes Numériques
Numériques
Abdoulkadri Chama

Enoncé
Enoncé de
de la
la Première
Première Session
Session
Exercice no 1 :
1. Soit f une fonction continue et dérivable sur R. Donner le développement de Taylor d’ordre un de f (x + h) et
f (x − h) en x.
2. En déduire l’approximation de la dérivée première de f par la méthode des différences finies centrées et avant.
3. En supposant que f soit une fonction à deux variables x et y, et en utilisant les différences centrées premières
∂f (x, y) ∂f (x, y) ∂ 2 f (x, y)
de et , déterminer l’approximation de la dérivée seconde de f , c’est-à-dire .
∂x ∂y ∂x2 2

Exercice no 2 :
La température u(x, y, t) dans une région bidimensionnelle évolue selon l’équation de la chaleur suivante :
 2
∂ u ∂2u

∂u
=α + 2 (1)
∂t ∂x2 ∂y
où :
— u(x, y, t) est la température en position (x, y) et au temps t.
— α = 0.1 est la constante de diffusivité thermique.
Les conditions aux limites spécifient la température aux bords de la région Ω, et la condition initiale définit la
distribution initiale de température au temps t = 0.
Les conditions aux limites spécifiées sont :
u(0, y, t) = u(a, y, t) = 0
u(x, 0, t) = u(x, b, t) = 0 .
La condition initiale est donnée par :
u(x, y, 0) = f (x, y) (2)
où f (x, y) = 50 sin(x · y) est la distribution initiale de température.
On considère les notations u(xi , yj , tn ) = unij où xi = ih, yj = jh et tn = n∆t. Avec ∆t = 0.01, h = 0.01, xi = ih,
yj = jh et tn = n∆t. i = 0, . . . , N .
1. Écrire l’approximation de l’équation de la chaleur en utilisant le résultat de l’exercice 1 et les notations précé-
dentes, ou plus précisément unij . Cette approximation est aussi appelée forme discrète de l’équation.
2. Calculer les premières valeurs de u, c’est-à-dire u1ij .

1
 Exercice no 3 :
Le code en C++ ci-joint implémente l’équation de la chaleur (1) décrite dans l’exercice 2. Veuillez compléter les
parties du code correspondant a la fonction "initialize" (pour initialiser la matrice de température) et la fonction
"apply_boundary_conditions" (pour appliquer les conditions aux limites).

1 # include < iostream >

2 # include < fstream >
3 # include < vector >
4 # include < cmath >
5 # include < iomanip >
6

8 const double alpha = 0.1;

9 const double dx = 0.1;
10 const double dy = 0.1;
11 const double dt = 0.01;
12 const int nx = 21;
13 const int ny = 31;
14 const double a = 2.0;
15 const double b = 3.0;
16 const int nt = 100;
17

18

19 void initialize ( std :: vector < std :: vector < double > >& u ) {
20 }
21

22

23 void a p p l y _ b o u n d a r y _c o n d i t i o n s ( std :: vector < std :: vector < double > >& u ) {
24 }
25

26

27 void update ( std :: vector < std :: vector < double > >& u , std :: vector < std :: vector < double > >& un ) {
28 for ( int i = 1; i < nx -1; ++ i ) {
29 for ( int j = 1; j < ny -1; ++ j ) {
30 double uxx = ( un [ i +1][ j ] - 2* un [ i ][ j ] + un [i -1][ j ]) / ( dx * dx ) ;
31 double uyy = ( un [ i ][ j +1] - 2* un [ i ][ j ] + un [ i ][ j -1]) / ( dy * dy ) ;
32 u [ i ][ j ] = un [ i ][ j ] + alpha * dt * ( uxx + uyy ) ;
33 }
34 }
35 }
36

37

38 int main () {
39 std :: vector < std :: vector < double > > u ( nx , std :: vector < double >( ny , 0) ) ;
40 std :: vector < std :: vector < double > > un ( nx , std :: vector < double >( ny , 0) ) ;
41

42 initialize ( u ) ;
43

44 std :: ofstream outfile ( " heat_data . csv " ) ;

45 outfile << std :: setprecision (6) ;
46

47

48 for ( int i = 0; i < nx ; ++ i ) {

49 double x = i * dx ;
50 for ( int j = 0; j < ny ; ++ j ) {
51 double y = j * dy ;
52 outfile << x << " , " << y << " , 0 , " << u [ i ][ j ] << " \ n " ;
53 }
54 }
55 outfile << " \ n " ;
56

57 for ( int t = 1; t < nt ; ++ t ) {

58 double time = t * dt ;
59 un = u ;
60 a p p l y _ b o u n d a r y _ co n d i t i o n s ( u ) ;
61 update (u , un ) ;
62

63

2
64 for ( int i = 0; i < nx ; ++ i ) {
65 double x = i * dx ;
66 for ( int j = 0; j < ny ; ++ j ) {
67 double y = j * dy ;
68 outfile << x << " , " << y << " , " << time << " , " << u [ i ][ j ] << " \ n " ;
69 }
70 }
71 outfile << " \ n " ;
72 }
73

74 outfile . close () ;
75 return 0;
76 }

Correction
Correction du
du Sujet
Sujet de
de la
la Première
Première Session
Session
Exercice no 1 :
1) On suppose que f est une fonction continue et de classe C n sur R, avec n ≥ 2.
Le développement de Taylor de f (x + h) autour de x d’ordre 1 est donné par :

f (x + h) = f (x) + hf ′ (x) + O(h2 )

De même, le développement de Taylor de f (x − h) autour de x d’ordre 1 est donné par :

f (x − h) = f (x) − hf ′ (x) + O(h2 )

Approximation de la dérivée première par la méthode des différences finies en avant :
La méthode des différences finies en avant pour approximer la dérivée première de f est donnée par :

f (x + h) − f (x)
f ′ (x) ≈
h
En utilisant le développement de Taylor de f (x + h), on a :

f (x + h) − f (x)
= f ′ (x) + O(h)
h
Ainsi, l’approximation de la dérivée première de f à l’aide des différences finies en avant est :

f (x + h) − f (x)
f ′ (x) ≈ avec une erreur de l’ordre de O(h).
h
2) Le développement de Taylor d’ordre 2 de f (x + h) autour de x est donné par :

h2 ′′
f (x + h) = f (x) + hf ′ (x) +
f (x) + O(h3 )
2
De même, le développement de f (x − h) d’ordre 2 autour de x est donné par :

h2 ′′
f (x − h) = f (x) − hf ′ (x) + f (x) + O(h3 )
2
En additionnant ces deux développements, on obtient :

f (x + h) + f (x − h) = 2f (x) + h2 f ′′ (x) + O(h3 )

Approximation de la dérivée seconde par la méthode des différences finies centrées :
L’approximation de la dérivée seconde de f à l’aide des différences finies centrées est donnée par :

f (x + h) − 2f (x) + f (x − h)
f ′′ (x) ≈
h2
En utilisant le développement de Taylor, on trouve que cette approximation a une erreur de l’ordre de O(h2 ).
Ainsi, la méthode des différences finies centrées pour approximer la dérivée seconde de f est :

f (x + h) − 2f (x) + f (x − h)
f ′′ (x) ≈
h2
avec une erreur d’approximation de l’ordre de O(h2 ).

3
 Exercice no 2 :
La température u(x, y, t) dans une région bidimensionnelle évolue selon l’équation de la chaleur suivante :
 2
∂ u ∂2u

∂u
=α +
∂t ∂x2 ∂y 2
où u(x, y, t) est la température en position (x, y) au temps t, et α = 0.1 est la constante de diffusivité thermique.
Les conditions aux limites sont :
u(0, y, t) = u(a, y, t) = 0, u(x, 0, t) = u(x, b, t) = 0
La condition initiale est donnée par :
u(x, y, 0) = f (x, y) = 50 sin(x · y)
On utilise les notations u(xi , yj , tn ) = unij , où xi = ih, yj = jh, tn = n∆t, avec ∆t = 0.01, h = 0.01, et i = 0, . . . , N .
1. L’approximation de l’équation de la chaleur en utilisant la méthode des différences finies est donnée par :

un+1 − unij uni+1,j − 2unij + uni−1,j uni,j+1 − 2unij + uni,j−1

 
ij
=α +
∆t h2 h2
Ce qui se réorganise en :
uni+1,j − 2unij + uni−1,j uni,j+1 − 2unij + uni,j−1
 
un+1
ij = unij + α∆t +
h2 h2
2. Les valeurs initiales de u au temps t = 0 sont données par :
u0ij = 50 sin(xi · yj )
Pour calculer u1ij , on utilise l’approximation de l’équation discrète. Par exemple, pour i = 1, j = 1, on calcule :

u021 − 2u011 + u001 u012 − 2u011 + u010

 
=u111 u011
+ α∆t +
h2 h2
Ainsi, en répétant ce processus, on peut calculer les valeurs de u pour t = ∆t.

1 # include < iostream >

2 # include < fstream >
3 # include < vector >
4 # include < cmath >
5 # include < iomanip >
6

7 // D f i n i t i o n des p a r a m t r e s du p r o b l m e
8 const double alpha = 0.1; // Coefficient de diffusion
9 const double dx = 0.01; // Pas spatial en x
10 const double dy = 0.01; // Pas spatial en y
11 const double dt = 0.01; // Pas temporel
12 const int nx = 201; // Nombre de points en x (0 2 avec un pas de 0.01)
13 const int ny = 301; // Nombre de points en y (0 3 avec un pas de 0.01)
14 const double a = 2.0; // Limite s u p r i e u r e en x
15 const double b = 3.0; // Limite s u p r i e u r e en y
16 const int nt = 100; // Nombre de pas de temps
17

18

19 void initialize ( std :: vector < std :: vector < double > >& u ) {
20 for ( int i = 0; i < nx ; ++ i ) {
21 for ( int j = 0; j < ny ; ++ j ) {
22 double x = i * dx ;
23 double y = j * dy ;
24 u [ i ][ j ] = 10* sin ( x * y ) ;
25 }
26 }
27 }
28

29

30 void a p p l y _ b o u n d a r y _c o n d i t i o n s ( std :: vector < std :: vector < double > >& u ) {
31 for ( int j = 0; j < ny ; ++ j ) {
32 u [0][ j ] = 0; // u (0 , y , t ) = 0
33 u [ nx -1][ j ] = 0; // u (a , y , t ) = 0
34 }
35 for ( int i = 0; i < nx ; ++ i ) {
36 u [ i ][0] = 0; // u (x , 0 , t ) = 0
37 u [ i ][ ny -1] = 0; // u (x , b , t ) = 0

4
38 }
39 }

Solution
Solution de
de la
la Deuxième
Deuxième Session
Session
Exercice no 1 :
1. Ecriture en terme de courant : (2pts)
Le courant dans le circuit est donné par :

dQ
I(t, Q) =
dt
En substituant cela dans l’équation différentielle :

1 V (t)
I(t, Q) + Q=
RC R
En simplifiant, on obtient :

V (t) 1
I(t, Q) = − Q
R RC
2. Méthode de Runge-Kutta d’ordre 2 : (2pts)

Soit l’équation différentielle de la forme y ′ = f (t, y). L’algorithme de la méthode de Runge-Kutta d’ordre 2 est
donné par :

k1 = hf (tn , yn )

k2 = hf (tn + h, yn + k1 )

1
yn+1 = yn + (k1 + k2 )
2
où h est le pas de temps, tn et yn sont les valeurs de temps et de la fonction au pas n.

3. Calcule des charges : (6pts)

La charge Q(0.1) : (3pts)
— Valeurs initiales :
t0 = 0, Q(0) = 0

5
 — Calcul de k1 :    
V Q0 10 0
k1 = h · − = 0.1 · − = 0.1 · 5 = 0.5
R R·C 2 2 · 0.5
— Calcul de k2 :    
V Q0 + k1 10 0 + 0.5
k2 = h · − = 0.1 · −
R R·C 2 2 · 0.5
k2 = 0.1 · (5 − 0.5) = 0.1 · 4.5 = 0.45
— Mise à jour de Q(0.1) :
1
Q(0.1) = Q0 + (k1 + k2) = 0 + 0.475 = 0.475 Coulombs
2

La charge Q(0.2) : (3pts)

On utilise Q(0.1) = 0.475 C comme condition initiale pour le pas suivant :

— Valeurs initiales :
t1 = 0.1 s, Q(0.1) = 0.475 C
— Calcul de k1 :    
V Q(0.1) 10 0.475
k1 = h · − = 0.1 · −
R
 R·C  2 2 · 0.5
0.475
k1 = 0.1 · 5 − = 0.1 · 4.525 = 0.4525
1
— Calcul de k2 :
 
V Q(0.1) + k1
k2 = h · −
R R·C
 
10 0.475 + 0.4525
= 0.1 · −
2 2 · 0.5
= 0.40725
— Mise à jour de Q(0.2 s) :
1
Q(0.2) = Q(0.1) + (k1 + k2 ) = 0.475 + 0.429875 = 0.904875 Coulombs
2

Exercice no 1 :
Le code en C++ ci-après implémente la méthode de Runge-Kutta d’ordre 2 pour résoudre l’équation différentielle
associée à la charge Q(t) dans le circuit RC, tel que décrit dans l’exercice précédent. Veuillez compléter les parties
du code correspondantes à la fonction I(·, ·)(pourcalculerladérivéedelacharge)etlaf onction

1 # include < iostream >

2 # include < fstream >
3 # include < cmath >
4

6 const double R = 2.0;

7 const double C = 0.5;
8 const double V = 10.0;
9 const double h = 0.1;
10 const int num_steps = 1000;
11

12

13

14 double I ( double t , double Q ) {

15 return ( V / R ) - (1.0 / ( R * C ) ) * Q ;
16 }
17

18

19

20

21 double runge_kutta ( double t , double Q , double h ) {

22 double k1 = I (t , Q ) ;
23 double k2 = I ( t + h , Q + h * k1 ) ;
24 return Q + ( h / 2.0) * ( k1 + k2 ) ;
25 }

6
26

27

28 int main () {
29 double t = 0.0;
30 double Q = 0.0;
31

32

33 std :: ofstream outfile ( " charge_Q . csv " ) ;

34 if (! outfile . is_open () ) {
35 std :: cerr << " Erreur d ’ ouverture du fichier " << std :: endl ;
36 return 1;
37 }
38

39 outfile << " Temps , Charge \ n " ;

40

41 for ( int i = 0; i < num_steps ; ++ i ) {

42

43 outfile << t << " ," << Q << " \ n " ;

44

45 Q = runge_kutta (t , Q , h ) ;
46 t += h ;
47 }
48

49

50 outfile . close () ;
51

52 std :: cout << " Les d o n n e s sont s a u v e g a r d e s dans charge_Q . csv pour la visualisation . "
<< std :: endl ;
53

54 return 0;
55 }

Visualisation
Le code en Python ci-dessous permet de lire les données générées par le programme en C++ et de tracer la courbe
du potentiel Q en fonction du temps t.
1 import matplotlib . pyplot as plt
2 import pandas as pd
3

4 data = pd . read_csv ( ’ charge_Q . csv ’)

5

7 plt . figure ( figsize =(10 , 6) )

8 plt . plot ( data [ ’ Temps ’] , data [ ’ Charge ’] , label = ’ Charge Q ( t ) ’ , color = ’ blue ’)
9 plt . xlabel ( ’ Time ( t ) ’)
10 plt . ylabel ( ’ Charge ( Q ) ’)
11 plt . title ( ’ Charge Q vs Time for an RC Circuit using Runge - Kutta 2 nd Order ’)
12 plt . legend ()
13 plt . grid ( True )
14 plt . savefig ( ’ charge_vs_time . png ’)
15 plt . show ()

Exécution des deux fichiers

Pour compiler les deux codes précédents, utilisez le script bash ci-dessous. Ce dernier permet d’automatiser la visua-
lisation des résultats. Une autorisation administrative du fichier est nécessaire avant l’exécution.
1 # !/ bin / bash
2

3 # Compile the C ++ code

4 g ++ -o RK2_RC RK2_RC_ENONCE . cpp
5

6 # Run the C ++ program to generate the CSV file

7 ./ RK2_RC
8

9 # Check if the data file was generated

10 if [ -f charge_Q . csv ]; then

7
11 echo " CSV data file generated successfully . "
12 else
13 echo " Error : CSV data file was not generated . "
14 exit 1
15 fi
16

17 # Run the Python script to visualize the data

18 python3 plot_charge . py
19

20 echo " Plotting complete . Check the generated ’ charge_vs_time . png ’ for the plot . "