This is the mail archive of the
gsl-discuss@sources.redhat.com
mailing list for the GSL project.
Questions on bsimp
- From: "Di Xiao" <dxiao at physics dot utexas dot edu>
- To: <gsl-discuss at sources dot redhat dot com>
- Date: Thu, 1 Aug 2002 11:19:14 -0500
- Subject: Questions on bsimp
Hi, there,
I just started using GSL for my research. I was wondering why for some
ODEs, rkf45 can give me very good results while the results from bsimp is
messy.
My code is as follows, the solution should be y1 = cos(x) and y2 = sin(x).
When I use bsimp as the stepper function, it seems it just jumps too each
step. I need bsimp because the speed is an important issue in my project.
The following code is just a test.
Thanks a lot!
Di
--------------------------------------------------------------
#include <stdio.h>
#include <math.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)
{
f[0] = -y[1];
f[1] = y[0];
return GSL_SUCCESS;
}
int jac (double t, const double y[], double *dfdy, double dfdt[],
void *params)
{
gsl_matrix_view dfdy_mat
= gsl_matrix_view_array (dfdy, 2, 2);
gsl_matrix * m = &dfdy_mat.matrix;
gsl_matrix_set (m, 0, 0, 0.0);
gsl_matrix_set (m, 0, 1, -1.0);
gsl_matrix_set (m, 1, 0, 1.0);
gsl_matrix_set (m, 1, 1, 0.0);
dfdt[0] = 0.0;
dfdt[1] = 0.0;
return GSL_SUCCESS;
}
int main (void)
{
const gsl_odeiv_step_type * T
= gsl_odeiv_step_rkf45;
gsl_odeiv_step * s
= gsl_odeiv_step_alloc (T, 2);
gsl_odeiv_control * c
= gsl_odeiv_control_standard_new(1e-30, 1e-10,
1.0, 1.0);
gsl_odeiv_evolve * e
= gsl_odeiv_evolve_alloc (2);
gsl_odeiv_system sys = {func, jac, 2, NULL};
double t = 0.0, t1 = 10.0;
double h = 1e-6;
double y[2] = {1.0, 0.0};
while (t < t1)
{
int status = gsl_odeiv_evolve_apply(e, c, s,
&sys,
&t, t1,
&h, y);
if (status != GSL_SUCCESS)
break;
printf ("%.5e %.5e %.5e\n", t, y[0], y[1]);
}
gsl_odeiv_evolve_free (e);
gsl_odeiv_control_free (c);
gsl_odeiv_step_free(s);
return 0;
}