Problem using ode "rkf"

Bouchra BENSIALI bouchra.bensiali at yahoo.fr
Wed Apr 28 15:42:33 CEST 2010


Hello world !

I have a problem using ode "rkf". I defined a function like that :
---
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
---
then I used : 
---
Sol = ode("rkf",X0,t0,t0:pas_temps:tf,F);
---

In this case, xdot and ydot don't depend on A=X(3), I obtain however a different result for x=X(1) and y=X(2) compared with the result obtained by integrating only x and y (same initial condition, etc.).

I noticed that it depended on how I define Adot. For example, doing this :
---
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);
---
with t0=0.17, pas_temps=0.001 et tf=t0+pas_temps

leads to :
---
t=0.170000
t=0.170250
t=0.170375
t=0.170923
t=0.171000
t=0.170500
t=0.171000
---

I verified that these times correspond to the intermediate times used by the "rkf" method : t0, t0+1/4*pas_temps, t0+3/8*pas_temps, t0+12/13*pas_temps, t0+pas_temps, t0+1/2*pas_temps, as said here http://math.fullerton.edu/mathews/n2003/RungeKuttaFehlbergMod.html.

Now if I use for example :
---
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);
---
 
then Scilab returns that :
---
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
---
and as a consequence that changes the results for x and y (they are supposed not to change because they don't depend on A). I checked the definition of v3 and don't see errors. For me, the number of intermediate times is supposed to be the same, am I right ?

An explanation of that mystery ?

Thank you in advance.


      
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.scilab.org/pipermail/users/attachments/20100428/3d6623b1/attachment.htm>


More information about the users mailing list