A Toy Runge-Kutta Differential Equation Solver

This post presents a simple, interactive differential equation solution graphing tool based on the classic Runge-Kutta method (which is really an algorithm).

It graphs solutions to one 1^st-order equation, one 2^nd-order equation, or a system of two 1^st-order equations. The right-hand side of the equation(s) must be entered using syntax similar to expressions in the C, Java, or JavaScript programming languages (syntax help). Since the solution of a differential equation depends on starting conditions, you can set the range of starting conditions which should be graphed and the “time” value at which the starting conditions apply. The tool will draw one line for each set of starting conditions.

Hover over a solution line to see what starting conditions it is based on. Roll your mouse wheel to zoom in and out; hold down the middle mouse button to pan. On mobile, use one finger to pan and pinch with two fingers to zoom.

Click these buttons to see a variety of samples: y' = y y' = cos(y) y'' = cos(t) y'' = cos(y) y' = ay - by² (the logistic equation) y'' = y(a - 2q cos(t)) (Mathieu's equation) y' = -y²t² + (2y / t) (a Riccati equation) y'' = -(a / t)y' - (b / t²)y (a Cauchy-Euler equation) y'' = -y(α² + ω² sign(cos(t))) (Meissner equation) y' = y - (0.01yz), z' = -z + (0.02yz) (a “predator-prey” system)

Try playing with the parameters for some of the samples; you can get pretty wild pictures!

Classic Runge-Kutta, also known as “RK4”, generates four different estimates of the rate of change of each variable at each time step, and takes a weighted sum of those four estimates as a final estimate which is (usually) more accurate than any of the four. Then, we use those estimated derivatives to adjust the values of each variable, bump “time” forward by one step, and repeat until we reach the ending time of the simulation. Here is a simple implementation of RK4 for a system with just one dependent variable:

// Trace out evolution of our system using classic Runge-Kutta (AKA "RK4")
// Store results in a packed array of floats
function rk4trace(y, t, Δt, fn, array, i, Δi, limit) {
  while (i !== limit) {
    const half_Δt = Δt / 2.0;
    const next_t = t + Δt;
    const half_t = t + half_Δt;

    const k_1 = fn(t, y); // Slope at starting point
    const k_2 = fn(half_t, y + (half_Δt * k_1)); // Estimated slope at mid-point
    const k_3 = fn(half_t, y + (half_Δt * k_2)); // Another estimate of slope at mid-point
    const k_4 = fn(next_t, y + (Δt * k_3)); // Estimated slope at endpoint

    const slope = (k_1 + 2*k_2 + 2*k_3 + k_4) / 6.0; // Weighted average of those four slopes

    y += Δt * slope;
    array[i] = y;
    t += Δt;
    i += Δi;
  }
}

// Apply RK4 to find phase lines for a system with one dependent variable
function rk4solve(y_0, t_0, t_start, t_end, Δt, fn) {
  const timeSteps = Math.floor(((t_end - t_start) / Δt) + 1);
  // Packed array of variable values at each time step:
  const array = new Float64Array(timeSteps);

  let t = t_0, y = y_0, i = Math.floor((t_0 - t_start) / Δt);
  array[i] = y_0;

  // Trace out phase line from starting point
  rk4trace(y_0, t_0, Δt, fn, array, i+1, 1, array.length);

  // Trace out phase line in the opposite direction from the starting point
  rk4trace(y_0, t_0, -Δt, fn, array, i-1, -1, -1);

  return array;
}

The above implementation stores the values of the dependent variable in a Float64Array instead of a regular JavaScript Array; this is for speed and memory efficiency. It doesn't store the time value for each entry in the array, since that can be easily rederived from t_start, t_end, and Δt.

Expression Syntax Help for this Tool

Special thanks to Gilbert Strang for his text “Differential Equations and Linear Algebra”, which inspired me to make this tool.