RungeKutta.cpp
b21d1b487b69master

RungeKutta.cpp
View Options


			#include <iostream>
			#include "RungeKutta.h"

			void RungeKutta::step(Eigen::VectorXd &y, double t) {
			switch (methodOrder) {
			case 1 : // simple Euler step
			y += stepSize * (A * y + g(t));
			break;
			case 2 : // midpoint method
			y += stepSize * (A * (y + 0.5 * stepSize * (A * y + g(t))) + g(t + 0.5 * stepSize));
			break;
			case 3 : {// Kutta's third order method
			Eigen::VectorXd k1 = stepSize * (A * y + g(t));
			Eigen::VectorXd k2 = stepSize * (A * (y + k1 / 2.) + g(t + stepSize / 2.));
			Eigen::VectorXd k3 = stepSize * (A * (y - k1 + k2 / 2.) + g(t + stepSize));
			y += 1./6. * (k1 + 4 * k2 + k3);
			break;
			}
			case 4 : {//classical RK-4 method
			Eigen::VectorXd k1 = stepSize * (A * y + g(t));
			Eigen::VectorXd k2 = stepSize * (A * (y + k1 / 2.) + g(t + stepSize / 2.));
			Eigen::VectorXd k3 = stepSize * (A * (y + k2 / 2.) + g(t + stepSize / 2.));
			Eigen::VectorXd k4 = stepSize * (A * (y + k3) + g(t + stepSize));
			y += 1./6. * (k1 + 2 * k2 + 2 * k3 + k4);
			break;
			}
			default:
			break;

			}
			}

			void RungeKutta::setMethodOrder(int n) {
			if (n < 1 \|\| n > 4){
			throw std::invalid_argument("Runge-Kutta method is only implemented for order from 1, 2, 3 or 4.");
			} else {
			methodOrder = n;
			}
			}