Problème de solver ode

Version imprimable

j'ai un systeme equation diff du type suivant:
Code:

1 2 3 dx(s_i,t)/dt=f(x(s_i,t), x(s_i,t)\sum_i x(s_i,t), y(s_i,t), y(s_i,t)\sum_i y(s_i,t), z(s_i,t), z(s_i,t)\sum_i z(s_i,t)) dy(s_i,t)/dt=g(x(s_i,t), x(s_i,t)\sum_i x(s_i,t), y(s_i,t), y(s_i,t)\sum_i y(s_i,t), z(s_i,t), z(s_i,t)\sum_i z(s_i,t)) dz(s_i,t)/dt=h(x(s_i,t), x(s_i,t)\sum_i x(s_i,t), y(s_i,t), y(s_i,t)\sum_i y(s_i,t), z(s_i,t), z(s_i,t)\sum_i z(s_i,t))
avec x, y, z: une fonction de s_i et de t et i=1:50 (s_i c'est une discretisation d'une deuxiemme variables).
f,g,h: des fonctions de type polynome
je peux vous envoyez les fonctions si cela est important

Je sais demontre mathematiquent que si la condition initiale dans [0,1] alors la solution reste dans [0,1].

Malheureusement qt je simule les solutions et pour certaine valeurs des parametres (intervalle asses petit), l'une des composante de la solution qui devrai tendre vers 0 commence par etre negative du type -0.8e-15 et ensuite cette composante devient clairement negative (~-1000).

Je me pose la question si c'est un probleme du solver ou un que j'ai une erreur dans l'ecriture de ma fonction

merci pour votre aide

16/10/2008, 14h15
FR119492

Salut!
Une explication possible serait que ton système différentiel est raide (stiff). Essaie soit de prendre un pas d'intégration plus petit, soit d'utiliser un sous-programme spécialement prévu pour ce cas (ode23s).
Jean-Marc Blanc
16/10/2008, 14h37
salseropom

Salut,
pour les problèmes raides je conseille un solveur qui s'appelle CVode. Il est possible d'appeler ce solveur dans Matlab.

A plusieurs reprises, les solveurs ode23s, ode15s, ode23tb etc... m'ont déçu car ils ne me donnaient jamais la bonne solution ... alors que CVode oui.

Citation:

Envoyé par slimane66

Je sais demontre mathematiquent que si la condition initiale dans [0,1] alors la solution reste dans [0,1].

Tu as une ODE dans R^3 et tu parles de conditions initiales dans [0, 1] donc dans R. N'y aurait-il pas une erreur de dimension dans tes conditions initiales.
16/10/2008, 14h46
slimane66

merci a tout deux
1. Effectivement, j'ai voulu dire [0,1]^{3}
2. qt tu dit raide (stiff) cela veux dire quoi exactement. La notion de raide fait reference a la courbe des solutions ?
16/10/2008, 16h05
FR119492

Salut!
Dit sommairement, un problème est dit raide si sa solution peut comporter des composantes à des échelles de temps très différentes.

Exemple: Tu veux étudier le comportement d'un moteur de locomotive démarrant en direction du Gothard. A l'enclenchement du moteur, tu as un régime transitoire avec des constantes de temps de l'ordre de quelques millisecondes; tu as aussi les équation du démarrage du train, dont la constante de temps est de l'ordre de la minute; tu as enfin l'échauffement progressif des bobinages de ton moteur et l'augmentation correspondante de leur résistance, dont les constantes de temps sont de l'ordre de l'heure. Si tu essaies de mettre tout ça ensemble, tu as un système raide.

Jean-Marc Blanc
16/10/2008, 20h19
slimane66

j'ai essaye le solver edo23s et le pb des nombre negatif dispare !!!!
j'ai installe sundials et je suis en train d'apprendre a utiliser cvode, j'ai encore un peut de mal en plus je n'ai pas trouver de help pour matlab !!! est ce que que vous en connaisez un
merci encore a vous deux

Citation:

Envoyé par FR119492

Salut!
Dit sommairement, un problème est dit raide si sa solution peut comporter des composantes à des échelles de temps très différentes.

Exemple: Tu veux étudier le comportement d'un moteur de locomotive démarrant en direction du Gothard. A l'enclenchement du moteur, tu as un régime transitoire avec des constantes de temps de l'ordre de quelques millisecondes; tu as aussi les équation du démarrage du train, dont la constante de temps est de l'ordre de la minute; tu as enfin l'échauffement progressif des bobinages de ton moteur et l'augmentation correspondante de leur résistance, dont les constantes de temps sont de l'ordre de l'heure. Si tu essaies de mettre tout ça ensemble, tu as un système raide.

Jean-Marc Blanc
16/10/2008, 21h56
salseropom

1 pièce(s) jointe(s)

Salut,
sundials fournit normalement des routines pour faire des mex files et ainsi tu peux appeler CVode directement à partir de matlab (sinon tu réécrits ton code en C... mais c'est un autre problème)

Pour continuer sur ce qu'à dit Jean-Marc : je t'envoie en pièce jointe une figure obtenue par un calcul. Auncun des solveurs de matlab ne trouvait la solution mais CVode a réussi (comme solution exacte j'ai pris la solution calculée par un logiciel industriel).

Les problèmes raides apparaissent tout le temps dans les problèmes de cinétique chimique

bonjour,
effectivement mes solutions ressemblent a celle que tu montres (pratiquement constante -- un saut --- constante).
Ce sont des equations issues de modeles biologiques qui a mon avis de sont pas tres eloigne de celles de la chimie.

J'aurai besoin d'une petite aide pour sundials sous matlab:
j'ai installe sundials et j'ai fait tourner les exemples qui existent dans l'archive, tout se passe tres bien.

J'ai commence par le script suivant:

Code:

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
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
%CVDX - CVODES example problem (serial, dense)
%   The following is a simple example problem, with the coding
%   needed for its solution by CVODES. The problem is from
%   chemical kinetics, and consists of the following three rate
%   equations:         
%      dy1/dt = -.04*y1 + 1.e4*y2*y3
%      dy2/dt = .04*y1 - 1.e4*y2*y3 - 3.e7*(y2)^2
%      dy3/dt = 3.e7*(y2)^2
%   on the interval from t = 0.0 to t = 4.e10, with initial
%   conditions: y1 = 1.0, y2 = y3 = 0. The problem is stiff.
%   While integrating the system, we also use the rootfinding
%   feature to find the points at which y1 = 1e-4 or at which
%   y3 = 0.01. This program solves the problem with the BDF method,
%   Newton iteration with the CVDENSE dense linear solver, and a
%   user-supplied Jacobian routine.
%   It uses a scalar relative tolerance and a vector absolute
%   tolerance. Output is printed in decades from t = .4 to t = 4.e10.
%   Run statistics (optional outputs) are printed at the end.
%
%   See also: cvdx_f, cvdx_g, cvdx_J
 
% Radu Serban <radu@llnl.gov>
% Copyright (c) 2005, The Regents of the University of California.
% $Revision: 1.3 $Date: 2006/03/15 19:31:26 $
clear all;
close all;
 
data.p = [0.04; 2.0e4; 3.0e7];
 
t0 = 0.0;
y0 = [1.0;0.0;0.0];
 
options = CVodeSetOptions('RelTol',1.e-8,...
                          'AbsTol',[1.e-8; 1.e-14; 1.e-6],...
                          'LinearSolver','Dense',...
                      'JacobianFn',@cvdx_J,...    
                      'RootsFn',@cvdx_g, 'NumRoots',2);
%                          
 
mondata.sol = true;
mondata.mode = 'text';
mondata.skip = 10;
mondata.update = 100;
options = CVodeSetOptions(options,'MonitorFn',@CVodeMonitor,'MonitorData',mondata);
 
CVodeMalloc(@cvdx_f,t0,y0,options,data);
 
t1 = 0.4;
tmult = 10.0;
nout = 10;
 
% fprintf('---------------------------------------------------\n');
% fprintf('     t        q       h\n');
% fprintf(['          y   '...
%          '          yd1'...
%          '          yd2\n']);
% fprintf('---------------------------------------------------\n');
 
iout = 0;
tout = t1;
while iout < nout
 
  [status,t,y] = CVode(tout,'Normal');
 
% Compute solution derivative in two ways:
 % [yd1, new_data] = cvdx_f(t, y, data);
%  yd2 = CVodeGet('DerivSolution', t, 1);
 
% Extract statistics
%  si = CVodeGetStats;
 
% Print output
 % fprintf('%0.4e    %1d  %0.4e',t, si.qlast, si.hlast);
 % if(status == 2)
%     fprintf(' ... Root found  %d   %d\n',si.RootInfo.roots(1), si.RootInfo.roots(2));
%   else
%     fprintf('\n');
%   end
%   for i = 1:3
%     fprintf('  %14.6e  %14.6e  %14.6e\n', y(i), yd1(i), yd2(i));
%   end
%   fprintf('---------------------------------------------------\n');
 
% Update output time
  if(status == 0)
    iout = iout+1;
    tout = tout*tmult;
  end
 
end
 
si = CVodeGetStats;
 
CVodeFree;

J'ai commente certaines partie pour comprendre comment le programme fonctionne.
Ma question est la suivante:
1. cvode ne donne pas les valeurs de y pour tout les pas de temps ? J'ai cru comprendre qu'il ne donne que les valeurs pour le temps final c'est-a-dire: tout
2. a quoi sert la partie du programme :

Code:

1
2
3
4
5
% Update output time
  if(status == 0)
    iout = iout+1;
    tout = tout*tmult;
  end

En fait je croit que je n'ai pas bien compris l'algorithme du programme.
merci d'avance

17/10/2008, 09h16
salseropom
Salut,

pour ma part, j'avais réécrit mon code en C. Ceci dit, dans l'exemple de CVode, effectivement, seul le dernier temps est conservé. Change cela de la manière suviante
Code:

1 2 3 4 5 6 while(t < tfin) 1) tu calcules ton EDO 2) tu écris sur le disque dur (ou bien tu stocke en mémoire) la solution de ton EDO (à l'instant t) 3) tu incrémentes le temps : t = t + dt end
puis tu post-traite tes résultats car tu les as à tous les temps

Dans ton message tu as écrits :
Code:

1 2 3 4 5 6 % Update output time if(status == 0) iout = iout+1; tout = tout*tmult; end
==> lorsque j'ai débuté avec CVode, je me suis posé la même question. CVode résout l'EDO entre t et tout (tout = t + dt). A la fin t vaut tout donc il faudra ensuite résoudre jusqu'à tout + dt

Mais dans l'exemple que tu as donné, tout = tout * tmult => suite géométrique, mais on s'en moque. L'idée à retenir est que tout est l'instant final d'intégration. Tu continues à intégréer jusqu'à ce que tout >= tfin.

Pour CVode, comme ton problème est raide, utilise la méthode BDF, Newton (ce sont des mots-clés à mettre, lit la doc). Enfin, met un ordre maximal de 2 (car à plus que 2 il peut y avoir des instabilité). Regarde les options de CVode.

Si tu veux, je peux te passer mon programme (mais il est écrit en C). Tu pourras t'en inspirer.
17/10/2008, 10h10
slimane66

rebonjour,
oui avec plaisir je souhaiterai avoir ton sript en C.

Arrete moi si je me trompe:
on est d'accord que la fonction CVode donne la valeur de y au temp tout. Pour calculer la valeur de y elle utilise differentes methodes d'integration, le choix du pas de temps est laisse au solveur (a mon avis on pourrai le changer en utilisant CVodeSetOption).
Si ce que je dit est juste le choix de "dt" ne doit pas choisi par rapport a la precision du pas d'integration mais juste comme points sortie du programme soit pour visualise de la courbe, soit pour utilise la valeur de y a d'autre fin !

is it like that ???
17/10/2008, 21h45
salseropom

Salut,

oui, CVode calcule la solution en t et tout (tout = t + dt). Donc a tout tu peux écrire la valeur de ta solution (ou la garder en memoire). Bref, à tous les dt, tu connais ta solution.

CVode calcule la solution entre t et t + dt. Entre ces 2 instants, CVode choisit lui-même son pas de temps. Toi, tu n'as une maitrise que sur la valeur de dt.

Ne prend pas un dt trop grand. Il n'y a pas de règle absolu. tout dépend de ton système d'EDO.
17/10/2008, 21h51
salseropom

2 pièce(s) jointe(s)

Voilà mon script (enfin une partie car j'ai effectué ma thèse dans le cadre d'un contrat industriel donc je ne peux pas tout divulguer). Les fonctions essentielles sont dans mon script.

Je joins aussi une figure que j'ai obtenu lors de mes calculs. On voit très bien que Matlab est dans les choux et que j'obtiens avec CVode la "bonne" solution (la référence étant celle en rouge, calculée avec Chemkin (Chemikal kinetic), un logiciel industriel).

Là on s'éloigne du forum "Matlab". Mais mon expérience (voir la figure en PJ) montre qu'il faut se méfier des routines toutes faites (idem qu'il faut se méfier de CVode je pense). Mais pour rester sur le forum Matlab, j'ai bien dit qu'il existe une possibilité pour brancher CVode à Matlab.

Je sais que CVode est un solveur très utilisé dans le monde de la recherche. Il a une limite (il parait) : lorsque la dynamique des variables est très différentes (donc on pondère les variables rapides par un epsilon très petit...)

N'hésite surtout pas à me poser des questions sur mon script (écrit en C) et sache que le forum C est là pour te répondre (je vais pas mal sur ce forum)
17/10/2008, 23h02
FR119492

Salut!
Regarde aussi www.unige.ch/~hairer/software.html. Je pense que Gerhard Wanner et Ernst Hairer sont parmi les meilleurs spécialistes de l'intégration des systèmes différentiels.
Jean-Marc Blanc
19/10/2008, 21h46
salseropom

Citation:

Envoyé par FR119492

Salut!
Regarde aussi www.unige.ch/~hairer/software.html. Je pense que Gerhard Wanner et Ernst Hairer sont parmi les meilleurs spécialistes de l'intégration des systèmes différentiels.
Jean-Marc Blanc

Salut,
super ton lien. Je le garde précieusement dans un coin.
Merci.

merci beaucoup a tous les deux, ça y est j'arrive depuis vendredi a faire tourner cvode:
par contre j'ai encore deux problème:

1. si je fais :

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
 
while tout < tmax
 
  [status,t,y] = CVode(tout,'Normal');
 
%Update output time
 if(status == 0)
   %iout = iout+1;
   tout = tout+dt;%
  %tout =tout*tmult;
  end
 
end

ou bien

Code:

1
2
3
4
 
tout=tmax
 
  [status,t,y] = CVode(tout,'Normal');

j'obtiens 2 choses différentes, c'est bizarre ????
Avez vous une idee du pourquoi ?

2. ci joint mon script

Code:

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
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
% Ce programme permete de lire les valeurs des parametres dans un fichier:
% Cb_ln_b<4.xls ou Ca_ln.xls etc etc. Pour chaque jeux de parametre il
% effectue la simulation de ODE dy/dt=fonctionSum(x,t,s,integral phi(s)ds,
% parametres) avec x\in\R3
% On discretise la variable s, s_i. Ce qui nous donne alors x\in\R3^nbre de s_i
 
%% Definition des parametres
clear all;
close all;
 
global nbs;
global ds s1 smax smin;
 
 
% parametre lies a l'integration
%Ca: cout des gamette, Cb: cout du chgt de sexe, Ceps: cout des combat
%entre male
paramEtude=2; % 1:Ca, 2:Cb et 3:Ceps
 
% type du parametre a etudie
axe=['a','b','c'];
xaxis=axe(paramEtude);
 
%% ouveture du fichier qui contiens les parametres
paramSet=xlsread('Cb_ln_b<4.xls');
%paramSet=xlsread(['C' xaxis '_ln.xls']); %
 
% Fichier contenant les parametres fixe pour notre edo: smin, smax etc etc
paramProg=xlsread('paramProg.xls');
 
nbreExp=size(paramSet,1); %Nbre de simultaion
 
%% Definition des parametre fixe
 
%Discretisation de s
smin=paramProg(1);
 smax=paramProg(2); %le chifre 11 correspond au plus petit telque max g(s)~1 avec g: ln
 ds=paramProg(3);
 s1=sort(smin+(smax-smin)*rand(101,1));%(smin:ds:smax)'; %smin+(smax-smin)*rand(1,101)
 nbs=size(s1,1);
 
% parametre du temps d'integration 
 tmin=paramProg(4);
% tmax=paramProg(5);
 %dt=paramProg(6);
 %T=[0 tmax];
 t0=tmin;
 
%% Definition de l'etat initiale
 
[Fnc0,Fc0,Mc0,Mnc0]=condinit(nbs);
y0=[Fnc0';Fc0';Mc0'];
 
%% Calcul des resultats
%gamma(s1,s1)
%fonction(1,etatinitiale)
 
%%%%%%
 
 
options = CVodeSetOptions('RelTol',1.e-5,...
                          'AbsTol',1.e-5,...
                          'LinearSolver','Dense',...  
                           'LMM','BDF',...
                      'NonlinearSolver','Newton',...
                      'MaxOrder',2);
%                 'JacobianFn',@cvdx_J,...
%                  'RootsFn',@cvdx_g, ...
%                  'NumRoots',2,...
 
%mondata.sol = true;
%mondata.mode = 'text';
%mondata.skip = 10;
%mondata.update = 100;
%options = CVodeSetOptions(options,'MonitorFn',@CVodeMonitor,'MonitorData',mondata);
t1=0.4;
tmult = 10.0;
nout = 5;
 
for k=1:nbreExp
    %Definition des parametre a etudie de mon edo
   param.Ca=paramSet(k,1);
   param.Cb=paramSet(k,2);
   param.Ceps=paramSet(k,3);
 
   CVodeMalloc(@fonctionSum,t0,y0,options,param);
 
   tout = t1;
   iout = 0;
 
while iout < nout
  [status,t,y] = CVode(tout,'Normal');
 
 % Extract statistics
 % si = CVodeGetStats;
 
 %Update output time
 if(status == 0)
   iout = iout+1;
   %tout = tout+dt;%
  tout =tout*tmult;
  end
 
end
 
si = CVodeGetStats;
CVodeFree;
 
% POur chaque jeux de parametre je stoke les solution
   dim_t=size(res,1);
   Fnc(k,1:nbs)=transpose(y(1:nbs));
   Fc(k,1:nbs)=transpose(y(nbs+1:2*nbs));
   Mc(k,1:nbs)=transpose(y(2*nbs+1:3*nbs));
end
 
%% Partie traitement des resultats
Je ne l'ai pas mis

mon problème est le suivant:
au bout de 5 boucles de k j'ai le message d'erreur suivant et le programme se met en boucle:

Citation:

Internal t=127.703 and h= 1.04795e-15 are such that t+h=t on the next step. The solver continue anyway
at t=127.703, nxstep steps taken before reaching tout
at t=127.703 and h=2.7 ... e-197, the corrector convergence test failed repeatedly

le h

Avez vous une idée pourquoi ?

amicalement

25/10/2008, 20h38
salseropom
Salut, pour répondre à tes questions

1) c'est normal. si tu fais
Code:

1 2 3 4 5 6 7 8 9 10 11 while tout < tmax [status,t,y] = CVode(tout,'Normal'); %Update output time if(status == 0) tout = tout+dt; end end
CVode va chercher un pas (adaptatif) entre t et tout. Il va ensuite "approximer" la solution pour que tu aies la valeur à exactement t = tout.

Si tu fais
Code:

1 2 3 4 tout=tmax [status,t,y] = CVode(tout,'Normal');
il n'y a plus cette approximation. C'est donc la meilleure chose à faire. Dans mon programme, j'étais obligé de faire la 1e solution (avec le tout = tout + dt) car je couplais différents modèles. Mais si tu veux résoudre une "brave" équation différentielle sur [0, tmax] fais la 2e méthode. Mais là tu n'auras que la valeur finale. La 1e méthode te permets de stocker les résultats intermédiaires et donc de tracer des courbes

2) c'est normal. quand tu résouts entre [t et t + dt] CVode utilise un pas adaptatif (qu'il appelle h). Ici ton h vaut 1e-15 donc il stagne. Quand j'avais ce problème, c'est que j'avais une erreur de programmation. Sinon, essaie d'augmenter ton pas de temps dt. Pour ma part, je te déconseille de faire un tout = tout * tmult mais fait plutôt un tout = tout + dt. Tu pourras ainsi stocker tous les résultats à un même incrément.

Autre chose, lorsque j'avais moi aussi un problème raide, pour minimiser les itérations je fixais tmin et tmax, avec tmin < tmax tel que la raideur qui se produit à l'instant t_stiff soit telle que tmin < t_stiff < tmax (bien sûr, comme tu ne connais pas t_stiff tu bornes large).
Ensuite si 0 < t < tmin ou tmax < t < tfin tu poses dt = dtmax
si tmin < t < tmax tu poses dt = dtmin.
Ainsi, tu auras une solution fine au niveau de la raideur, et "grossière" ailleur. Mais fait des essais sur la valeur de tmin, tmax, dtmin et dtmax car il ne faut pas non plus trop s'éloigner de la solution réelle (que tu ne connais pas mais que tu peux calculer avec dt ultra petit)