/* Runge-Kutta method of order 4 for first order ODEs
*
* Implemented in ANSI C.
*
* Author: Michael Mazack,
*
* Date: March 29th, 2006
*
* License: Public Domain. Redistribution and modification without
* restriction is granted.
*/
#include
/* Prototype for our ODE, dx/dt = f(t, x) */
double ode1(double, double);
/* Prototype for Runge-Kutta method of order 4 */
double rk4(double, double, double, double (*)(double, double));
int main()
{
/* Initial conditions:
* x(0) = 1
*/
double t = 0;
double x = 1;
/* Step size */
double h = 0.001;
/* Here we find all RK4 estimates on the interval [0,1] for t in
* intervals of h.
*/
for(t=0;t<=1;t+=h)
{
x = rk4(t, x, h, &ode1);
printf("%lf\t%lf\n", (t + h), x);
}
return 0;
}
/* Note: Any first order ODE can be used here.
* This function is only a sample.
*/
double ode1(double t, double x)
{
/* dx/dt = -3*x^2 + 2*t */
return (-3*x*x + 2*t);
}
/* Implementation of the Runge-Kutta method of order 4 for
* estimating differential equations.
*
* First argument: t (independent variable of the differential equation)
* Second argument: x (initial condition or previous evaluation of rk4())
* Third argument: h (step size)
* Fourth argument: (*f) (a function pointer to dx/dt = f(t, x))
*/
double rk4(double t, double x, double h, double (*f)(double, double))
{
/* Since this function will be called many many times, it
* is probably best to declare these as static.
*/
static double k1, k2, k3, k4;
/* This is one of the most elegant things I've ever seen. */
k1 = h*f(t, x);
k2 = h*f(t + 0.5*h, x + 0.5*k1);
k3 = h*f(t + 0.5*h, x + 0.5*k2);
k4 = h*f(t + h, x + k3);
return (x + (1.0/6.0)*(k1 + 2*k2 + 2*k3 + k4));
}