• https://me.yahoo.com
COST (GBP)
6.00
0.00
0

# runge

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

C++

## Runge

 std::vectorrunge( double (*f)(double, double)[function pointer] double y0 double a double b double h )
Consider $\inline&space;\displaystyle&space;f:I&space;\times&space;\mathbb{R}&space;\rightarrow&space;\mathbb{R}$ a continuous function, where $\inline&space;\displaystyle&space;I$ is an interval and let $\inline&space;\displaystyle&space;y_0&space;\in&space;\mathbb{R}$ be the initial value in the Cauchy problem:

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

We need to find the function $\inline&space;\displaystyle&space;y:I&space;\rightarrow&space;\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}&space;=&space;y_k&space;+&space;\frac{1}{6}(\gamma_0&space;+&space;2\gamma_1&space;+&space;2\gamma_2&space;+&space;\gamma_3),&space;\qquad&space;k&space;\geq&space;0$

where

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

and $\inline&space;\displaystyle&space;x_{k+1}&space;-&space;x_k&space;=&space;h&space;>&space;0$. Hence the abscissas $\inline&space;x_0&space;<&space;x_1&space;<&space;x_2&space;<&space;\ldots&space;<&space;x_n$ divide a given interval $\inline&space;[a,&space;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&space;[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&space;[0,1]$ using a step of $\inline&space;h&space;=&space;0.1$:
$\left\{&space;\begin{array}{rcl}&space;\frac{dy}{dt}&space;&=&&space;g&space;-&space;1.3&space;y\\&space;\\&space;y(0)&space;&=&&space;3&space;\end{array}&space;\right&space;.$
which has the exact solution $\inline&space;\displaystyle&space;y(t)&space;=&space;\frac{g}{1.3}&space;+&space;(3&space;-&space;\frac{g}{1.3}){\rm&space;e}^{-1.3&space;t}$, with $\inline&space;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&space;y(0)&space;=&space;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&space;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.