YANE-Framework Tutorial 1.1.0

Simple example of an ordinary differential equation

Introduction

This example solves the ordinary differential equation

\begin{eqnarray*} \ddot{x} & = & \lambda \cdot (1 - x ^2) \cdot \dot{x} - x \end{eqnarray*}

which is equivalent to the differential system

\begin{eqnarray*} \dot{x}_1(t) & = & x_2(t) \\ \dot{x}_2(t) & = & \lambda \cdot (1 - x_1(t) ^2) \cdot x_2(t) - x_1(t) \end{eqnarray*}

Lambda is a constant which changes the stiffness of the equation.

In this example we use the btmb::OdeSol3::Radau and btmb::OdeSol3::DoPri solver to solve the equation. We solve the equation with 3 different configurations to show the pros and cons of Radau and DoPri.

Code

First the code is explained to solve an equation with Radau or DoPri and after that we'll use some examples to explain the different results and runtimes.

Libraries

The first thing to do is to define the libraries which we require. These libraries form the head of our main program:

 #include <odesol3.h>
 #include <btmbutils2/exception.h>
 #include <btmbutils2/rtclock.h>

 #include <iostream>
 #include <cmath>
 #include <iomanip>
 #include <string>
 #include <cstdlib>
 #include <cstring>
 #include <fstream>

Namespaces

Having defined the libraries which we are going to use, we now have to define the namespaces in which the classes of these libraries operate:

 using namespace btmb::OdeSol3;
 using namespace std;

The equation

Every solver in the OdeSol class needs an equation. The solver calls this function by a function pointer to evaluate the equation. rpar and ipar are arrays that e.g. can be used to set and get constants. In this example we use rpar[0] as the lambda constant.

 void dgl ( int * n, double * t, double * x, double * dx, double * rpar, int * ipar )
 {
 dx[0] = x[1];
 dx[1] = rpar[0] * ( 1 - x[0] * x[0] ) * x[1] - x[0];
 }

Jacobian matrix

Radau needs the jacobian matrix to solve the equation. The radau solver is able to compute the matrix by its own but sometimes (e.g. very difficult equations) its better to write a function that sets the jacobian matrix.

 void jacobi ( int * n, double * t, double * x, double * dfx, int * ldfx, double * rpar, int* ipar )
 {
 dfx[0] = 0;
 dfx[1] = 1;
 dfx[2] = 2 * rpar[0] * x[0] - 1;
 dfx[3] = rpar[0] * ( 1 - x[0] * x[0] );
 }

Function to solve the equation

This function (for both, Radau and DoPri) iterates some steps to solve the equation. To solve an equation we need a solver and a function and a configuration for the solver. First we need to apply the function and the config to the solver by resetting the current data. Then we initialiase the solver with t0 (the initial time) and x0 (the initial value).

 void solve(OdeConfig * odeconfig, OdeSolveFirst * odesolver, OdeFunction * odefunc, double * x, double * rpar)
 {
 odesolver->reset ( odefunc, odeconfig );

 odesolver->init ( 0.0, x );

 cout << setw ( 7 ) << "t"  << " | " << setw ( 20 ) << "x" << endl;
 cout << "-------------------------------------" << fixed << setprecision ( 5 ) << endl;

 // Iteration loop
 }

After that we iterate some steps to solve the equation. An exception is thrown if the solver is unable to solve the equation, therefore we need try and catch. First we refresh the config and after that we call calc with the end time (i.e. the time up to which the equation is calculated) and the data array rpar.

 for ( ... )
 {
 try
 {
 odesolver->refreshConfig();
 odesolver->calc ( 0.2 * ( double ) i, rpar );
 cout << 0.2 * ( double ) i << " | ( " << setw ( 8 ) << x[0] << ", " << setw ( 8 ) << x[1] << " )" << endl;
 }
 catch ( btmb::Utils2::Exception e )
 {
 cout << "Unable to solve equation." << endl;
 }
 }

Compare function

In the main we create and delete a Radau and a DoPri solver, function and config. First we create a new clock object to . After that we create an odeconfig, and odesolver and an odefunc with a pointer to the equation above and the degree of the equation (or to be specific the dimension of the system of equations with degree one).

 btmb::Utils2::RTClock * clock = new btmb::Utils2::RTClock();
 clock->reset();

 OdeConfig * odeconfig = new DoPriConfig ( );
 OdeSolveFirst * odesolver = new DoPri5 ( );
 OdeFunction * odefunc = new OdeFunction ( dgl, 2 );

The solver overwrites the current x0 with the calculated x, therefore we need to create a temporary array. After that solve is called.

 double * tmp_x = new double[2];
 tmp_x[0] = x0[0]; tmp_x[1] = x0[1];

 solve(odeconfig, odesolver, odefunc, tmp_x, rpar);

 cout << "runtime using DoPri: " << clock->elapsedSeconds() << endl << endl;

After that we clean up.

 delete odeconfig;
 delete odesolver;
 delete odefunc;

Now the same with Radau, additionaly to the code above we set the Jacobi-Matrix.

 clock->reset();

 odeconfig = new RadauConfig ( );
 odesolver = new Radau5 ( );
 odefunc = new OdeFunction ( dgl, 2 );
 odefunc->setJacobiMatrix ( jacobi );

 ...

Now we do the same things we did above with Dopri the calculate to solution of the equation.

Main Program

In the main function we just need to create Arrays for lambda (rpar[0]) and x0, set the values and call compare.

 double * rpar = new double[1];
 double * x0 = new double[2];

 rpar[0] = 1.0;
 x0[0] = 2.0; x0[1] = 0.0;

 cout << "Minimal non-stiff example" << endl << endl;
 compare(x0, rpar);

Pros and cons

In this part we set various values of lambda.

Settings

In the first example we use lambda = 1 to get a non stiff problem. DoPri has problems with stiff problems, therefore we choose a much higher value (e.g. lambda = 5000) as second example. In the third example Dopri isn't able to solve the equation because the equation is too stiff (e.g. lambda = 10000).

Result

In the first example DoPri is a litte bit faster than Radau. In the second example the equation is quite stiff, but DoPri is able to calculate the solution but needs much longer than Radau. In the third example Dopri fails and Radau solves the equation without any problems.

Conclusion