COST (GBP)

0.00

Maths › Calculus › Ode ›

runge

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

Controller: CodeCogs

Contents

Interface
Runge
Page Comments

Interface

C++

Runge

std::vector<double>runge(	double	(`f`)(double, double)[function pointer]*
	double	`y0`
	double	`a`
	double	`b`
	double	`h`	)

Consider $\inline \displaystyle f:I \times \mathbb{R} \rightarrow \mathbb{R}$ a continuous function, where $\inline \displaystyle I$ is an interval and let $\inline \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 .$

(2)

We need to find the function $\inline \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$

(3)

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

(4)

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

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 $\inline [a,b]$ interval. However this module generally produces better estimates compared to Maths/Calculus/ODE/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 $\inline [0,1]$ using a step of $\inline h = 0.1$ :

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

(5)

which has the exact solution $\inline \displaystyle y(t) = \frac{g}{1.3} + (3 - \frac{g}{1.3}){\rm e}^{-1.3 t}$ , with $\inline 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 $\inline 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 $\inline 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

References

Mihai Postolache - "Metode Numerice", Editura Sirius

Parameters

f	the function which describes the Cauchy problem
y0	the initial value
a	the inferior limit of the interval
b	the superior limit of the interval
h	the precision constant (the step)

Returns

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

Authors

Lucian Bentea (September 2006)

Source Code

Source code is available when you buy a Commercial licence.

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