ode.hpp

Go to the documentation of this file.

#pragma once
#include <dsp/common.hpp>
 
 
namespace rack {
namespace dsp {
 
 
template <typename T, typename F>
void stepEuler(T t, T dt, T x[], int len, F f) {
    T k[len];
 
    f(t, x, k);
    for (int i = 0; i < len; i++) {
        x[i] += dt * k[i];
    }
}
 
template <typename T, typename F>
void stepRK2(T t, T dt, T x[], int len, F f) {
    T k1[len];
    T k2[len];
    T yi[len];
 
    f(t, x, k1);
 
    for (int i = 0; i < len; i++) {
        yi[i] = x[i] + k1[i] * dt / T(2);
    }
    f(t + dt / T(2), yi, k2);
 
    for (int i = 0; i < len; i++) {
        x[i] += dt * k2[i];
    }
}
 
template <typename T, typename F>
void stepRK4(T t, T dt, T x[], int len, F f) {
    T k1[len];
    T k2[len];
    T k3[len];
    T k4[len];
    T yi[len];
 
    f(t, x, k1);
 
    for (int i = 0; i < len; i++) {
        yi[i] = x[i] + k1[i] * dt / T(2);
    }
    f(t + dt / T(2), yi, k2);
 
    for (int i = 0; i < len; i++) {
        yi[i] = x[i] + k2[i] * dt / T(2);
    }
    f(t + dt / T(2), yi, k3);
 
    for (int i = 0; i < len; i++) {
        yi[i] = x[i] + k3[i] * dt;
    }
    f(t + dt, yi, k4);
 
    for (int i = 0; i < len; i++) {
        x[i] += dt * (k1[i] + T(2) * k2[i] + T(2) * k3[i] + k4[i]) / T(6);
    }
}
 
 
} // namespace dsp
} // namespace rack