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), что у меня активно используется в программе.

Comments

[kosiakk.livejournal.com]

Там ещё можно считать ожидаемую ошибку вычисления и уменьшать dt если эта ошибка кажется слишком большой (и рассчитывать заново).

Тут запрограммирован RK-2,3 теоретически даёт качество третьего порядка o(dt^3) при количестве вычислений второго порядка.

Плюс ещё можно оптимизировать, что `a` на следующем шаге получает значение `d` из предыдущего успешного.

[thedeemon.livejournal.com]

Тут RK4 же.