Problem using ode "rkf"

Bouchra BENSIALI bouchra.bensiali at yahoo.fr
Thu Apr 29 12:06:23 CEST 2010


Hi,

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.).

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);
---
with t0=0.17, pas_temps=0.001 et tf=t0+pas_temps

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 less than that: the intermediate times used by the "rkf" method are 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.

So I defined myself a function rkf. After testing it, it returned the same results as given by integrating only x and y using the function ode.

An explanation of that mystery?
Where can I find the code of the Scilab function ode?

Thank you in advance.


      
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.scilab.org/pipermail/dev/attachments/20100429/4ee49f26/attachment.htm>


More information about the dev mailing list