RungeKutta4.hpp

Go to the documentation of this file.

#ifndef QUANTUMDYNAMICS_SOLVERS_RUNGEKUTTA4_HPP
#define QUANTUMDYNAMICS_SOLVERS_RUNGEKUTTA4_HPP
#include <armadillo>

#include "../QuantumDynamics.hpp"
#include "QuantumDynamicSolver.hpp"

namespace solid
{

template <template <typename> class T1, typename T2, typename T3>
class RK4 : public IQuantumDynamicSolver<T1, T2, T3>
{
public:

    inline static constexpr auto label = "rk4";

    std::string getLabel() override
    {
        return label;
    }

    void Propagate(double time, double dtime, QuantumDynamics<T1, T2, T3> &qDynamics)
    {

        arma::cx_double I(0, 1);
        arma::cx_vec k1, k2, k3, k4;

        qDynamics.quantumSystem.Fill();
        k1 = -I * dtime * qDynamics.quantumSystem.hamiltonian.matrixElements * qDynamics.quantumState.vector;
        time += 0.5 * dtime;
        qDynamics.LoadParameters();
        qDynamics.quantumSystem.Fill();
        k2 = -I * dtime * qDynamics.quantumSystem.hamiltonian.matrixElements * (qDynamics.quantumState.vector + 0.5 * k1);
        k3 = -I * dtime * qDynamics.quantumSystem.hamiltonian.matrixElements * (qDynamics.quantumState.vector + 0.5 * k2);
        time += 0.5 * dtime;
        qDynamics.LoadParameters();
        qDynamics.quantumSystem.Fill();
        k4 = -I * dtime * qDynamics.quantumSystem.hamiltonian.matrixElements * (qDynamics.quantumState.vector + k3);
        qDynamics.quantumState.vector += 1. / 6. * (k1 + 2 * k2 + 2 * k3 + k4);
    }
};

} // namespace solid

#endif