Maths › Calculus › Ode ›

runge sys

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

Controller: CodeCogs

Contents

Interface
Runge Sys
Page Comments

Private project under development, to help contact the author:

Interface

C++

Runge Sys

std::vector<D_VECTOR>runge_sys(	D_VECTOR	(`f`)(double, D_VECTOR)[function pointer]*
	D_VECTOR	`Y0`
	double	`a`
	double	`b`
	double	`h`	)

Consider $\inline 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}$

(2)

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

If F is written as $\inline F = (F_0, F_1, \ldots, F_{p-1})$ being composed of p separate functions $\inline F_k:I \times \mathbb{R}^p \rightarrow \mathbb{R}$ , $\inline 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}$

(3)

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$

(4)

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}$

(5)

and $\inline x_{k+1} - x_k = h > 0$ . Hence the abscissas $\inline x_0 < x_1 < x_2 < \ldots < x_n$ divide a given interval $\inline [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 $\inline [0,1]$ using a step of $\inline 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}$

(6)

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}$

(7)

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

f	the function which describes the system of equations
Y0	the initial values vector
a	the inferior limit of the interval
b	the superior limit of the interval
h	the precision constant (the step)

Returns

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

Source Code

This module is private, for owner's use only.

Not a member, then Register with CodeCogs. Already a Member, then Login.