Решатели неавтономных линейных ODE или решатели ODE лиевых групп

Решатели неавтономных линейных уравнений ODE используются для уравнений в общем виде

и применяют структуру группы Ли решения для ускорения численных методов и записи определенных свойств в процессе решения. Одно из распространенных упрощений состоит в том, что решателям требуются независимые от состояния операторы, что предполагает следующий вид:

Другой тип решателя необходим, когда операторы зависят от состояния, т. е.

Для других специально требуется линейность, т. е.

где является константным оператором.

Стандартные интеграторы ODE

Стандартные интеграторы обыкновенных дифференциальных уравнений (ODE) будут работать с неавтономными линейными задачами ODE путем автоматической трансформации в ODE первого порядка. Более подробную информацию см. на странице решателей ODE.

Специализированные интеграторы OrdinaryDiffEq.jl

Если не указано иное, все алгоритмы OrdinaryDiffEq поставляются с эрмитовой полиномиальной интерполяцией третьего порядка. Алгоритмам, обозначенным как имеющие «свободную» интерполяцию, не требуются дополнительные шаги для интерполяции. Для функций более высокого порядка с несвободной интерполяцией дополнительные шаги выполняются «ленивым» способом (то есть не во время решения).

Обратите внимание, что все эти методы имеют фиксированный временной шаг, если не указано иное.

Экспоненциальные методы для линейных и аффинных задач

Эти методы требуют, чтобы был постоянным.

LinearExponential — точная формула решения линейных, не зависящих от времени задач.

Параметры:

krylov — символ. Один из следующих
- :off (по умолчанию) — предварительно кэширует оператор. Требуется метод Matrix(A), определенный для оператора A.
- :simple — использует простые аппроксимации Крылова с фиксированным размером подпространства m.
- :adaptive — использует адаптивные аппроксимации Крылова с внутренним переходом по шагам.
m — целое число, значение по умолчанию: 30. Управляет размером подпространства Крылова, если krylov=:simple, и размером начального подпространства, если krylov=:adaptive.
iop — целое число, значение по умолчанию: 0. Если не равно нулю, определяет длину неполной процедуры ортогонализации (IOP) ^[1]]. Заметим, что если линейный оператор или якобиан является эрмитовым, всегда будет использоваться алгоритм Ланцоша, а параметр IOP игнорируется.

using DifferentialEquations
_A = [2 -1; -3 -5] / 5
A = DiffEqArrayOperator(_A)
prob = ODEProblem(A, [1.0, -1.0], (1.0, 6.0))
sol = solve(prob, LinearExponential())

retcode: Success
Interpolation: 3rd order Hermite
t: 2-element Vector{Float64}:
 1.0
 6.0
u: 2-element Vector{Vector{Float64}}:
 [1.0, -1.0]
 [11.923375744981003, -4.833128734161089]

LinearExponential является точным и поэтому по умолчанию использует dt=tspan[2]-tspan[1]. Используемая интерполяция является неточной (эрмитовой третьего порядка). Таким образом, значения, создаваемые в ходе интерполяции (с помощью sol(t) или saveat), будут неточными с возрастающей погрешностью по мере увеличения размера временного интервала. Во избежание этой ситуации напрямую задайте dt или используйте tstops вместо saveat. Дополнительные сведения см. в описании этой проблемы.

Независимые от состояния решатели

Эти методы требуют, чтобы зависел только от независимой переменной, т. е. .

MagnusMidpoint — метод Магнуса средней точки второго порядка.
MagnusLeapfrog — метод Магнуса с перешагиванием второго порядка.
MagnusGauss4 — метод Магнуса четвертого порядка, аппроксимированный с помощью двухступенчатой квадратуры Гаусса.
MagnusGL4 — метод Магнуса четвертого порядка, аппроксимированный с помощью квадратуры Гаусса-Лежандра.
MagnusNC6 — метод Магнуса шестого порядка, аппроксимированный с помощью квадратуры Ньютона-Котеса.
MagnusGL6 — метод Магнуса шестого порядка, аппроксимированный с помощью квадратуры Гаусса-Лежандра.
MagnusNC8 — метод Магнуса восьмого порядка, аппроксимированный с помощью квадратуры Ньютона-Котеса.
MagnusGL8 — метод Магнуса восьмого порядка, аппроксимированный с помощью квадратуры Гаусса-Лежандра.

Пример:

function update_func(A, u, p, t)
    A[1, 1] = cos(t)
    A[2, 1] = sin(t)
    A[1, 2] = -sin(t)
    A[2, 2] = cos(t)
end
A = DiffEqArrayOperator(ones(2, 2), update_func = update_func)
prob = ODEProblem(A, ones(2), (1.0, 6.0))
sol = solve(prob, MagnusGL6(), dt = 1 / 10)

retcode: Success
Interpolation: 3rd order Hermite
t: 51-element Vector{Float64}:
 1.0
 1.1
 1.2000000000000002
 1.3000000000000003
 1.4000000000000004
 1.5000000000000004
 1.6000000000000005
 1.7000000000000006
 1.8000000000000007
 1.9000000000000008
 ⋮
 5.199999999999998
 5.299999999999998
 5.399999999999998
 5.499999999999997
 5.599999999999997
 5.699999999999997
 5.799999999999996
 5.899999999999996
 6.0
u: 51-element Vector{Vector{Float64}}:
 [1.0, 1.0]
 [0.9560322602141308, 1.138059339025732]
 [0.8837222756549167, 1.2712953135220313]
 [0.7836511349460551, 1.3924887039087297]
 [0.6585912321071351, 1.4945430086037101]
 [0.5134387044387002, 1.5712483084024886]
 [0.35484956094894243, 1.6179942444961344]
 [0.190617494783206, 1.6323079788211143]
 [0.028885769775000115, 1.6141203517260592]
 [-0.12267917859524063, 1.5657147744869344]
 ⋮
 [0.1649471518131617, 0.19050761906534755]
 [0.19016510949800752, 0.18488701621917145]
 [0.21691662630888853, 0.17939687440026209]
 [0.24554653627671552, 0.17419950752424784]
 [0.27642938801560896, 0.16953149318836871]
 [0.3099648774418323, 0.16572504624198953]
 [0.3465668460745136, 0.16323354028339193]
 [0.3866431012103112, 0.1626603704342466]
 [0.4305628420420221, 0.16478922184291447]

Начальные значения для в этом и подобных случаях не имеют значения, так как update_func сразу же их перезаписывает. С помощью ones(2,2) можно получить изменяемую матрицу 2x2.

Зависимые от состояния решатели

Эти методы можно использовать, когда зависит от переменных состояния, т. е. .

CayleyEuler — метод первого порядка, использующий преобразования Кэли.
LieEuler — метод Ли Эйлера первого порядка.
RKMK2 — метод Рунге — Кутты — Мунте-Кааса второго порядка.
RKMK4 — метод Рунге — Кутты — Мунте-Кааса четвертого порядка.
LieRK4 — метод Ли Рунге — Кутты четвертого порядка.
CG2 — метод Крауч — Гроссмана второго порядка.
MagnusAdapt4 — адаптивный метод Магнуса четвертого порядка.

Пример:

function update_func(A, u, p, t)
    A[1, 1] = 0
    A[2, 1] = sin(u[1])
    A[1, 2] = -1
    A[2, 2] = 0
end
A = DiffEqArrayOperator(ones(2, 2), update_func = update_func)
prob = ODEProblem(A, ones(2), (0, 30.0))
sol = solve(prob, LieRK4(), dt = 1 / 4)

retcode: Success
Interpolation: 3rd order Hermite
t: 121-element Vector{Float64}:
  0.0
  0.25
  0.5
  0.75
  1.0
  1.25
  1.5
  1.75
  2.0
  2.25
  ⋮
 28.0
 28.25
 28.5
 28.75
 29.0
 29.25
 29.5
 29.75
 30.0
u: 121-element Vector{Vector{Float64}}:
 [1.0, 1.0]
 [0.7273968579551914, 1.165794872747704]
 [0.4243104556165744, 1.2459210126860043]
 [0.10953160297897346, 1.2654872693089263]
 [-0.2070137554757648, 1.2681541526375768]
 [-0.5272224527692214, 1.3028146975278558]
 [-0.8648991354908154, 1.4149945285244698]
 [-1.244399456687535, 1.642767157085863]
 [-1.6975604074674677, 2.0038463754183535]
 [-2.2537711317460407, 2.449176910606262]
 ⋮
 [-2.6663888753291105, 2.6954837342889535]
 [-3.3611719087644736, 2.7811209839661157]
 [-4.008187561635523, 2.2967303664591197]
 [-4.470316530373341, 1.3459658518794193]
 [-4.66626299411022, 0.2081263142861722]
 [-4.57289678174422, -0.9494989022112654]
 [-4.200152270252206, -1.99585038589234]
 [-3.6052140155232593, -2.6767378307341736]
 [-2.9109618325260116, -2.781287635056373]

В приведенном примере решается нежесткое неавтономное линейное уравнение ODE с оператором, зависящим от состояния, с использованием метода LieRK4. Аналогично, жесткое неавтономное линейное уравнение ODE с операторами, зависящими от состояния, может быть решено с помощью специализированных адаптивных алгоритмов, например MagnusAdapt4.

Пример:

function update_func(A, u, p, t)
    A[1, 1] = 0
    A[2, 1] = 1
    A[1, 2] = -2 * (1 - cos(u[2]) - u[2] * sin(u[2]))
    A[2, 2] = 0
end
A = DiffEqArrayOperator(ones(2, 2), update_func = update_func)
prob = ODEProblem(A, ones(2), (30, 150.0))
sol = solve(prob, MagnusAdapt4())

retcode: Success
Interpolation: 3rd order Hermite
t: 592-element Vector{Float64}:
  30.0
  30.051554410975918
  30.14327232699745
  30.247524724052145
  30.375696909344892
  30.510855783890896
  30.65090569337366
  30.788186122035587
  30.91952329587291
  31.04078864675185
   ⋮
 147.96273686915632
 148.1386523776747
 148.34620600715417
 148.59559211110502
 148.8847559565441
 149.1635807149984
 149.44955972217943
 149.73645683847397
 150.0
u: 592-element Vector{Vector{Float64}}:
 [1.0, 1.0]
 [1.0418828143197263, 1.0526121015605732]
 [1.1297712835913454, 1.1520649699813994]
 [1.2527691803606753, 1.276033498183453]
 [1.4404455264543194, 1.4481983217458168]
 [1.6784196223822325, 1.6585827542112963]
 [1.941152341880844, 1.9121745307485782]
 [2.1365107433203687, 2.1934881055190636]
 [2.1250337632314107, 2.476421374656319]
 [1.7927262923353107, 2.7176941957849485]
 ⋮
 [1.3382785710099845, -1.3478235487981378]
 [1.1224319158310865, -1.132482377423453]
 [0.9547031935838253, -0.918358952842659]
 [0.8437842085519446, -0.6957609582428903]
 [0.7910831565788868, -0.4608017890713062]
 [0.7783134037066098, -0.24250843704692368]
 [0.7772133962093504, -0.020168014665712365]
 [0.7777533743732511, 0.20284549420577422]
 [0.7858794742201944, 0.4085646081595665]

Зависимые от времени и состояния операторы

Эти методы можно использовать, когда зависит от переменных времени и состояния, т. е. .

CG3 — метод Крауч — Гроссмана третьего порядка.

1. Описание IOP см. в этом документе.