Discussion :

[bouwaha]

Bonjour,

J'ai un problème en utilisant ode "rkf" dans Scilab.
J'ai une fonction

Code :

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

1
2
3
4
5
6
7
8
function Xdot = F(t,X)
  [v1,v2,v3] = fonc(t,X(1),X(2));
  xdot = v1;
  ydot = v2;
  Adot = v3;
 
  Xdot = [xdot; ydot; Adot];
endfunction

Et j'utilise :

Code :

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

1
2
 
Sol = ode("rkf",X0,t0,t0:pas_temps:tf,F);

Dans mon cas xdot et ydot ne dépendent pas de A=X(3), j'obtiens cependant un résultat différent de ce que me donne l'intégration uniquement de x=X(1) et y=X(2) en utilisant les mêmes conditions initiales etc.

En fait, j'ai fait un test sur un seul pas de temps, et si je prends par exemple Adot = v1; au lieu de Adot = v3; je n'ai pas ce problème. J'ai remarqué en incluant un printf("%f",X(1)) dans la focntion F, que le nombre de points intermédiaires change en fonction de ce que je définis pour Adot, est-ce normal ? je dirais que non vu que la méthode "rkf" Runge-Kutta-Fehlberg est une méthode à 6 pas.

Merci d'avance de votre aide et de vos éclaircissements.
_
Bou

[bouwaha]

Bonjour,

J'avoue que mon problème n'est pas très clair, mais comme la fonction fonc fait appel à d'autres fonctions etc., je ne peux pas tout inclure.

Mais par exemple en faisant ceci :

Code :

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

1
2
3
4
5
6
7
8
9
10
11
12
function Xdot = F(t,X)
  printf("t=%f \n",t);
 
  [v1,v2,v3] = fonc(t,X(1),X(2));
  xdot = v1;
  ydot = v2;
  Adot = 1;
 
  Xdot = [xdot; ydot; Adot];
endfunction
 
Sol = ode("rkf",X0,t0,t0:pas_temps:tf,F);

avec t0=0.17, pas_temps=0.001 et tf=t0+pas_temps;

j'obtiens :

Code :

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

1
2
3
4
5
6
7
t=0.170000
t=0.170250
t=0.170375
t=0.170923
t=0.171000
t=0.170500
t=0.171000

J'ai vérifié que ça correspondait bien aux temps intermédiaires utilisés par "rkf" à savoir : t0, t0+1/4*pas_temps, t0+3/8*pas_temps, t0+12/13*pas_temps, t0+pas_temps, t0+1/2*pas_temps, comme expliqué ici.

En revanche faire cela :

Code :

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

1
2
3
4
5
6
7
8
9
10
11
12
function Xdot = F(t,X)
  printf("t=%f \n",t);
 
  [v1,v2,v3] = fonc(t,X(1),X(2));
  xdot = v1;
  ydot = v2;
  Adot = v3;
 
  Xdot = [xdot; ydot; Adot];
endfunction
 
Sol = ode("rkf",X0,t0,t0:pas_temps:tf,F);

donne :

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
t=0.170000
t=0.170000 
t=0.170250 
t=0.170375 
t=0.170923 
t=0.171000 
t=0.170500 
t=0.170218 
t=0.170326 
t=0.170803 
t=0.170870 
t=0.170435 
t=0.170870 
t=0.170903 
t=0.170919 
t=0.170990 
t=0.171000 
t=0.170935 
t=0.171000

et du coup, je n'obtiens pas du tout le même résultat pour x, y qui, je le rappelle, ne dépendent pas de A.

Je ne comprends pas pourquoi ça change le nombre de points intermédiaires. Qu'est-ce qui pourrait être à l'origine de cela ? sacahant que j'ai vérifié la définition de v3 dans la fonction fonc, et je ne vois pas d'erreur.

D'avance merci.

[bouwaha]

Bonjour,

Désolée d'insister, mais ça me casse la tête depuis quelques jours ...

Je vais simplifier le problème.

Supposons que j'aie cela :

Code :

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

1
2
3
4
5
6
7
8
 
function Xdot=F(t,X)
 
  //printf("t=%f \n",t)
 
  Xdot = ... ;
 
endfunction;

Je fais une première intégration :

Code :

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

1
2
3
4
5
6
7
8
 
pas_temps = 0.001;
t0 = 170*pas_temps;
tf = 180*pas_temps;
 
X0 = ... ;
 
Trajectoire1 = ode("rkf",X0,t0,t0:pas_temps:tf,F);

Et une deuxième en partant du résultat obtenu par la première intégration à l'instant t = 0.178 :

Code :

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

1
2
3
4
5
6
7
 
t0 = 178*pas_temps;
tf = 180*pas_temps;
 
X0 = Trajectoire1(:,9);
 
Trajectoire2 = ode("rkf",X0,t0,t0:pas_temps:tf,F);

Est-ce que je suis censée avoir Trajectoire1(:,9:11)==Trajectoire2 ? Parce qu'il m'arrive d'avoir des résultats différents ...

J'espère que c'est plus clair.
Merci.

[bouwaha]

Bonjour,

Je reviens avec un exemple concret. J'espère que vous pourrez m'aider.

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
 
function Xdot=F(t,X)
  printf("t=%f \n ",t)
 
  Xdot = ((1-t)*X-X^2)/(10^(-4));
 
endfunction
 
//
pas_temps=0.001;
t0=170*pas_temps;
tf=180*pas_temps;
 
X0=1;
 
Trajectoire1=ode("rkf",X0,t0,t0:pas_temps:tf,F);
 
//
printf("\n \n");
 
t0=178*pas_temps;
tf=180*pas_temps;
 
X0=Trajectoire1(:,9);
 
Trajectoire2=ode("rkf",X0,t0,t0:pas_temps:tf,F);
 
//comparaison
Dif = Trajectoire1(:,9:$)-Trajectoire2

On trouve que les temps intermédiaires utilisés entre l'instant t=0.178 et t=0.179 (pour prédire la solution à t=0.179) ne sont pas les mêmes entre la première et la seconde intégration. Pourquoi ces instants ne sont pas traités de la même façon ? Une dépendance par rapport à ce qui s'est passé avant ?

__
Mon problème initial est que intégrer :

Code :

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

1
2
3
 
xdot = f(t,x)
x(t0) = x0

ou intégrer :

Code :

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

1
2
3
4
5
 
xdot = f(t,x)
ydot = g(t,y)
x(t0) = x0
y(t0) = y0

ne donnent pas le même résultat pour x.

Un exemple qui illustre cela :

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
 
function Xdot=F1(t,X)
  printf("t=%f \n ",t)
 
  Xdot = X/10^(-3);
 
endfunction
 
//
pas_temps=0.001;
t0=170*pas_temps;
tf=172*pas_temps;
 
X0=1;
 
Trajectoire1=ode("rkf",X0,t0,t0:pas_temps:tf,F1);
 
//
printf("\n \n");
 
function Xdot=F2(t,X)
  printf("t=%f \n ",t)
 
  Xdot = [X(1)/10^(-3);((1-t)*X(2)-X(2)^2)/(10^(-4))];
 
endfunction
 
//
pas_temps=0.001;
t0=170*pas_temps;
tf=172*pas_temps;
 
X0=[1;1];
 
Trajectoire2=ode("rkf",X0,t0,t0:pas_temps:tf,F2);
 
//comparaison
Dif = Trajectoire1 - Trajectoire2(1,:)

D'après ce que j'ai compris, le code essaie de réduire le pas en temps si l'erreur estimée dépasse l'erreur tolérée. Dans la seconde intégration, l'erreur prend en compte x et y, alors que dans la première, il n'y a que x qui intervient. Ceci explique peut-être cette différence dans les résultats. (?)

Mais je trouve quand-même bizarre qu'une intégration de façon vectorielle d'équations complètement découplées donne un résulat différent de l'intégration de chacune des équations à part.

Merci d'avance de votre aide.

PS : ode45 est un peu l'équivalent Matlab de ode "rkf" dans Scilab.