GSL - GNU Scientific Library#

The GNU Scientific Library (GSL) is a numerical library for C and C++ programmers. It is free software under the GNU General Public License.

Link

Web site

Language

C

License

GPL

Documentation

documentation

Lorenz equations#

We would like to solve the Lorenz equations, a classical chaotic system given by:

\[\begin{split} \begin{cases} \dot{x} &= \sigma (y - x) \\ \dot{y} &= \rho x - y - x z \\ \dot{z} &= x y - \beta z \end{cases} \end{split}\]

with parameter \(\sigma=10\), \(\rho = 28\) and \(\beta = \frac{8}{3}\), and the initial state \((x_0, y_0, z_0) = (1,1,1)\). We solve it with a classical Runge-Kutta method of order 4, given by gsl_odeiv2_step_rk4 method.

To solve a problem with GSL we first write our equation as a ODE of the form:

\[ \dot{u} = f(t, u) \]

and the user provide the function \(f\) as:

  int f(double t, const double u[], double du[], void* params);

where u the current state of the function, t the current time and output du \(f(t,u)\). To use GSL you should convert your state type in C-style array.

 7int lorenz(double t, const double y[], double dy[], void *params)
 8{
 9  double sigma = ((double *)params)[0];
10  double rho = ((double *)params)[1];
11  double beta = ((double *)params)[2];
12
13  dy[0] = sigma * (y[1] - y[0]);
14  dy[1] = y[0] * (rho - y[2]) - y[1];
15  dy[2] = y[0] * y[1] - beta * y[2];
16
17  return GSL_SUCCESS;
18}

Now we should define the system and the driver for the chosen integrator

29  gsl_odeiv2_system lorenz_pb = {lorenz, NULL, 3, &params};
30  gsl_odeiv2_driver * meth = gsl_odeiv2_driver_alloc_y_new (&lorenz_pb, gsl_odeiv2_step_rk4, dt, 1e-6, 0.0);

now we can call the integrator by

53    int status = gsl_odeiv2_driver_apply(meth, &tn, tnp1, yn);

The complet time loop is

45  while( tn < tf )
46  {
47    sol[n][0] = tn;
48    sol[n][1] = yn[0];
49    sol[n][2] = yn[1];
50    sol[n][3] = yn[2];
51
52    double tnp1 = tn + dt;
53    int status = gsl_odeiv2_driver_apply(meth, &tn, tnp1, yn);
54    if (status != GSL_SUCCESS)
55    {
56      printf ("error, return value=%d\n", status);
57      break;
58    }
59
60    ++n;
61  }

where sol is an array where we save all states.

For the complet example, see lorenz.c source file.

Transport equation#

In this example we would like to solve the following PDE:

\[ \partial_t u + a \partial_x u = 0 \]

with \(t>0\), on the torus \(x\in[0, 1)\), a velocity \(a=1\) and the initial condition given by a hat function:

\[\begin{split} u(0, x) = \begin{cases} x - 0.25 & \text{if } x\in[0.25, 0.5[ \\ -x + 0.75 & \text{if } x\in[0.5, 0.75[ \\ 0 & \text{else} \end{cases} \end{split}\]

We choose a first order up-wind scheme to estimate the \(x\) derivative and a forward Euler method for the time discretization.

To define the explicit Euler method, we define a euler_state_t data structure to store current state and intermediate step:

typedef struct
{
    double* k;
    double* y0;
    double* y_onestep;
} euler_state_t;

and add a new gsl_odeiv2_step_type with following lines:

static gsl_odeiv2_step_type const euler_type = { "euler", // name
    1,                                                    // can use dydt_in
    1,                                                    // gives exact dydt_out
    &euler_alloc,
    &euler_apply,
    &stepper_set_driver_null,
    &euler_reset,
    &euler_order,
    &euler_free };

gsl_odeiv2_step_type const* gsl_odeiv2_step_euler = &euler_type;

where we should define how to allocate a euler_state_t (euler_alloc), increment it with euler_apply, reset data with euler_reset, a function to return the order euler_order and a deallocator euler_free.

We define the up-wind scheme as:

169int
170upwind( double t, double const y[], double dy[], void* params )
171{
172    double a  = ( (double*)params )[0];
173    double dx = ( (double*)params )[1];
174    int n_x   = ( (double*)params )[2];
175
176    dy[0] = -( fmax( a, 0. ) * ( y[0] - y[n_x - 1] ) + fmin( a, 0. ) * ( y[1] - y[0] ) ) / dx;
177
178    for ( int i = 1; i < n_x - 1; ++i )
179    {
180        dy[i] = -( fmax( a, 0. ) * ( y[i] - y[i - 1] ) + fmin( a, 0. ) * ( y[i + 1] - y[i] ) ) / dx;
181    }
182
183    dy[n_x - 1] = -( fmax( a, 0. ) * ( y[n_x - 1] - y[n_x - 2] ) + fmin( a, 0. ) * ( y[0] - y[n_x - 1] ) ) / dx;
184
185    return GSL_SUCCESS;
186}

The time loop is the same as for Lorenz equation.

For the complet example, see transport.c source file.

Arenstorf orbit#

The Arenstorf orbit problem is a classical problem to test adaptive time step methods:

\[\begin{split} \begin{cases} \ddot{x} &= x + 2\dot{y} - \frac{1-\mu}{r_1^3}(x+\mu) - \frac{\mu}{r_2^3}(x-1+\mu) \\ \ddot{y} &= y - 2\dot{x} - \frac{}{1-\mu}{r_1^3}y - \frac{\mu}{r_2^3}y \end{cases} \end{split}\]

with initial condition \((x,\dot{x},y,\dot{y})=(0.994, 0, 0, -2.001585106)\), \(r_1\) and \(r_2\) given by

\[ r_1 = \sqrt{(x+\mu)^2 + y^2},\quad r_2 = \sqrt{(x-1+\mu)^2 + y^2} \]

and with parameter \(\mu = 0.012277471\).

First of all, we need to rewrite this problem into a first order derivative equation in time

\[\begin{split} \begin{pmatrix} y_1 \\ y_2 \\ y_3 \\ y_4 \end{pmatrix} = \begin{pmatrix} x \\ y \\ \dot{x} \\ \dot{y} \end{pmatrix}, \qquad \begin{cases} \dot{y}_1 = y_3 \\ \dot{y}_2 = y_4 \\ \dot{y}_3 = y_1 + 2y_4 - \frac{1-\mu}{r_1^3}(y_1 + \mu) - \frac{\mu}{r_2^3}(y_1-1+\mu) \\ \dot{y}_4 = y_2 - 2y_3 - \frac{1-\mu}{r_1^3}y_2 - \frac{\mu}{r_2^3}y_2 \\ \end{cases} \end{split}\]

We define this system as:

 8int arenstorf(double t, const double y[], double dy[], void *params)
 9{
10  double mu = *(double *)params;
11
12  double const y1 = y[0];
13  double const y2 = y[1];
14  double const y3 = y[2];
15  double const y4 = y[3];
16
17  double const r1 = sqrt( ( y1 + mu ) * ( y1 + mu ) + y2 * y2 );
18  double const r2 = sqrt( ( y1 - 1. + mu ) * ( y1 - 1. + mu ) + y2 * y2 );
19
20  dy[0] = y3;
21  dy[1] = y4;
22  dy[2] = y1 + 2. * y4 - ( 1. - mu ) * ( y1 + mu ) / ( r1 * r1 * r1 ) - mu * ( y1 - 1. + mu ) / ( r2 * r2 * r2 );
23  dy[3] = y2 - 2. * y3 - ( 1. - mu ) * y2 / ( r1 * r1 * r1 ) - mu * y2 / ( r2 * r2 * r2 );
24
25  return GSL_SUCCESS;
26}

We solve this example with given method gsl_odeiv2_step_rk8pd which is the method RK8(7) 13M in [PD81] (mainly call DOPRI8). To do this we define our system to GSL and also a stepper, and some objects to control adaptive time step method:

37  gsl_odeiv2_system arenstorf_pb = {arenstorf, NULL, 4, &mu};
38  const gsl_odeiv2_step_type * meth = gsl_odeiv2_step_rk8pd;
39
40  gsl_odeiv2_step * s = gsl_odeiv2_step_alloc(meth, 4);
41  gsl_odeiv2_control * c = gsl_odeiv2_control_y_new(1e-5, 1e-5);
42  gsl_odeiv2_evolve * e = gsl_odeiv2_evolve_alloc(4);

then the time loop becomes

57  while( tn < tf )
58  {
59    sol[n][0] = tn;
60    sol[n][1] = yn[0];
61    sol[n][2] = yn[1];
62    sol[n][3] = yn[2];
63    sol[n][4] = yn[3];
64    sol[n][5] = dt;
65
66    int status = gsl_odeiv2_evolve_apply(e, c, s, &arenstorf_pb, &tn, tf, &dt, yn);
67
68    if (status != GSL_SUCCESS)
69    {
70      printf ("error, return value=%d\n", status);
71      break;
72    }
73
74    ++n;
75  }

we should call gsl_odeiv2_evolve_apply to only make a step and not call the method from initial time to final time.

For the complet example, see arenstorf.c source file.