résolution d'un système d'équation différentielle, schéma implicite

**autocadNUL** · 19/04/2011, 10h12

Bonjour
je doit résoudre un système d'équation différentielle de la forme :

dT/dt=f[T,r,q(t)]
dr/dt=g(T)

J'ai d'abord pensé à pde ou ode, mais ces deux outils son plutôt mal adaptées, la première servant a résoudre des équations conique et la seconde une équadiff d'une seule variable. J'ai donc opté pour une résolution direct, en mettant à rude épreuve mes connaissance en simulation numérique.

Le schéma de discrétisation est implicite et j'utilise la méthode du point fixe à chaque pas de temps pour déterminer T (boucle de type x=f(x)). Voici le coeur de mon code (attention T noté Tp):

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
%boucle temporelle
for i=2:N
 
    temp0=Tp(i-1,j);
    temp1=0;
    temp2=Tp(i-1,j);
    rayon0=r(i-1,j);
    test=0;
    rayon=r(i-1,j);
    %boucle de résolution de l'équation différentielle
    while (abs(temp2-temp1)>prec)&(test==0)
        %%%%calcul des caractéristique non constantes%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
        if rayon<0
            rayon=0;
            temp2=Tg;
            test=1;
        else
            %flux de masse [kg/m².s]
            mu=((M_v/(2*pi*R_GP*temp2))^0.5)*Pref*exp(delta_H*(temp2-Tref)/(R_GP*temp2*Tref));
            rayon=rayon0-mu/rho_s;
 
            %capacacité calorifique des suies [J/kg.K]
            Cp_p=(16.86+0.00477*temp2-85400/(temp2*temp2))/0.012;
 
            %emmisivité
            eps=48*pi*rayon*n_r*n_i/(lambda*(4*n_r*n_r*n_i*n_i+(2+n_r*n_r-n_i*n_i)*(2+n_r*n_r-n_i*n_i)));
 
            %capacité calorifique gaz
            somme2=0;
            dTemp=(temp2-Tg)/Ni;
            for ka=1:Ni
                tampon=28.58+0.00377*temp2-(5e4)/(temp2*temp2);
                tampon=tampon-R_GP;
                somme2=somme2+tampon*dTemp;
            end
            integrale=somme2/M_v;
 
            %coefficient de transfert radiatif
            h_r=sigma*(Te*Te+temp2*temp2)*(Te+temp2);
 
            %%%%calcul du temps d'integration suivant%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
            temp1=temp2;
            tampon=dt/(3*rayon*rho_s*Cp_p);
            temp2=(temp0+tampon*(eps*(q(i,j)/4+h_r*Te)-mu*(delta_H/M_s+alpha*(-R_GP*Tg/(2*M_v)+integrale))))/(1+tampon*(eps*h_r+alpha*mu*R_GP/(2*M_v)));
        end
    end
    Tp(i,j)=temp2;
    r(i,j)=rayon;

Voilà le problème c'est que rien ne marche, tout se passe comme si je n'avais aucune relaxation... Pardonnez moi mon code brouillon et mal adapté à matlab, mais je suis plus abitué au fortran, et je ne vois pas quoi faire d'autre. Peut être auriez vous un outil spécifique ou bien dénicherez vous l'erreur qui viens surement de la méthode numérique utilisé...

Merci d'avance

**autocadNUL** · 19/04/2011, 16h15

Rebonjour
Il me semble que j'ai trouvé un moyen d'utiliser l'outil ode pour résoudre mon système. C'est pas gagné, mais le schmilblick avance

. Par contre je ne suis pas convaincu par ode45 qui est un mélange de schéma implicite et explicite (Runge Kutta), et ayant déja essayé un schéma explicite, je peux vous dire que j'ai de l'instabilité

. Connaîtriez vous un ode qui marche sur un schéma purement implicite?
encore merci

**FR119492** · 20/04/2011, 00h54

Salut!
Dans la série des sous-programmes ODExx, il y en a d'autres que ODE45 qui sont spécialement destinés à l'intégration des équations différentielles raides (angl: stiff).
Jean-Marc Blanc

**autocadNUL** · 20/04/2011, 09h00

bonjour
oui tu as raison, apparement ode23s correspond à un schéma implicite.
Par contre a tu jeté un coup d'oeil a mon code? sa me turlupine de ne pas voire où sa bug...

**autocadNUL** · 20/04/2011, 16h04

rebonjour
voilà, j'ai un nouveau problème avec ode23s... Je n'ai pas assez de point. Il en fraudais 1000, et il ne m'en sort que 23 au maximum, du coup mes résultat sont complètement faux! y a il un moyen de "forcer" ode a accepter une échelle de temps?
merci d'avance

**autocadNUL** · 20/04/2011, 16h37

bon, j'ai trouvé tout seul, désolé d'avoir abusé du flood...
Je passe juste pour livrer ma solution :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
ti=[0:dt:tmax];
for j=1:M
    q_r=qr(j);
    [t,Y]=ode23s('T_fonction',ti,[T0 r0]);
    T(:,j)=Y(:,1);
    r(:,j)=Y(:,2);
end

et la fonction :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
function dY=T_fonction(t,Y)
 
...
global q_r
...
 
T=Y(1);
r=Y(2);
 
...
 
dT=tampon*((eps*q/4+sigma*(Te^4-T^4))-mu*(delta_H/M_s+alpha*((R_GP*(T-Tg)/(2*M_v)+somme2))));
dr=-mu/rho_s;
dY=[dT;dr];
 
return

pour passer mes argument j'utilise des variables globales, c'est dangereux mais j'ai pas vu d'autre solution...
pour "forcer" le pas de temps, j'utilise directement ti pour l'argument tspan.
pour résoudre un système au lieu d'une seule équation, Y est le système lui même. Enfin l'usage de ode23s permet d'éviter les divergence due à un schéma explicite instable. Bon courage