CodeCogs

Computes an approximate solution to the Cauchy problem using the 4th order Runge-Kutta method.

Interface

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

using namespace Maths::Calculus::Ode;

std::vector<double> runge(double (*f)(double, double), double y0, double a, double b, double h)

Consider $\displaystyle f:I \times \mathbb{R} \rightarrow \mathbb{R}$ a continuous function, where $\displaystyle I$ is an interval and let $\displaystyle y_0 \in \mathbb{R}$ be the initial value in the Cauchy problem:

$\left\{ \begin{array}{rcl} \frac{dy}{dx} &=& f(x,y)\\ \\ y(x_0) &=& y_0, \qquad x_0 \in I. \end{array}\right .$

(1)

We need to find the function $\displaystyle y:I \rightarrow \mathbb{R}$ which satisfies the above conditions. It may happen that one cannot be able to compute the solution analytically, and in those cases numerical methods can give an approximate answer. This module approximates the solution to the Cauchy problem 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$

(2)

where

$\begin{array}{rcl} \gamma_0 &=& hf(x_k, y_k) \\ \\ \gamma_1 &=& hf\left(x_k + \frac{h}{2}, y_k + \frac{\gamma_0}{2}\right)\\ \\ \gamma_2 &=& hf\left(x_k + \frac{h}{2}, y_k + \frac{\gamma_1}{2}\right)\\ \\ \gamma_3 &=& hf(x_k + h, y_k + \gamma_2) \end{array}$

(3)

and $\displaystyle 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.

You may notice that the error increases with the index of the term in the recurrence relation, or as we get closer to the superior limit of the $[a,b]$ interval. However this module generally produces better estimates compared to euler.

Example 1

Next we give an example of how to use this function and display the absolute error (e) from the exact solution (Y). An approximate solution is found to the following Cauchy problem on the interval $[0,1]$ using a step of $h = 0.1$ :

$\left\{ \begin{array}{rcl} \frac{dy}{dt} &=& g - 1.3 y\\ \\ y(0) &=& 3 \end{array} \right .$

(4)

which has the exact solution $\displaystyle y(t) = \frac{g}{1.3} + (3 - \frac{g}{1.3}){\rm e}^{-1.3 t}$ , with $g$ the constant of gravitational acceleration. This is a mathematical model of the following situation: A person jumps out of an airplane at a speed of $y(0) = 3$ meters per second. Obviously he will feel the effects of gravity and will be accelerated downwards. At the same time he will feel the wind resistance which decreases his acceleration by a factor of $1.3$ up to the point where it becomes zero, when the velocity will become constant. The solution of the above differential equation then gives the velocity of the person at a precise moment of time in his free fall.

#include <codecogs/maths/calculus/ode/runge.h>
#include <stdio.h>
#include <math.h>
 
// precision constant
#define H  0.1
 
// initial value of the problem
#define Y0 3.0
 
// limits of the approximation interval
#define A  0.0
#define B  1.0
 
// the given function
double f(double x, double y)
{
  return 9.806650 - 1.3*y;
}

// the exact solution 
double exact(double x)
{
  return 9.80665/1.3 + (Y0 - 9.80665/1.3)*exp(-1.3*x);
}

int main()
{
  // compute the approximate solution
  std::vector<double> sol = Maths::Calculus::ODE::runge(f, Y0, A, B, H);
 
  // display the problem data
  printf("\n");
  printf("f(x, y) = g - 1.3*y\n");
  printf("     y0 = %.10lf\n\n", Y0);
  printf("      a = %.3lf\n", A);
  printf("      b = %.3lf\n", B);
  printf("      h = %.3lf\n\n", H);
 
  // display the results, including error estimation
  printf("Point       Approximation     Actual value     Error\n\n");

  // display the result
  for (int i = 0; i < sol.size(); i++)
    printf("x = %.1lf     %.11lf     %.11lf    %.11lf\n", 
    H*i + A, sol[i], exact(H*i + A), fabs(sol[i] - exact(H*i+A)));
 
  return 0;
}

Output

f(x, y) = g - 1.3*y
     y0 = 3.0000000000

      a = 0.000
      b = 1.000
      h = 0.100

Point       Approximation     Actual value     Error

x = 0.0     3.00000000000     3.00000000000    0.00000000000
x = 0.1     3.55388141096     3.55388278689    0.00000137593
x = 0.2     4.04024231492     4.04024473132    0.00000241639
x = 0.3     4.46731374976     4.46731693250    0.00000318274
x = 0.4     4.84232335469     4.84232708102    0.00000372633
x = 0.5     5.17161768890     5.17162177900    0.00000409009
x = 0.6     5.46076963892     5.46077394871    0.00000430979
x = 0.7     5.71467273264     5.71467714778    0.00000441514
x = 0.8     5.93762395601     5.93762838678    0.00000443076
x = 0.9     6.13339647409     6.13340085106    0.00000437696
x = 1.0     6.30530348701     6.30530775744    0.00000427043

Parameters

the function which describes the Cauchy problem

the initial value

the inferior limit of the interval

the superior limit of the interval

the precision constant (the step)

Returns

a vector containing approximate values of the solution at equally spaced abscissas

References

Mihai Postolache - "Metode Numerice", Editura Sirius