differential_equation.cpp

#include <cmath>
#include <functional>
#include <iomanip>
#include <iostream>

using Point = std::pair<double, double>;
using Target = std::function<double(Point)>;

struct Solver {
  std::function<double(Point, double)> method;
  char const *name;
};

Solver euler(Target f) {
  return Solver{
      [=](auto p, double h) {
        auto [x, y] = p;
        return y + h * f(p);
      },
      "euler",
  };
}

Solver improved_euler(Target f) {
  return Solver{
      [=](auto p, double h) {
        auto [x, y] = p;
        return y + h * (f(p) + f({x + h, y + h * f(p)})) / 2;
      },
      "improved euler",
  };
}

Solver runge_kutta(Target f) {
  return Solver{
      [=](auto p, double h) {
        auto [x, y] = p;
        auto k1 = f(p), k2 = f({x + h / 2, y + h * k1 / 2}),
             k3 = f({x + h / 2, y + h * k2 / 2}), k4 = f({x + h, y + h * k3});
        return y + h * (k1 + 2 * k2 + 2 * k3 + k4) / 6;
      },
      "runge-kutta",
  };
}

template <class Solution> Solver exact(Solution solution) {
  return Solver{
      [=](auto p, double h) {
        auto [x, y] = p;
        return solution(x + h);
      },
      "exact",
  };
}

void solve(Solver &solver, Point init_p, double h) {
  Point p = init_p;
  while (p.first <= 3.0) {
    auto [x, y] = p;
    std::cout << x << "," << y << "\n";
    p.second = solver.method(p, h);
    p.first += h;
  }
}

int main() {
  auto const hs = {0.125, 0.5};
  Point const init_p = {1, 1};

  /*
    y' = 2y / x (x = 1, y = 1)
    y = x^2
  */
  Target target = [](auto p) {
    auto [x, y] = p;
    return 2 * y / x;
  };
  auto const solvers = {
      exact([](auto x) { return x * x; }),
      euler(target),
      improved_euler(target),
      runge_kutta(target),
  };

  std::cout << std::scientific;
  for (auto h : hs) {
    std::cout << "h=" << h << "\n";
    for (auto solver : solvers) {
      std::cout << solver.name << " method:\n";
      solve(solver, init_p, h);
    }
  }
}

error.cpp

#include <cmath>
#include <functional>
#include <iomanip>
#include <iostream>

using Point = std::pair<double, double>;
using Target = std::function<double(Point)>;

struct Solver {
  std::function<double(Point, double)> method;
  char const *name;
};

Solver euler(Target f) {
  return Solver{
      [=](auto p, double h) {
        auto [x, y] = p;
        return y + h * f(p);
      },
      "euler",
  };
}

Solver improved_euler(Target f) {
  return Solver{
      [=](auto p, double h) {
        auto [x, y] = p;
        return y + h * (f(p) + f({x + h, y + h * f(p)})) / 2;
      },
      "improved euler",
  };
}

Solver runge_kutta(Target f) {
  return Solver{
      [=](auto p, double h) {
        auto [x, y] = p;
        auto k1 = f(p), k2 = f({x + h / 2, y + h * k1 / 2}),
             k3 = f({x + h / 2, y + h * k2 / 2}), k4 = f({x + h, y + h * k3});
        return y + h * (k1 + 2 * k2 + 2 * k3 + k4) / 6;
      },
      "runge-kutta",
  };
}

double solve(Solver &solver, Point init_p, double h) {
  Point p = init_p;
  while (true) {
    auto [x, y] = p;
    if (3.0 < p.first + h) {
      return y;
    }
    p.second = solver.method(p, h);
    p.first += h;
  }
}

int main() {
  auto const hs = {0.125, 0.25, 0.5, 1.0};
  Point const init_p = {1, 1};

  /*
    y' = 2y / x (x = 1, y = 1)
    y = x^2
  */
  Target const target = [](auto p) {
    auto [x, y] = p;
    return 2 * y / x;
  };
  auto const exact = 3.0 * 3.0;
  auto const solvers = {
      euler(target),
      improved_euler(target),
      runge_kutta(target),
  };

  std::cout << std::scientific;
  for (auto solver : solvers) {
    std::cout << solver.name << " method:\n";
    std::cout << "h,error\n";
    for (auto h : hs) {
      std::cout << h << ",";
      auto approximate = solve(solver, init_p, h);
      std::cout << exact - approximate << "\n";
    }
  }
}