[scilab-Users] Re : [scilab-Users] A question about ode "rkf"
Julio Gonzalez-Saenz
julio.gonzalez at ymail.com
Mon May 3 22:25:23 CEST 2010
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/20100503/c6c9d1c6/attachment.htm>
More information about the users
mailing list