// Two-tank Mixing Problem
// x'(t) = f(t, x, y)
// y'(t) = g(t, x, y)
// x(0) = 5, y(0) = 1
// See “Ordinary Differential Equations
// from Calculus to Dynamical Systems”
// by Virginia W. Noonburg for the exact
// solution, view pages 165 – 167. Also
// see “Elementary Numerical Analysis:
// An Algorithmic Approach Third
// Edition” by S. D. Conte and Carl de
// Boor pages 398 – 491 for the Runge-
// Kutta-4 equations.

#include "RK4System.h"
void RK4System::Solve(
real t0, real t1,
real x0, real z0,
real(*f)(real, real, real),
real(*g)(real, real, real),
int nSteps, vector<real>& tv,
vector<real>& xv, vector<real>& yv)
{
real k1, k2, k3, k4, l1, l2, l3, l4;
real h, tn, xn, yn, xn1, yn1;
h = (t1 - t0) / nSteps;
tn = t0;
xn = x0;
yn = z0;
tv.push_back(tn);
xv.push_back(xn);
yv.push_back(yn);
for (unsigned int n = 1; n <= nSteps; n++)
{
tn = t0 + n * h;
k1 = h * f(tn, xn, yn);
l1 = h * g(tn, xn, yn);
k2 = h * f(tn + 0.5 * h, xn + 0.5 * k1, yn + 0.5 * l1);
l2 = h * g(tn + 0.5 * h, xn + 0.5 * k1, yn + 0.5 * l1);
k3 = h * f(tn + 0.5 * h, xn + 0.5 * k2, yn + 0.5 * l2);
l3 = h * g(tn + 0.5 * h, xn + 0.5 * k2, yn + 0.5 * l2);
k4 = h * f(tn + h, xn + k3, yn + l3);
l4 = h * g(tn + h, xn + k3, yn + l3);
xn1 = xn + (k1 + 2 * k2 + 2 * k3 + k4) / 6.0;
yn1 = yn + (l1 + 2 * l2 + 2 * l3 + l4) / 6.0;
xn = xn1;
yn = yn1;
tv.push_back(tn);
xv.push_back(xn);
yv.push_back(yn);
}
}
#pragma once
// Solves the following system of first order
// ordinary differential equations. Formulas
// are from "Elementary Numerical Analysis:
// An Algorithmic Approach" by S. D. Conte &
// Carl de Boor (c) 1980 8.12 page 398.
// x' = f(t, x, y)
// y' = g(t, x, y)
#include <vector>
using namespace std;
typedef long double real;
class RK4System
{
public:
static void Solve(
real t0, real t1,
real x0, real y0,
real(*f)(real, real, real),
real(*g)(real, real, real),
int nSteps, vector<real>& tv,
vector<real>& xv, vector<real>& yv);
};