Runge-Kutta метод

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

например, состоянием может быть набор координат и скоростей точек, тогда сложение и умножение на число делается с координатами и скоростями поэлементно, а производная по времени ставит скорости вместо координат и вычисленные из координат и скоростей ускорения - вместо скоростей. у меня состоянием был массив комплексных чисел - волновая ф-я

весь код получения состояния системы через шаг времени dt выглядит (на Elm'e) так:

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 ]

возможно операцию вычисления производной стоит не вносить в тайпкласс, а передавать явно. тогда можно для одного типа состояния считать его эволюцию при разных внешних условиях (разных методах вычисления сил/ускорений/операторов/whatever)


comments

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

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

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

тут RK4 же

Габриель Гонзалес, помнится, в своё время предлагал вместо тайпклассов носить руками словари

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

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

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