This is the mail archive of the
gsl-discuss@sources.redhat.com
mailing list for the GSL project.
Re: ODEs - runge kutta
- From: Slaven Peles <peles at cns dot physics dot gatech dot edu>
- To: "Sam Halliday" <plendily at hotmail dot com>,gsl-discuss at sources dot redhat dot com
- Date: Thu, 21 Feb 2002 11:52:06 -0500
- Subject: Re: ODEs - runge kutta
- References: <F177K8r6beUzmRf6Ci600001d70@hotmail.com>
On Thursday 21 February 2002 11:19 am, Sam Halliday wrote:
> hello,
>
> i am currently trying to numerically solve the lane-emden equation, which
> is a 2nd order differential equation, and experience has shown that the 4th
> order runge kutta is the best way to solve this.
>
> i was interrested in using the GSL for solving this problem, although the
> examples section for ODE's scared me too much... must i really define so
> much, (jacobians and so forth), is there not a way in which i can just send
> my function, parameter (which actually is the time) and time step length to
> a library and let it do all the work? does anyone have a simple example i
> can follow so that i can use this library for this task? teh one in the
> documenattion seems to be too crypic, and tries to do too much for my
> needs...
>
It looks scary, indeed. But, at this point it seems to me that results
obtained with gsl are much more accurate than those obtained using Numerical
Recipes, for example. So I guess it is worth to put a little effort in
learning how to use gsl. Below is a working example for solving Van der Pol
equation using RK4 fixed stepsize, without defining Jacobian. This is the
best what I could do at the moment, so if anybody has some better example,
please send it to the list.
Cheers,
Slaven
#include <stdio.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_odeiv.h>
int
func (double t, const double y[], double f[],
void *params)
{
double mu = *(double *)params;
f[0] = y[1];
f[1] = -y[0] - mu*y[1]*(y[0]*y[0] - 1);
return GSL_SUCCESS;
}
int
jac (double t, const double y[], double *dfdy,
double dfdt[], void *params)
{
double mu = *(double *)params;
gsl_matrix_view dfdy_mat
= gsl_matrix_view_array (dfdy=NULL, 2, 2);
dfdt[0] = 0.0;
dfdt[1] = 0.0;
return GSL_SUCCESS;
}
int
main (void)
{
const gsl_odeiv_step_type * T
= gsl_odeiv_step_rk4;
gsl_odeiv_step * s
= gsl_odeiv_step_alloc (T, 2);
double mu = 10;
gsl_odeiv_system sys = {func, jac, 2, &mu};
double t = 0.0, t1 = 100.0;
double h = 1e-2;
double y[2] = { 1.0, 0.0 }, y_err[2], dfdy[4], dydt_in[2], dydt_out[2];
/* initialise dydt_in */
GSL_ODEIV_JA_EVAL(&sys, t, y, dfdy, dydt_in);
while (t < t1)
{
int status = gsl_odeiv_step_apply (s, t, h, y, y_err,
dydt_in, dydt_out, &sys);
if (status != GSL_SUCCESS)
break;
dydt_in[0] = dydt_out[0];
dydt_in[1] = dydt_out[1];
t += h;
printf("%.5e %.5e %.5e\n", t, y[0], y[1]);
}
gsl_odeiv_step_free(s);
return 0;
}