Forgotten Password?

Or login with:

  • Facebook
  • Google
  • Yahoo
this unit 6.00
sub units 0.00


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




std::vector<double>rungedouble(*f)(double, double)[function pointer]
doubleh )
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:

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:


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: 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("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;
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


Mihai Postolache - "Metode Numerice", Editura Sirius


fthe function which describes the Cauchy problem
y0the initial value
athe inferior limit of the interval
bthe superior limit of the interval
hthe precision constant (the step)


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


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.

Last Modified: 3 Apr 10 @ 08:31     Page Rendered: 2022-03-14 15:58:32