Comparison of Euler Method, Runge-Kutta 2, Runge-Kutta 4, and Taylor Series Solutions of a Single Ordinary Differential Equation

#include "RungeKutta4.h"
#include <iomanip>
#include <iostream>
using namespace std;

int main()
{
	int nSteps;
	vector<Point> solution1, solution2;
	vector<Point> solution3, solution4;

	while (true)
	{
		cout << "n-steps (zero to exit app) = ";
		cin >> nSteps;

		if (nSteps == 0)
			break;

		cout << "n\t\txn\tyn2p(xn)\tyn4p(xn)\tEuler\t\tTaylor\t\tDFDX2\t\tDFDX4\r\n";
		RungeKutta4::ComputeEuler(nSteps, 1, 2, -1, solution1);
		RungeKutta4::ComputeRK2(nSteps, 1, 2, -1, solution2);
		RungeKutta4::ComputeRK4(nSteps, 1, 2, -1, solution3);
		RungeKutta4::ComputeTaylor(nSteps, 1, 2, -1, solution4);

		for (int i = 0; i <= nSteps; i++)
		{
			POINT pt1 = solution1[i];
			POINT pt2 = solution2[i];
			POINT pt3 = solution3[i];
			POINT pt4 = solution4[i];
			cout << setw(2) << pt1.n << "\t";
			cout << setprecision(8) << fixed << pt1.x << "\t";
			cout << setprecision(8) << fixed << pt3.y << "\t";
			cout << setprecision(8) << fixed << pt4.y << "\t";
			cout << setprecision(8) << fixed << pt1.y << "\t";
			cout << setprecision(8) << fixed << pt2.y << "\t";
			cout << setprecision(8) << fixed << pt2.derivative << "\t";
			cout << setprecision(8) << fixed << pt3.derivative << "\r\n";
		}
	}
}

#pragma once
#include <vector>
using namespace std;

typedef struct Point
{
	int n;
	long double derivative, x, y;
} POINT, *PPOINT;

class RungeKutta4
{
private:
	static long double FX(
		long double x,
		long double y);
	static long double DFDX(
		long double x,
		long double y,
		long double yp);
public:
	static void ComputeEuler(
		int nSteps,
		long double x0,
		long double x1,
		long double y0,
		vector<Point>& pts);
	static void ComputeRK2(
		int nSteps,
		long double x0,
		long double x1,
		long double y0,
		vector<Point>& pts);
	static void ComputeRK4(
		int nSteps,
		long double x0,
		long double x1,
		long double y0,
		vector<Point>& pts);
	static void ComputeTaylor(
		int nSteps,
		long double x0,
		long double x1,
		long double y0,
		vector<Point>& pts);
};

// Translated from FORTRAN 77 source code found
// in "Elementary Numerical Analysis:
// An Algorithmic Approach" by S. D. Conte and Carl de
// Boor. Translator: James Pate Williams, Jr. (c)
// August 20, 2023

#include "RungeKutta4.h"

long double RungeKutta4::FX(
	long double x, long double y)
{
	long double x2 = x * x;
	long double term1 = 1.0 / x2;
	long double term2 = -y / x;
	long double term3 = -y * y;
	return term1 + term2 + term3;
}

long double RungeKutta4::DFDX(
	long double x,
	long double y,
	long double yp)
{
	long double x3 = x * x * x;
	long double term1 = -2.0 / x3;
	long double term2 = -yp / x;
	long double term3 = y / x / x;
	long double term4 = -2.0 * y * yp;
	return term1 + term2 + term3 + term4;
}

void RungeKutta4::ComputeEuler(
	int nSteps,
	long double x0,
	long double x1,
	long double y0,
	vector<Point>& pts)
{
	Point pt0 = {};
	pt0.n = 0;
	pt0.x = x0;
	pt0.y = y0;
	double derivative = FX(x0, y0);
	pt0.derivative = derivative;
	pts.push_back(pt0);

	if (nSteps == 1)
		return;

	long double h = (x1 - x0) / nSteps;
	long double xn = x0;
	long double yn = y0;

	for (int n = 1; n <= nSteps; n++)
	{
		xn = x0 + n * h;
		long double yn = y0 + h * FX(xn, y0);
		derivative = FX(xn, y0);
		pt0.n = n;
		pt0.derivative = derivative;
		pt0.x = xn;
		pt0.y = yn;
		pts.push_back(pt0);
		y0 = yn;
	}
}

void RungeKutta4::ComputeRK2(
	int nSteps,
	long double x0,
	long double x1,
	long double y0,
	vector<Point>& pts)
{
	Point pt0 = {};
	pt0.n = 0;
	pt0.x = x0;
	pt0.y = y0;
	double derivative = FX(x0, y0);
	pt0.derivative = derivative;
	pts.push_back(pt0);

	if (nSteps == 1)
		return;

	long double h = (x1 - x0) / nSteps;
	long double xn = x0;
	long double yn = y0;

	for (int n = 1; n <= nSteps; n++)
	{
		long double k1 = h * FX(xn, yn);
		long double k2 = h * FX(xn + h, yn + k1);
		yn += 0.5 * (k1 + k2);
		xn = x0 + n * h;
		derivative = FX(xn, yn);
		pt0.n = n;
		pt0.derivative = derivative;
		pt0.x = xn;
		pt0.y = yn;
		pts.push_back(pt0);
	}
}

void RungeKutta4::ComputeRK4(
	int nSteps,
	long double x0,
	long double x1,
	long double y0,
	vector<Point>& pts)
{
	Point pt0 = {};
	pt0.n = 0;
	pt0.x = x0;
	pt0.y = y0;
	double derivative = FX(x0, y0);
	pt0.derivative = derivative;
	pts.push_back(pt0);

	if (nSteps == 1)
		return;

	long double h = (x1 - x0) / nSteps;
	long double xn = x0;
	long double yn = y0;

	for (int n = 1; n <= nSteps; n++)
	{
		long double k1 = h * FX(xn, yn);
		long double k2 = h * FX(xn + h / 2, yn + k1 / 2);
		long double k3 = h * FX(xn + h / 2, yn + k2 / 2);
		long double k4 = h * FX(xn + h, yn + k3);
		xn = x0 + n * h;
		yn = yn + (k1 + k2 + k3 + k4) / 6.0;
		derivative = FX(xn, yn);
		pt0.n = n;
		pt0.derivative = derivative;
		pt0.x = xn;
		pt0.y = yn;
		pts.push_back(pt0);
	}
}

void RungeKutta4::ComputeTaylor(
	int nSteps,
	long double x0,
	long double x1,
	long double y0,
	vector<Point>& pts)
{
	Point pt0 = {};
	pt0.n = 0;
	pt0.x = x0;
	pt0.y = y0;
	double derivative = FX(x0, y0);
	pt0.derivative = derivative;
	pts.push_back(pt0);

	if (nSteps == 1)
		return;

	long double h = (x1 - x0) / nSteps;
	long double xn = x0;
	long double yn = y0;
	long double yp = FX(x0, y0);

	for (int n = 1; n <= nSteps; n++)
	{
		xn = x0 + n * h;
		yp = FX(xn, y0);
		long double yn = y0 + h * (FX(xn, y0) + h * DFDX(xn, y0, yp) / 2);
		derivative = FX(xn, y0);
		pt0.n = n;
		pt0.derivative = derivative;
		pt0.x = xn;
		pt0.y = yn;
		pts.push_back(pt0);
		y0 = yn;
	}
}