CodeCogs

Computes an approximate solution to a system of first-order ODE's using the 4th order Runge-Kutta method.

Interface

#include <codecogs/maths/calculus/ode/runge_sys.h>

using namespace Maths::Calculus::Ode;

std::vector<D_VECTOR> runge_sys(D_VECTOR (*f)(double, D_VECTOR), D_VECTOR Y0, double a, double b, double h)

FUNCTION

runge_sys

Consider $F:I \times \mathbb{R}^p \rightarrow \mathbb{R}^p$ a continuous function, where I is a real interval and p is a positive integer. This module solves the following system of ordinary differential equations:

$\left\{ \begin{matrix}{rcl} \frac{dY}{dx} &=& F(x,Y)\\ \\ Y(x_0) &=& y_0\\ \end{array}$

(1)

where $x_0 \in I$ , $y_0 \in \mathbb{R}^p$ and the unknown is the function $Y:\mathbb{R} \rightarrow \mathbb{R}^p$ , $Y = (Y_0, Y_1, \ldots, Y_{p-1})$ being composed of $p$ separate functions $Y_k:I \rightarrow \mathbb{R}$ , $0 \leq k \leq p - 1$ .

If F is written as $F = (F_0, F_1, \ldots, F_{p-1})$ being composed of p separate functions $F_k:I \times \mathbb{R}^p \rightarrow \mathbb{R}$ , $0 \leq k \leq p-1$ , then the above system can be rewritten equivalently as follows:

$\left\{ \begin{matrix}{rcl} \frac{dY_0}{dx} &=& F_0(x, Y)\\ \\ \frac{dY_1}{dx} &=& F_1(x, Y)\\ \\ \ldots & \ldots & \ldots\\ \\ \frac{dY_{p-1}}{dx} &=& F_{p-1}(x, Y)\\ \\ (Y_0, Y_1, \ldots, Y_{p-1})(0) &=& y_0 \end{array}$

(2)

This module approximates the solution of the above system at equally spaced abscissas using the recurrence relation:

$y_{k+1} = y_k + \frac{1}{6}(\Gamma_0 + 2\Gamma_1 + 2\Gamma_2 + \Gamma_3), \qquad k \geq 0$

(3)

where

$\begin{matrix}{rcl} \Gamma_0 &=& hF(x_k, y_k)\\ \\ \Gamma_1 &=& hF\left(x_k + \frac{1}{2}h, y_k + \frac{1}{2}\Gamma_0\right)\\ \\ \Gamma_2 &=& hF\left(x_k + \frac{1}{2}h, y_k + \frac{1}{2}\Gamma_1\right)\\ \\ \Gamma_3 &=& hF(x_k + h, y_k + \Gamma_2) \end{array}$

(4)

and $x_{k+1} - x_k = h > 0$ . Hence the abscissas $x_0 < x_1 < x_2 < \ldots < x_n$ divide a given interval $[a, b]$ in equal segments. The previous recurrence is known as the 4th order Runge-Kutta method for solving a system of first-order ODE's.

References

C-XSC, A C++ Class Library for Extended Scientific Computing, http://www.rz.uni-karlsruhe.de/~iam/html/language/cxsc/node14.html

Example 1

Next we give a fully documented example in which the solution of the following system is approximated on the interval $[0,1]$ using a step of $h = 0.1$ :

$\left\{ \begin{matrix}{rcl} \frac{dY}{dx} &=& (Y_1 Y_2, -Y_0 Y_2, -0.522\, Y_0 Y_1)\\ \\ Y(0) &=& (0,1,1) \end{array}$

(5)

which can be rewritten equivalently as:

$\left\{ \begin{matrix}{rcl} \frac{dY_0}{dx} &=& Y_1 Y_2\\ \\ \frac{dY_1}{dx} &=& -Y_0 Y_2\\ \\ \frac{dY_2}{dx} &=& -0.522\, Y_0 Y_1\\ \\ Y_0(0) &=& 0\\ Y_1(0) &=& 1\\ Y_2(0) &=& 1 \end{array}$

(6)

Notice that we have defined a new data type D_VECTOR to make it easier to declare variables or specify functions which are needed by this module.

#include <codecogs/maths/calculus/ode/runge_sys.h>
#include <stdio.h>
#include <math.h>

// precision constant
#define H  0.1

// limits of the approximation interval
#define A  0.0
#define B  1.0

// the given function
D_VECTOR f(double x, D_VECTOR Y)
{
  // in this case f is independent of x

  D_VECTOR Z;

  // the system of differential equations
  Z.push_back( Y[1]*Y[2] );
  Z.push_back( -Y[0]*Y[2] );
  Z.push_back( -0.522*Y[0]*Y[1] );

  return Z;
}

int main()
{
  // initial values vector
  D_VECTOR Y0;
  Y0.push_back( 0.0 );
  Y0.push_back( 1.0 );
  Y0.push_back( 1.0 );

  // compute the approximate solution
  std::vector<D_VECTOR> sol = 
  Maths::Calculus::ODE::runge_sys(f, Y0, A, B, H);

  // display the results
  for (int i = 0; i < sol.size(); i++)
    printf("x = %.2lf   Y[0] = %.10lf   Y[1] = %.10lf   Y[2] = %.10lf\n",
    H*i + A, sol[i][0], sol[i][1], sol[i][2]);

  return 0;
}

Output

x = 0.00   Y[0] = 0.0000000000   Y[1] = 1.0000000000   Y[2] = 1.0000000000
x = 0.10   Y[0] = 0.0997468763   Y[1] = 0.9950128339   Y[2] = 0.9973998103
x = 0.20   Y[0] = 0.1979929460   Y[1] = 0.9802034182   Y[2] = 0.9897155790
x = 0.30   Y[0] = 0.2933197040   Y[1] = 0.9560143620   Y[2] = 0.9772864890
x = 0.40   Y[0] = 0.3844633096   Y[1] = 0.9231402010   Y[2] = 0.9606466378
x = 0.50   Y[0] = 0.4703690538   Y[1] = 0.8824697046   Y[2] = 0.9404832898
x = 0.60   Y[0] = 0.5502233243   Y[1] = 0.8350174154   Y[2] = 0.9175872973
x = 0.70   Y[0] = 0.6234622825   Y[1] = 0.7818532761   Y[2] = 0.8928020595
x = 0.80   Y[0] = 0.6897598251   Y[1] = 0.7240380715   Y[2] = 0.8669765655
x = 0.90   Y[0] = 0.7489996600   Y[1] = 0.6625701827   Y[2] = 0.8409264404
x = 1.00   Y[0] = 0.8012372763   Y[1] = 0.5983465271   Y[2] = 0.8154050061

Parameters

the function which describes the system of equations

the initial values vector

the inferior limit of the interval

the superior limit of the interval

the precision constant (the step)

Returns

A matrix containing approximate values of the solution at equally spaced abscissas.