generic fourth-order Runge-Kutta method

[thedeemon]

Когда понадобилось численно просчитать эволюцию системы из ее производной по времени, и наивное решение не заработало, вспомнил слова "Рунге-Кутта", которые со студенческой скамьи не использовал. Быстро находимые описания в интернетах не супер понятные, а конкретные примеры порой слишком частные, но оказалось, что алгоритм можно очень просто выразить в довольно общем виде. Он работает для любого типа (представления состояния системы), для которого даны всего три операции: сложение, умножение на вещественное число и вычисление производной по времени. Такой вот тайпкласс. Например, состоянием может быть набор координат и скоростей точек, тогда сложение и умножение на число делается с координатами и скоростями поэлементно, а производная по времени ставит скорости вместо координат и вычисленные из координат и скоростей ускорения - вместо скоростей. У меня состоянием был массив комплексных чисел - волновая ф-я. Весь код получения состояния системы через шаг времени dt выглядит так:

type alias DiffVecSpace v = {
timeDerivative : v -> v,
mulByFloat : Float -> v -> v,
add : v -> v -> v
}

evolveRK : DiffVecSpace v -> Float -> v -> v
evolveRK ops dt state =
let a = ops.timeDerivative state
b = ops.add state (ops.mulByFloat (dt/2) a) |> ops.timeDerivative
c = ops.add state (ops.mulByFloat (dt/2) b) |> ops.timeDerivative
d = ops.add state (ops.mulByFloat dt c) |> ops.timeDerivative
in List.foldl ops.add state [ops.mulByFloat (dt/6) a,
ops.mulByFloat (dt/3) b,
ops.mulByFloat (dt/3) c,
ops.mulByFloat (dt/6) d]

(тут на Elm'e)
Если выразить на языке с нормальными тайпклассами или ООЯ с определяемыми операторами, получится еще проще. С другой стороны, возможно операцию вычисления производной стоит не вносить в тайпкласс, а передавать явно. Тогда можно для одного типа состояния считать его эволюцию при разных внешних условиях (разных методах вычисления сил/ускорений/операторов/whatever), что у меня активно используется в программе.