[scilab-Users] Re : [scilab-Users] A question about ode "rkf"

[Bouchra Bensiali]

Thank you Julio for your confirmation,

I checked that when tending the "pas_temps" to zero (t0, tf unchanged), the results tend to be the same.

What about the first example?

If a Scilab developer is near here, could (s)he give an explanation please? 

Thanks.




________________________________
De : Julio Gonzalez-Saenz <julio.gonzalez at ymail.com>
À : users at lists.scilab.org
Envoyé le : Lun 3 mai 2010, 22 h 25 min 23 s
Objet : Re: [scilab-Users] Re : [scilab-Users] A question about ode "rkf"


Hi Bouchra,

I am not a SCILAB developer, but normally an ODE uses a variable step size based on the relative error, in the case (x,y) it uses the you have a x-component of the error and a y-component of the error.. so is in the first case, just x-component.. Therefore if you start where the error is big, you have small step size, while if you have small error, you can increase the step size..

I would like to hear an answer from the development team about this, but 90% sure your approach is correct...

Regards,

Julio




________________________________
From: Bouchra Bensiali <bouchra.bensiali at yahoo.fr>
To: users at lists.scilab.org
Sent: Mon, 3 May, 2010 19:50:30
Subject: [scilab-Users] Re : [scilab-Users] A question about ode "rkf"


Hi,

I come with concrete examples, and wish I can find some help.

After some research www.math.univ-metz.fr/~sallet/ODE_Scilab.pdf, it seems that the methods used by ode choose their steps to control errors. They are not fixed step solvers.

But have some questions:

1/ Consider this:
=========================================================================
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
=============================================================================
If you try this example, you'll find that the intermediate times used between t=0.178 and t=0.179 (to predict the solution at t=0.179) are not the same in the first and second integration. Why the time t=0.178 is not handeled the same way? A depandance to what happened before maybe?


2/ My initial problem was that integrating xdot = f(t,x), x(t0)=x0, or integrationg together xdot=f(t,x), ydot=g(t,y), x(t0)=x0, y(t0)=y0, don't give the same result for x. An example that illustrates this is the following:
==============================================================================
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,:)
===================================================================================
From what I understood, the "rkf" method, as the ode45 (Matlab) is a variable step solver: if the estimated error is bigger than the tolerance, the stepsize is reduced. In the second integration, the error is related to x and y, while in the first integration, it's only related to x, so that can explain the difference in the results. Confirmation?

Thanks.




________________________________
De : Bouchra Bensiali <bouchra.bensiali at yahoo.fr>
À : users at lists.scilab.org
Envoyé le : Ven 30 avril 2010, 16 h 21 min 41 s
Objet : [scilab-Users] A question about ode "rkf"


Hi,

I wonder if the "rkf" method uses a fixed number of intermediate times as said here http://math.fullerton.edu/mathews/n2003/RungeKuttaFehlbergMod.html.

For example, if I define a function F:
--
function Xdot=F(t,X)

  //printf("t=%f \n",t)
  
  Xdot = ... ;
  
endfunction;
--
and make a first integration:
--
pas_temps = 0.001;
t0 = 170*pas_temps;
tf = 180*pas_temps;

X0 = ... ;

Trajectoire1 = ode("rkf",X0,t0,t0:pas_temps:tf,F);
--
and a second integration beginning with the result obtained at time t=0.178 (for example):
--
t0 = 178*pas_temps;
tf = 180*pas_temps;

X0 = Trajectoire1(:,9);

Trajectoire2 = ode("rkf",X0,t0,t0:pas_temps:tf,F);
--
then is Trajectoire2 supposed to be the same as Trajectoire1(:,9:11)?

Because I have an example that gives different results and found (by using printf("t=%f \n",t)) that the number of intermediate times changed from one integration to another.

Thank you in advance for your help. I'm really stuck.


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