1. Introduction

2. Equations différentielles du premier ordre

3. Equations différentielles du second ordre

4. Méthode de Runge-Kutta
    4.1. Méthode de Runge-Kutta du second ordre
    4.2. Méthode de Runge-Kutta à l'ordre 4

5. Méthode Matricielle avec des "Conditions aux Limites"

6. Conversion de coordonnées
    6.1. Coordonnées polaires
    6.2. Coordonnées cylindriques
    6.3. Coordonnées sphériques

7. Problèmes en Coordonnées Cylindriques

8. Discrétisation de l'équation de la Conduction en régime instationnaire
 

1. Introduction

Le comportement dynamique des systèmes est un sujet très important en physique. Un système mécanique par exemple, se traduit par des déplacements, des vitesses et des accélérations. Un système électrique ou électronique, se traduit par des tensions, des intensités et des dérivées temporelles sur ces quantités. Un système thermique se traduit par des températures, des gradients de température, de coefficients d’échanges thermiques, etc. En général, les équations utilisées pour décrire de tels comportements dynamiques, incluent des quantités inconnues représentant les fonctions recherchées et leurs dérivées.

Une équation qui comporte une ou plusieurs dérivées de la fonction inconnue est appelée ‘équation différentielle’, qui est représentée dans MATLAB par l’abréviation ‘ODE’. L’ordre de cette équation est déterminé par l’ordre du degré le plus élevé de la dérivation.

Les équations différentielles peuvent être classée en deux catégories : les équations différentielles avec des conditions initiales et les équations différentielles avec des conditions aux limites.

2. Equations différentielles du 1^er ordre

Les équations du 1^er ordre à conditions initiales peuvent être écrites sous la forme :

où  est une fonction de  et de , et la seconde équation représente la condition initiale sans laquelle la solution de l’équation différentielle ne peut être évaluée.

Exemple 1 :

Soit à résoudre l’équation différentielle suivante :

où :

Déterminer numériquement . On choisit un pas de temps , et on donne comme condition initiale . Tracer la solution de cette équation différentielle pour .

Solution :

L’équation  est équivalente à  où  et .

La solution de cette équation peut être effectuée par la méthode d’Euler, qui consiste à écrire :

Le programme MATLAB suivant (eul1.m) permet la résolution de cette équation :

while t<=20
    n=n+1 ;
    V=V+h*(-C/M*V.^2+g);
    t=t+h;
    t_R(n+1)=t;
    V_R(n+1)=V;
end

plot(t_R,V_R);
xlabel('Temps (en s)');
ylabel('Vitesse (en m/s)');
grid on;

On obtient ainsi un graphique (fig. 1) qui donne l’évolution de la solution  en fonction du temps.

Fig. 1 : évolution de la vitesse en fonction du temps

Exemple 2 :

On veut résoudre l’équation différentielle du 1^er ordre suivante :

en utilisant les formules d’Adams ouvertes à l’ordre 1 et les formules d’Adams (prédicteur-correcteur) pour un pas . On donne la limite de convergence .
 

• Calculer  par les deux méthodes pour 

• Ecrire un programme dans Matlab permettant de calculer et de tracer les résultats de ces deux méthodes sur le même graphe,

• Conclusion.

Solution :
* Les formules d’Adams (prédicteur-correcteur) consiste à écrire la solution sous la forme :

Dans ce cas, on a :

or  est inconnue ; donc  est aussi inconnue.

Dans ce cas, on utilise la méthode itérative suivante :

où  est la k^ième valeur itérée de , et est une valeur arbitraire initiale de . La résolution doit se faire de façon à ce que :

Si , on obtient :

On arrête les calculs si .

** Formules d’Adams ouvertes à l’ordre 1 (Méthode d’Euler) :

Pour la méthode d’Euler, la solution  est donnée par l’expression :

b) Le programme suivant ‘eul2.m’ permet d’inclure les 2 méthodes :

while t(n)<1.5
    n=n+1;
    t(n)=t(n-1)+h;
    y1(n)=y1(n-1)+h*(-y1(n-1).^1.5+1);
end

% Méthode du prédicteur-correcteur
ym(1)=10.5;t(1)=0;h=0.1;n=1;e=1e-4;n1=10000;

while t(n)<1.5
    n=n+1;
    t(n)=t(n-1)+h;
    ym(n)=ym(n-1)+h*(-ym(n-1).^1.5+1);
    % Itération succéssives et substitution
    for k=1:n1
        ymb=ym(n)+0.5*h*(-(ym(n)).^1.5-(ym(n-1)).^1.5+2);
        if abs(ym(n)-ymb)<=e;
            break;
            printf('Erreur de (%d)=%f\n',k,abs(ym(n)-ymb));
        end
        if k==n1
            fprintf('Pas de convergence pour t=%f\n',t(n));
        end
    end
 end

fprintf('t=%f\t y1=%f\t ym=%f\n',t,y1,ym);
grid on;
hold on;
plot(t,y1,t,y1,'*');
plot(t,ym,'-r');
xlabel('t');ylabel('y');
gtext('***Méthode d''Euler (ordre 1)');
gtext('-Prédicteur-correcteur');

Après exécution du programme ci-dessus, on obtient le graphique de la figure 2 ci-dessous dans lequel il y a les 2 solutions (méthode d'Euler et Méthode d'Adams).

Fig. 2 : Comparaison des deux méthodes de résolution numérique

3. Équations différentielles du second ordre

Ces équations sont du type :

où ,  et  sont des constantes ou des fonctions de .
et  sont des conditions initiales.

Avant d'appliquer la méthode d'Euler par exemple, on peut rendre l'équation précédente à une équation du type 1^er ordre. Dans ce cas, on pose .

D'où :

Les conditions initiales deviennent :

où : 

En faisant le calcul pas par pas, on obtient :

pour ,

pour ,

Ainsi, dans Matlab, le calcul pour chaque pas de temps peut être décrit par les matrices ci-dessous.

On définit au départ  et  par :

et 

ensuite, on écrit , et on résout l'équation :

Exemple :

Résoudre l'équation différentielle suivante :

             avec : 

On donne le pas .
 

•  Calculer  pour en utilisant la méthode d'Euler à l'ordre 1 (calcul manuel).

• Calculer, en utilisant un algorithme dans Matlab, la solution  de l'équation différentielle précédente pour .

• Tracer l'évolution de  pour .

Solution :
 

• a. Le système précédent peut être écrit sous la forme :

         avec : 

Pour , on a :

,

,

• b. Sous une forme vectorielle, l'équation précédente s'écrit :

     où 

La méthode d'Euler donne :

Le programme 'eul3.m', listé ci-dessous permet de calculer la solution  de l'équation différentielle proposée :

axis([0 5 -1 1]);
plot(t,y(1,:),t,y(2,:),':r');
xlabel('Temps (en s)');ylabel('v(t) et u(t)');
gtext('Solution : v(t)');
L=length(t);
text(t(L-2),y(2,L-5),'Solution : u(t)');

La fonction 'def_f(y,t)' est un sous programme dont la liste est dressée ci-dessous :

La figure 3 ci-après donne l'évolution des solutions  et . La solution  oscille autour de 0 et tend à s'amortir au fur et à mesure que le temps t augmente. On pourrait éventuellement utiliser une autre méthode de résolution numérique pour chercher .

À titre d'exemple, essayer d'utiliser la méthode itérative d'Adams (Prédicteur-Correcteur) pour retrouver la solution .
Comparer cette solution à celle d'Euler à l'ordre 1.

Fig. 3 : évolution des solutions  et  en fonction du temps

4. Méthode de Runge-Kutta

4.1. Méthode de Runge-Kutta du second ordre

Elle est de la forme :

ou bien sa forme standard peut être écrite ainsi :

Cette méthode est équivalente à celle d'Euler à deux itérations seulement.

Exemple :

Une plaque métallique épaisse à la température de 200°C (ou 473°K) est soudainement placée dans une chambre de 25°K, où la plaque est refroidie à la fois par la convection naturelle et le transfert radiatif de chaleur. On donne les constantes physiques suivantes :

En supposant que la distribution de température dans le métal est uniforme, l'équation donnant la température en fonction du temps est :

Résoudre cette équation différentielle par la méthode de Runge-Kutta à l'ordre 2 pour  et .

Solution :

La liste du programme 'RK2.m' est la suivante :

while t(i)<180
    k1=h*fonc_(T(i),Av,Epsg,hc);
    k2=h*fonc_(T(i)+k1,Av,Epsg,hc);
    T(i+1)=T(i)+0.5*(k1+k2);
    t(i+1)=t(i)+h;
    i=i+1;
end
plot(t,T);grid on;xlabel('Temps (en s)');
ylabel('Température (en °K)');

La fonction 'fonc_.m' est un sous programme dont la liste est :

La solution  obtenue de l'équation différentielle précédente est donnée par la figure 4 ci-après. Ainsi, la température obtenue, décroît en fonction du temps d'une façon hyperbolique.

4.2. Méthode de Runge-Kutta à l'ordre 4

Cette méthode s'exprime sous la forme :

et la solution  est donnée par :

Exemple 1 :

Résoudre le système suivant par la méthode de Runge-Kutta à l'ordre 4 :

avec les conditions initiales suivantes :

On donne :

Les constantes sont égales à :

Solution :

Si on définit : ;  et , le système précédent devient :

Ce système peut s'écrire sous la forme matricielle suivante :

où 

et :

Le programme suivant 'RK4.m' permet la résolution du système précédent :

while t<=30
    k1=h*f1m(y(:,i),C,F);
    k2=h*f1m(y(:,i)+k1/2,C,F);
    k3=h*f1m(y(:,i)+k2/2,C,F);
    k4=h*f1m(y(:,i)+k3,C,F);
    y(:,i+1)=y(:,i)+(1/6)*(k1+2*k2+2*k3+k4);
    t(i+1)=i*h;
    i=i+1;
end

plot(t,y(1:3,:));
grid on;
text(t(70),y(1,70),'y1');
text(t(70),y(2,70),'y2');
text(t(70),y(3,70),'y3');
xlabel('Temps (en s)');
ylabel('Solutions,y1, y2, y3');

La fonction 'f1m.m' est le sous programme listé ci-dessous :

Après exécution du programme 'RK4.m', on obtient le graphe suivant donnant les solutions  () recherchées.

Fig. 5 : évolution des trois solutions , et 

Exemple 2 :

Une baguette de longueur 0,2 m () est placée dans un écoulement d’air à la température de 293°K. La température à gauche de la baguette (en x=0) est maintenue égale à 493°K, mais le côté droit est isolé. Le transfert de chaleur dans la baguette se fait par convection. Déterminer la distribution de la température le long de l’axe de la baguette, sachant que les caractéristiques du matériau constituant cette baguette sont :

Solution :
L'équation de la conduction de la chaleur dans la direction de l'axe  de la baguette est donnée par :

avec les conditions aux limites suivantes :

La température de l'air est à 293 °K.

Ce problème est à conditions aux limites en  et . Pour le résoudre, on doit rendre les conditions aux limites en conditions initiales.

On définit :

Ainsi, l'équation différentielle (1^er ordre) ci-dessus à conditions aux limites peut être rendue en une équation différentielle du 1^er ordre à conditions initiales, en posant :

Seule, une condition initiale est connue à partir des conditions aux limites, mais la seconde condition initiale qui doit être  reste inconnue.

Cependant, on peut résoudre l'équation en faisant une transformation de la condition aux limites connue  en une condition initiale

.

Ainsi, on cherche la valeur  de telle sorte que  permet de donner .

Le programme suivant 'RK4pcl.m' permet de calculer  pour chaque valeur arbitraire de  introduite par l'intermédiaire du clavier. Une fois cette valeur  satisfait la condition , on peut arrêter les calculs en tapant la valeur –99999.

    A=0.0001;P=0.01; hc=120;k=60;b=293;l=0.2 ;
    a=P*l*hc/A/k;
    i=1;x(1)=0;h=0.01;
    y(:,1)=[493;y2];

    while x<=0.3
        k1=h*f2m(y(:,i),x(i),a,b);
        k2=h*f2m(y(:,i)+k1/2,x(i)+h/2,a,b);
        k3=h*f2m(y(:,i)+k2/2,x(i)+h/2,a,b);
        k4=h*f2m(y(:,i)+k3,x(i)+h,a,b);
        y(:,i+1)=y(:,i)+(1/6)*(k1+2*k2+2*k3+k4);
        x(i+1)=i*h;

        if (x(i)-0.2001)*(x(i)-0.1999)<0
            y2_fin=y(2,i+1);
            break;
        end
        i=i+1;
    end

    % y2_fin =(y(1,n+1)-y(1,n))/h;
    plot(x,y(1,:),'-',x,y(2,:)/10,':');
    xlabel('x (en m)');ylabel('y:- et v/10:...');
    text(0.15,-200,['y2(0.2)=',num2str(y2_fin)]);
    text(0.02,-200,['Val. arb. y2(0)=',num2str(y2)]);
    text(x(10),y(1,10)-20,'y1(x)');
    text(x(10),y(2,10)/10-20,'y2(x)/10');
    axis([0 0.2 -300 500]);
end

La fonction 'f2m.m' est un sous programme permettant l'obtention de la fonction f :

Après exécution du programme 'RK4pcl.m'; on obtient la figure 6 sur laquelle se trouvent les solutions  et .

5. Méthode Matricielle avec des "Conditions aux Limites"

Considérons l'équation de la conduction de la chaleur selon la direction  :

où la condition initiale est :

et les conditions aux limites sont :

En utilisant les différences finies, l'équation différentielle ci-dessus s'écrit :

où ,  est le nombre d'intervalles sur ,  et .

Dans ce cas, les conditions aux limites s'écrivent :

On obtient ainsi le système suivant :

où  et  sont des vecteurs colonnes et est une matrice carrée d'ordre , définis ci-après :

; et 

Ce système peut être résolu par la méthode de Runge-Kutta par exemple.

Exemple 1 :

La température d'une barre de fer de longueur  est initialement à la température de . La température du côté gauche est soudainement réduite à  à , mais la température du côté droit est maintenue constante à . Tracer la l'évolution de la température dans la barre à chaque intervalle de temps égal à , et ceci jusqu'à  (c'est à dire :, , , , , ).
Les propriétés thermiques du matériau sont :

Solution :

La diffusivité thermique est calculée par :

Pour chercher l'évolution de la température dans la barre, on divise l'intervalle  en 49 petits intervalles (par exemple). L'équation différentielle à résoudre est celle de la chaleur définie précédemment.

Les valeurs des conditions aux limites sont  et  telles que :

Le programme suivant ('RK4cc.m'), utilise la méthode de Runge-Kutta à l'ordre 4 pour la résolution de l'équation aux dérivées partielles de la conduction de la chaleur suivant la direction et en fonction du temps.

%*******************************
% Construction de la matrice A *
%*******************************
A(1,1:2)=[-2 1];A(1,3:50)=0;
A(50,1:48)=0;A(50,9:10)=[1 -2];

for i=2:49
    for j=1:50
        if i<j-1 & j>i+1
            A(i,j)=0;
        end

        if i==j
            A(i,j)=-2;
            A(i,j-1)=1;
            A(i,j+1)=1;
        end
    end
end

M=A*alpha*1/dx^2;
S(1)=TL;S(50)=TR;S(2:49)=0;S=S*alpha/dx^2;S=S';
T(1:50)=40;
T=200*ones(T);
T=T';
n=0;t=0;h=20;m=0;
axis([0 10 0 220]);
j=[0,1:length(T),length(T)+1];
T_p=[TL,T',TR];
plot(j,T_p);
text(j(2),T_p(2),['t=',int2str(t),'s']);
xlabel('i : Nombre de points');
ylabel('T (en °C)');

for k=1:5
    for m=1:50
        n=n+1;
        k1=h*(A*T+S);
        k2=h*(A*(T+k1/2)+S);
        k3=h*(A*(T+k2/2)+S);
        k4=h*(A*(T+k3)+S);
        T=T+(k1+2*k2+2*k3+k4)/6;
        t=h*n;
    end

    hold on;
    j=[0,1:length(T),length(T)+1];
    T_p=[TL,T',TR];
    plot(j,T_p);
    text(j(k+1),T_p(k+1),int2str(t));
end

Exemple 2 :

Étude de la répartition de la température dans les parois un four

On considère un four dans le plan de la section droite représentés ci-dessous. On appelle  la surface interne du four et  la surface externe.
En régime permanent, la température  en un point  de la paroi vérifie l'équation de Laplace :

avec les conditions aux limites suivantes :

Fig. 7 : Géométrie et maillage des parois du four

- Ecrire un programme calculant la température  de la paroi du four aux points de la grille définie sur la figure 7 ci-dessus (figure en haut).
On suppose que les dimensions , ,  et  permettent ce quadrillage.

Application Numérique :

     et : 
 

Solution :

Équation de Laplace :

avec :

Le développement en série de Taylor de la fonction  autour de  où  s'écrit :

En additionnant ces deux relations et en divisant par , on trouve :

En opérant de la même manière pour la variable , on trouve :

Dans notre cas, on a . Donc, le Laplacien bidimensionnel s'écrit :

On pose :

Pour simplifier la notation, l'expression précédente s'écrit sous la forme :

              (*)

Cette équation peut être représentée par la forme moléculaire suivante :

Pour  et  par exemple, on obtient 26 nœuds à l'intérieur du domaine (fig. 7). Donc, , où  est l'ensemble des points intérieurs du domaine (et ;  sur  et  sur ) .

Écrivons les premières et dernières lignes du système (*) :

Remarquons que :

et que :

On a le système suivant alors :

On a donc un système de 26 équations à 26 inconnues. Donc, ce système (ou l'équation de Laplace discrétisée) peut être résolu(e) par la méthode de Gauss-Seidel (méthode explicite) où l'équation de Laplace donne le point central de chaque forme moléculaire par l'expression :

                     (**)

Remarque :
Avec l'algorithme précédent (**), il est inutile de mémoriser explicitement la matrice du système (*). Il suffit de mémoriser la température  de chacun des points de la grille.

Définissons les valeurs particulières des indices correspondants aux parois de la figure 8 suivante :

Le programme suivant appelé 'laplace.m' permet de résoudre le système (**) par la méthode explicite. Pour cela, on se donne une valeur (distribution) arbitraire initiale , qui portée dans l'équation (**) au second membre pour chaque couple , donne une nouvelle valeur , et ainsi de suite. L'arrêt des calculs se fait quand  où  est la limite de convergence que l'on se donne.

% Initialisation de la température dans le four
for i=1:n
    for j=1:m
        T(i,j)=Thetaint;
    end
end

% Température de la paroi externe
for i=1:n
    T(i,1)=Thetaext;
    T(i,m)=Thetaext;
end

for j=1:m
    T(1,j)=Thetaext;
    T(n,j)=Thetaext;
end

% Température de la paroi interne
for i=n1:n2
    T(i,m1)=Thetaint;
    T(i,m2)=Thetaint;
end

for j=m1:m2
    T(n1,j)=Thetaint;
    T(n2,j)=Thetaint;
end

% Méthode de Gauss-Seidel (Itérations)
for k=1:k1
    for i=2:n-1
        for j=2:m1-1
            T(i,j)=0.25*(T(i-1,j)+T(i+1,j)+T(i,j-1)+T(i,j+1));
        end
    end

    for i=2:n-1
        for j=m2+1:m-1
            T(i,j)=0.25*(T(i-1,j)+T(i+1,j)+T(i,j-1)+T(i,j+1));
        end
    end

    for i=2:n1-1
        for j=m1:m2
            T(i,j)=0.25*(T(i-1,j)+T(i+1,j)+T(i,j-1)+T(i,j+1));
        end
    end

    for i=n2+1:n-1
        for j=m1:m2
            T(i,j)=0.25*(T(i-1,j)+T(i+1,j)+T(i,j-1)+T(i,j+1));
        end
    end

    if abs(T(n-1,m-1)-T(2,2))<=eps
        fprintf('\n \n');
        fprintf('Températures après "%d itérations\n',k);
        fprintf('\n \n');
        break;
    end
end

for i=1:n
    fprintf('%5.0f\t',T(i,1:m));
    fprintf('\n');
end

% Tracé des courbes de température en 3D
hold off;
if n==m
    figure(1);
    i=1:n;
    j=1:m;
    grid on;
    [x,y]=meshgrid((i-1)*dx,(j-1)*dx);
    mesh(x,y,T);
    Title('Evolution de la température dans le four');
    xlabel('i');ylabel('j');zlabel('T (en °C)');
end

% Tracé des courbes isothermes en 2D
figure(2);
i=1:n;j=1:m;
grid on;
contour(i,j,T(i,j),15);
title('Lignes isothermes dans les parois du four');
xlabel('i');ylabel('j');
t_mis=toc
Nb_opt=flops

On exécute ce programme (laplace.m), en rentrant les valeurs de , , , ,  et enfin  et .

>>Introduire la valeur de L1 :
60
Introduire la valeur de L2 :
60
Introduire la valeur de L3 :
20
Introduire la valeur de L4 :
20
Introduire la valeur du pas dx :
1
Introduire la valeur de Theta interne :
1220
Introduire la valeur de Theta externe :
35
t_mis=
113.9700
Nb_opt=
7362728

6. Conversion de coordonnées

Dans Matlab, il existe plusieurs commandes pour la conversion des coordonnées. Nous en donnons ci-après un bref apperçu.
 

6.1. Coordonnées polaires

Les fonctions 'cart2pol' et 'pol2cart' permettent respectivement le passage des coordonnées cartésiennes en coordonnées polaires et inversement. Les syntaxes sont :

>>[theta,r]=cart2pol(x,y)

>>[x,y]=pol2cart(theta,r)

et  doivent être des vecteurs de même taille représentant les coordonnées des points considérés,

r et  représentent les coordonnées polaires et sont des vecteurs de mêmes dimensions que  et , où  est exprimée en radians.
 

6.2. Coordonnées cylindriques

Les fonctions 'cart2pol' et 'pol2cart' permettent respectivement le passage des coordonnées cartésiennes en coordonnées cylindriques et inversement. Les syntaxes sont :

>>[theta,r,z]=cart2pol(x,y,z)
>> [x,y,z]=pol2cart(theta,r,z)

,  et  sont des vecteurs de même taille,

et  sont des vecteurs de mêmes dimensions que ,  et , et  est exprimée en radians.

6.3. Coordonnées sphériques

Les fonctions 'cart2sph' et 'sph2cart' permettent respectivement le passage des coordonnées cartésiennes en coordonnées sphériques et inversement. Les syntaxes sont :

>>[Az,Elv,r]=cart2sph(x,y,z)

>>[x,y,z]=sph2cart(Az,Elv,r)

r, ,  et  sont des vecteurs de même taille,  et  sont respectivement l'azimut et l'élévation exprimées en radians, de mêmes dimensions que ,  et .

7. Problèmes en Coordonnées Cylindriques

Considérons l'équation de la conduction de la chaleur exprimée en coordonnées cartésiennes :

En coordonnées polaires, cette équation s'écrit sous la forme :

sur la domaine

Sur ce domaine , on construit un maillage en coordonnées polaires  comme le montre la figure 11 ci-dessous, avec  et  où et  sont des entiers.

Ainsi, la température  et la fonction  deviennent au point  :

Pour des valeurs non nulles de , les dérivées secondes de  par rapport à  et  s'écrivent sous la forme disrétisée suivante :

En utilisant les différences centrées,  peut s'écrire au point  sous la forme :

Ainsi, l'équation de la chaleur discrétisée est :

Soit, après regroupement des différents termes :

C'est l'équation de conduction de la chaleur discrétisée (différences finies) en coordonnées cylindriques, obtenue pour des valeurs de  non nulles (). Les indices i et  sont des entiers (commençant par 1 dans Matlab).

A l'origine (), l'équation de la conduction en coordonnées polaires présente une singularité. Cette dernière doit être éliminée. Pour cela, on utilise le Laplacien de l'équation de la conduction non pas en coordonnées cylindriques, mais en coordonnées cartésiennes :

         quand 

On construit ensuite un cercle de rayon , centré en . Considérons  la température à , et , ,  et  sont les températures sur le cercle aux 4 nœuds (intersection avec les axes  et oy). Ainsi, l'équation précédente (en coordonnées cartésiennes) s'écrit sur le cercle de rayon  :

La rotation des axes  et oy autour de  conduit aux mêmes résultats (équation de la chaleur en coordonnées cartésiennes discrétisée). En considérant  comme la moyenne arithmétique des températures autour du cercle de rayon , l'équation précédente devient :

à 

où  est la moyenne arithmétique des valeurs de  autour du cercle de rayon  et de centre  et  est la valeur de la température à .

Les coordonnées polaires  -deux dimensions- peuvent être extrapolées en coordonnées cylindriques (trois dimensions) pour l'obtention de l'équation de la conduction de la chaleur discrétisée.

Pour le problème bidimensionnel en  évoqué précédemment, compte tenu de la symétrie axiale, l'équation de la conduction de la chaleur peut être réduite à :

Pour , l'équation précédente discrétisée s'écrit sous la forme :

où , , et  est un entier positif.

Au centre (), en utilisant la règle de l'Hopital nous obtenons :

Ainsi l'équation de la conduction de la chaleur en deux dimensions  s'écrit compte tenu de la symétrie axiale du problème :

en 

Soit en différences finies :

pour 

En coordonnées cylindriques, l'équation de la conduction de la chaleur est donnée par l'expression suivante :

Les coordonnées  sont représentées par :

où et  sont des entiers.

La température  au nœud  est notée par :

et les différentes dérivées partielles deviennent :

L'équation de la chaleur discrétisée en coordonnées cylindriques au nœud  devient :

pour des valeurs de  non nulles.

Pour , on a :

et l'équation de conduction en  devient :

en 

Soit en utilisant les différences finies au nœud  :

8. Discrétisation de l'équation de la Conduction en régime instationnaire

Dans le système de coordonnées cartésiennes (trois dimensions), l'équation de conduction de la chaleur en régime instionnaire (temporel) s'écrit (si ) :

Pour résoudre cette équation aux dérivées partielles en utilisant la méthode des différences finies, on utilise un maillage cubique (mailles de côtés  avec :

où ,  et  sont entiers positifs, et le domaine temporel est divisé en petits intervalles de temps  de telle sorte que : 

Dans ce cas, la température  au point P(x,y,z) à un temps donné  est représentée par :

En utilisant la méthode des différences finies, les différents termes de l'équation ci-dessus au dérivées partielles s'écrivent au point P(x,y,z) :

et l'équation de la conduction discrétisée en trois dimensions devient :

En regroupant les différents termes de cette équation, on obtient :

Cette équation ne peut être résolue que par une méthode itérative (méthode de Gauss-Seidel par exemple).

Si on pose :

alors, la condition de convergence de la solution recherchée dépend essentiellement du signe de la quantité. Cette quantité doit être strictement positive.

Ainsi, une condition doit relier le pas de temps  et les autres pas spatiaux (, et ) :

Pour amorcer le calcul itératif, on ensemence le domaine de frontière  par des valeurs arbitraires ( par exemple), et on arrête les calculs quand la condition suivante sera réalisée :

où  est la précision que l’on se fixera.