[Top][All Lists]
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: [Help-gsl] ODE solver: oscillation around expected value
From: |
Jeroen Meidam |
Subject: |
Re: [Help-gsl] ODE solver: oscillation around expected value |
Date: |
Mon, 01 Aug 2011 13:50:35 +0100 |
User-agent: |
Mozilla/5.0 (Windows; U; Windows NT 6.1; en-GB; rv:1.9.2.18) Gecko/20110616 Lightning/1.0b2 Thunderbird/3.1.11 |
In python I used the scipy.odeint solver.
I can give an example of what I am doing in the case of an adaptive stepper.
In the main I create the vector with initial conditions, calculated
before then I calculate an integration time, I give it a sampling rate
and from that I calculate the number of calculations that the solver needs.
double y[4] = { r0, phi0, pr0, pphi0 } ;
calculate_all(y, var, con) ; /*calculate all variables at t=0.*/
con->tfs = 25.0 * pow( ( con->M/(2.8*con->MO) ), -5.0/3.0) *
pow(con->f0_GW/40.0), -8.0/3.0) ; /*Integration time in seconds*/
con->tf = con->tfs/con->M ; /*Integration time*/
con->SR = 8192.0 ; /*sample rate [s^-1]*/
con->N = (int)(con->SR * con->tfs) ; /*number of steps*/
con->delta = con->tf/con->N ; /*stepsize [t/M]*/
Next I call this function which contains the solver:
void ode_solver_adaptivestep(double* y, struct variables* var, struct
constants* con, FILE *fp) {
void *params = var ;
gsl_odeiv2_driver *d ;
double err = con->accuracy ;
double hstart = con->delta ; /* The initial stepsize */
gsl_odeiv2_system sys = {func, jac, 4, params};
switch (con->simulation_type) { /* To easily access different stepper
functions*/
case 5:
printf ("Starting simulation: Adaptive timestep, rk2...\n") ;
d = gsl_odeiv2_driver_alloc_y_new (&sys, gsl_odeiv2_step_rk2, hstart,
err, err);
break ;
case 6:
printf ("Starting simulation: Adaptive timestep, rk4...\n") ;
d = gsl_odeiv2_driver_alloc_y_new (&sys, gsl_odeiv2_step_rk4, hstart,
/* And there are more options, but I want to save space */
}
int i,ti;
double t = 0.0 ;
double amplitude_h, omg ;
int evaluations ;
evaluations = con->N;
for (i = 1; i <= evaluations; i++) {
ti = i * con->tf / evaluations;
calculate_all(y, var, con) ; /* This function calculates all
variables needed for each step*/
amplitude_h = pow(var->omega,1.0/3.0)*cos(2.0*y[1]) ;
omg = var->omega ;
int status = gsl_odeiv2_driver_apply (d, &t, ti, y);
if (status != GSL_SUCCESS)
{
printf ("error, return value=%d\n", status);
break;
}
if (y[0] < con->rLR) {
printf ("Simulation has reached light ring at time %.3f and
iteration %d.\n", t, i) ;
break ;
}
fprintf(fp, "%.5f;%.5f;%.5f;%.5f;%.5f;%.5f;%.5f\n", t, amplitude_h,
y[0], y[1], y[2], y[3], omg);
}
gsl_odeiv2_driver_free (d);
}
Does that give any insight into what might be a problem?
On 1-8-2011 13:32, Benjamin wrote:
On 08/01/2011 02:08 PM, Jeroen Meidam wrote:
Hi,
I'm using the gsl ode solver to solve a system of 4 coupled
differential equations. I've done the exact same thing in python, but
needed to rewrite it in c. There is one component of the plotted
solutions that shows an unexpected oscillatory behaviour around the
value that python gives as aresult.
I have no idea what could cause this. I have tried all the available
stepping functions and none appear to get rid of this oscillation.
Also, I would expect more differences between the different solving
methods. If I plot for example the rk2 , rk4, rkf45 and rk8pd
solutions in one graph, I can't see the slightest bit of difference.
Is that to be expected?
_______________________________________________
Help-gsl mailing list
address@hidden
https://lists.gnu.org/mailman/listinfo/help-gsl
Hi,
I am not an expert on this, but if you use the adaptive stepping
algorithm from odeiv the different solving methods could give very
similar results because of the adaptive control of errors.
I would suspect that if you would plot their respective step sizes you
should see a difference.
As for the oscillation I have no idea but I think it would in general
a good idea if you'd also post your source code/a working example so
that it will be easier to suggest something.
Maybe you can also mention how you implemented the solution in python,
did you write your own solver routine or did you for example use the
scipy.odeint solver?