rk4.cpp

// rk4.cpp
// Classical 4th-order Runge-Kutta ODE integrator.
//
// This file supplies the RungeKutta4() routine that Q4 (q4.cpp, q4_3body.cpp,
// trappist.cpp) call but which was missing from the project, so those programs
// could never link. It integrates the system y' = rhs(t, y) from t0 to tf using
// a fixed step h. Each Q4 source provides its own global rhs(t, y) (declared in
// cw2.hpp), which this integrator calls — so link exactly one Q4 source together
// with this file, e.g.:  g++ -I.. -I../Q1 Q4/q4.cpp rk4.cpp -o q4
//
// The Q4 mains call it as RungeKutta4(t, y, h, t + h), i.e. one step per call;
// the loop below also handles tf - t0 spanning several steps in general.

#include <valarray>

#include "cw2.hpp"

valarray<double> RungeKutta4(double t0, valarray<double> y0, double h, double tf) {
    valarray<double> y = y0;
    double t = t0;

    // March from t0 to tf in steps of h, trimming the final step so we land on
    // tf exactly and never overshoot.
    while (t < tf - 1.0e-12) {
        double step = h;
        if (t + step > tf) {
            step = tf - t;
        }

        valarray<double> k1 = rhs(t, y);
        valarray<double> k2 = rhs(t + 0.5 * step, y + 0.5 * step * k1);
        valarray<double> k3 = rhs(t + 0.5 * step, y + 0.5 * step * k2);
        valarray<double> k4 = rhs(t + step, y + step * k3);

        y += (step / 6.0) * (k1 + 2.0 * k2 + 2.0 * k3 + k4);
        t += step;
    }

    return y;
}

// printVec() is declared in cw2.hpp; provide it here for completeness so any
// translation unit that uses it links cleanly.
void printVec(const valarray<double> v) {
    cout << "[ ";
    for (size_t i = 0; i < v.size(); ++i) {
        cout << v[i] << ' ';
    }
    cout << "]" << endl;
}