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

[juan-gandhi.livejournal.com]

А действительно ж. Спасибо.

[zeit-raffer.livejournal.com]

лет 12 назад учил этому студентов. с тех пор не пользовался. хотя все могло сложиться и иначе.

[kosiakk.livejournal.com]

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

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

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

[thedeemon.livejournal.com]

Тут RK4 же.

[swizard.livejournal.com]

> Если выразить на языке с нормальными тайпклассами

Габриель Гонзалес, помнится, в своё время, наоборот, предлагал вместо тайпклассов носить руками словари: http://www.haskellforall.com/2012/05/scrap-your-type-classes.html

[thedeemon.livejournal.com]

А как он их в операторах вроде +, - и * будет передавать?
Хочется же умножение и сложение как + и * записывать, для чего и хочется хотя бы свои инстансы для number уметь определять.
А еще хочу умножать разные типы (число на вектор, матрицу на вектор и т.п.), т.е. наивных тайпклассов тоже недостаточно для щастья.

А Рунге с Куттом разве не муж и жена?

[nponeccop.livejournal.com]

> А еще хочу умножать

У Страуструпа эти мотивы тоже проглядывались в его ранних книгах. А до чего они его довели!

а не проще счтать в переменных мпульс-плотность?

[Сергей Макаров (from livejournal.com)]

почему расчет именно в переменных волновой функци, а не плотность мпульс? чем так лучше?

Re: а не проще счтать в переменных мпульс-плотность?

[thedeemon.livejournal.com]

Нагляднее просто, проще представить визуально. Мои глаза как-то не привыкли смотреть в пространстве импульсов. :)

Re: а не проще счтать в переменных мпульс-плотность?

[Сергей Макаров (from livejournal.com)]

я о навье стокс представлении.
пространство импульсов-плотностей....