Начало работы с дифференциальными уравнениями в Julia
В этом руководстве содержатся сведения о функциональных возможностях для решения обыкновенных дифференциальных уравнений (ОДУ). Кроме того, этот материал пошагово разбирается в видеоруководстве.
Пример 1. Решение скалярных уравнений
В этом примере мы будем решать уравнение
на временном интервале ], где . Здесь — это переменная текущего состояния, — переменная параметра (выражающая, например, скорость реакции или постоянную силы тяжести), а — текущее время.
(Известно, что аналитическое решение данного уравнения имеет вид u(t)=u₀\exp(αt)
, но мы используем DifferentialEquations.jl для решения данной задачи численным методом, которым решаются задачи, аналитическое решение которых неизвестно.)
В общем виде процесс состоит из определения задачи, ее решения и анализа результата. Полный код для решения этой задачи выглядит так:
using DifferentialEquations
f(u, p, t) = 1.01 * u
u0 = 1 / 2
tspan = (0.0, 1.0)
prob = ODEProblem(f, u0, tspan)
sol = solve(prob, Tsit5(), reltol = 1e-8, abstol = 1e-8)
using Plots
plot(sol, linewidth = 5, title = "Solution to the linear ODE with a thick line",
xaxis = "Time (t)", yaxis = "u(t) (in μm)", label = "My Thick Line!") # legend=false
plot!(sol.t, t -> 0.5 * exp(1.01t), lw = 3, ls = :dash, label = "True Solution!")
Его отдельные части рассматриваются ниже.
Шаг 1. Определение задачи
Для численного решения мы определяем тип задачи, задавая уравнение, начальное условие и временной интервал, на котором ищется решение:
using DifferentialEquations
f(u, p, t) = 1.01 * u
u0 = 1 / 2
tspan = (0.0, 1.0)
prob = ODEProblem(f, u0, tspan)
ODEProblem with uType Float64 and tType Float64. In-place: false
timespan: (0.0, 1.0)
u0: 0.5
Обратите внимание, что DifferentialEquations.jl выбирает типы для задачи на основе типов, которые использовались для определения задачи. В нашем примере u0
имеет тип Float64, поэтому при решении зависимые переменные будут иметь тип Float64. Так как tspan = (0.0,1.0)
— это кортеж значений типа Float64, при решении независимые переменные будут иметь тип Float64 (обратите внимание, что начальное и конечное время должно быть одного типа). Таким образом, для решения можно выбрать числа произвольной точности, числа с единицами измерения и т. д. Дополнительные примеры см. в руководствах на основе блокнотов.
Виды задач имеют и много других особенностей, включая возможность определять матрицы масс и содержать обратные вызовы для событий. Для каждого вида задачи имеется отдельная страница, на которой представлены конструктор и доступные поля. Для ОДУ страница находится здесь. Кроме того, пользователь может указать дополнительные связанные функции, чтобы ускорить работу решателей. Они подробно описываются на странице, посвященной перегрузкам производительности.
Шаг 2. Решение задачи
Управление решателями
После определения задачи она решается с помощью функции solve
.
sol = solve(prob)
retcode: Success
Interpolation: specialized 4th order "free" interpolation, specialized 2nd order "free" stiffness-aware interpolation
t: 5-element Vector{Float64}:
0.0
0.09964258706516003
0.3457024247583422
0.6776921908052249
1.0
u: 5-element Vector{Float64}:
0.5
0.552938681151017
0.7089376245893467
0.9913594502399238
1.3728004409033037
Решателями можно управлять с помощью параметров, описанных на странице руководства, посвященной общим параметрам решателей. Например, с помощью команды reltol
можно уменьшить относительную погрешность (чтобы получить более точный результат за счет увеличения количества временных шагов):
sol = solve(prob, reltol = 1e-6)
retcode: Success
Interpolation: specialized 7th order lazy interpolation, specialized 4rd order "free" stiffness-aware interpolation
t: 6-element Vector{Float64}:
0.0
0.08395921283331977
0.29459485215346937
0.5701065300378511
0.8686054589334371
1.0
u: 6-element Vector{Float64}:
0.5
0.5442490221301345
0.6732716570549966
0.889283154476355
1.2021893261622514
1.3728005075542231
Существует также много способов управлять выходными данными. Например, если задать saveat=0.1
, решатель будет сохранять каждую 0.1
временную точку. В сочетании с настройкой погрешности это выглядит так:
sol = solve(prob, reltol = 1e-6, saveat = 0.1)
retcode: Success
Interpolation: 1st order linear
t: 11-element Vector{Float64}:
0.0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1.0
u: 11-element Vector{Float64}:
0.5
0.553138320881742
0.6119240040558909
0.6769572322220972
0.7489019734808645
0.8284927602387266
0.9165421885254866
1.0139492143572488
1.1217083317951
1.240919726340709
1.3728005075542231
В общем случае значением параметра saveat
может быть любая коллекция временных точек, которые нужно сохранять. Обратите внимание, что при этом используются интерполяции для неограниченного временного шага, что позволяет ускорить решение. Кроме того, если важна только конечная точка, можно вообще отключить сохранение промежуточных конечных точек:
sol = solve(prob, reltol = 1e-6, save_everystep = false)
retcode: Success
Interpolation: 1st order linear
t: 2-element Vector{Float64}:
0.0
1.0
u: 2-element Vector{Float64}:
0.5
1.3728005075542231
В результате будет сохранена только конечная временная точка.
Выбор алгоритма решателя
В DifferentialEquations.jl есть метод для выбора алгоритма решателя по умолчанию. Он подбирает эффективный способ решения вашей задачи. Чтобы помочь пользователям в выборе подходящего алгоритма, DifferentialEquations.jl предлагает метод для выбора алгоритмов посредством подсказок. Метод выбора по умолчанию выбирает алгоритм исходя из точности числовых типов и именованных аргументов (например, погрешностей). Кроме того, можно указать alg_hints
для выбора подходящих значений по умолчанию на основе особенностей задачи и необходимых функций для решения. Например, если имеется жесткая задача и требуется решение высокой точности, но оптимальный жесткий алгоритм для нее неизвестен, можно использовать такое выражение:
sol = solve(prob, alg_hints = [:stiff], reltol = 1e-8, abstol = 1e-8)
retcode: Success
Interpolation: specialized 4rd order "free" stiffness-aware interpolation
t: 14-element Vector{Float64}:
0.0
0.012407826196308189
0.04943554396660186
0.09925889132467855
0.15847071468568555
0.2285503638356638
0.30734723422488125
0.3945147615835462
0.48838261679111594
0.5880567114707754
0.6923232042782961
0.800231609063282
0.9109623745010172
1.0
u: 14-element Vector{Float64}:
0.5
0.5063053789114713
0.5255987021154772
0.552724440835384
0.5867879554585146
0.62982623626939
0.6819994183534772
0.7447644602267912
0.8188284026934183
0.9055526107898983
1.006117894091279
1.1219707582414566
1.25473553046225
1.372800507495401
Используемый алгоритм можно также выбрать явным образом. DifferentialEquations.jl предлагает гораздо более широкий ассортимент алгоритмов решателей, чем традиционные библиотеки для решения дифференциальных уравнений. Многие из этих алгоритмов основаны на результатах недавних исследований и показывают более высокую эффективность по сравнению с «обычными» алгоритмами. Например, можно выбрать метод Цитураса 5-го порядка:
sol = solve(prob, Tsit5())
retcode: Success
Interpolation: specialized 4th order "free" interpolation
t: 5-element Vector{Float64}:
0.0
0.09964258706516003
0.3457024247583422
0.6776921908052249
1.0
u: 5-element Vector{Float64}:
0.5
0.552938681151017
0.7089376245893467
0.9913594502399238
1.3728004409033037
Обратите внимание, что параметры решателя можно использовать в сочетании с выбором алгоритма. Например, решить задачу с помощью Tsit5()
с более низкой погрешностью можно так:
sol = solve(prob, Tsit5(), reltol = 1e-8, abstol = 1e-8)
retcode: Success
Interpolation: specialized 4th order "free" interpolation
t: 17-element Vector{Float64}:
0.0
0.012407826196308189
0.042501278333560696
0.0817804940926822
0.12887385246570865
0.18409797482126694
0.24627458684331965
0.31479299098506436
0.38859636482100657
0.46686178626184555
0.5487161506239446
0.633434723668386
0.7203630271543068
0.8089580249850934
0.8987655385980757
0.9894162242718207
1.0
u: 17-element Vector{Float64}:
0.5
0.5063053789114713
0.5219304750950854
0.5430527156531716
0.5695067808068426
0.6021743690740426
0.6412025714747711
0.6871475333575261
0.7403258498418478
0.8012223528078949
0.8702768771165198
0.9480214886604926
1.0350186821189897
1.1319031558272872
1.239373504325514
1.3582039555461738
1.3728005076225747
Вот некоторые распространенные алгоритмы для решения ОДУ в DifferentialEquations.jl:
-
AutoTsit5(Rosenbrock23())
подходит для решения как жестких, так и нежестких уравнений. Это неплохой алгоритм в том случае, если об уравнении ничего не известно. -
AutoVern7(Rodas5())
подходит для эффективного решения как жестких, так и нежестких уравнений с высокой точностью. -
Tsit5()
подходит для стандартных нежестких уравнений. Это первый алгоритм, который следует попробовать в большинстве случаев. -
BS3()
подходит для быстрого решения нежестких уравнений с низкой точностью. -
Vern7()
подходит для решения нежестких уравнений с высокой точностью. -
Rodas4()
илиRodas5()
подходят для небольших жестких уравнений с определенными в Julia типами, событиями и т. д. -
KenCarp4()
илиTRBDF2()
подходят для решения жестких уравнений среднего размера (100—2000 ОДУ). -
RadauIIA5()
подходит для решения жестких уравнений с очень высокой точностью. -
QNDF()
подходит для решения больших жестких уравнений.
Исчерпывающий список доступных алгоритмов и подробные рекомендации см. в документации по решателям. Каждому виду задач посвящена отдельная страница, на которой представлены все подходящие решатели.
Шаг 3. Анализ решения
Обработка типа решения
Результатом функции solve
является объект решения. К 5-му значению решения можно получить доступ так:
sol[5]
0.5695067808068426
А получить время 8-го временного шага можно следующим образом:
sol.t[8]
0.31479299098506436
Также имеются вспомогательные функции. Мы можем построить массив путем включения по кортежам решений:
[t + u for (u, t) in tuples(sol)]
17-element Vector{Float64}:
0.5
0.5187132051077795
0.564431753428646
0.6248332097458538
0.6983806332725513
0.7862723438953096
0.8874771583180907
1.0019405243425905
1.1289222146628544
1.2680841390697404
1.4189930277404645
1.5814562123288787
1.7553817092732964
1.9408611808123806
2.13813904292359
2.3476201798179943
2.372800507622575
Либо в более общем виде:
[t + 2u for (u, t) in zip(sol.u, sol.t)]
17-element Vector{Float64}:
1.0
1.0250185840192507
1.0863622285237315
1.1678859253990255
1.2678874140793939
1.3884467129693523
1.5286797297928618
1.6890880577001166
1.8692480645047023
2.069306491877635
2.289269904856984
2.5294777009893714
2.790400391392286
3.0727643366396675
3.377512547249104
3.7058241353641685
3.7456010152451493
Такой вариант позволяет использовать больше частей типа решения. Возвращаемый объект по умолчанию представляет собой непрерывное решение благодаря интерполяции. К интерполированным значениям можно обращаться, используя sol
как функцию, например:
sol(0.45) # Значение решения при t=0,45
0.7876927465687831
Обратите внимание на следующее различие: при индексировании с помощью [i]
мы получаем значение на i
-м шаге, а при обращении по (t)
— интерполяцию для времени t
!
Если в решателе задан параметр dense=true
(это вариант по умолчанию, если не используется saveat
), то данная интерполяция представляет собой интерполяцию высокого порядка и поэтому обычно соответствует погрешности временных точек решения. Интерполяции, связанные с каждым решателем, подробно описываются на странице алгоритма решателя. Если задан параметр dense=false
(если он не указан явным образом, то действует только при save_everystep=false
или использовании saveat
), по умолчанию возвращается линейная интерполяция.
Дополнительные сведения об управлении выходными данными см. на странице, посвященной обработке решений.
Решения построения графиков
Хотя с помощью указанных выше инструментов можно напрямую строить графики временных точек решений, также определены вспомогательные команды посредством шаблонов для Plots.jl. Чтобы построить график объекта решения, просто вызовите функцию plot:
#]add Plots # Перед первым использованием пакета Plots.jl его необходимо установить!
using Plots
#plotly() # Вы также можете выбрать бэкенд построения графиков
plot(sol)
Результат функции plot можно форматировать с помощью атрибутов, доступных в Plots.jl. Дополнительные параметры, относящиеся к DiffEq, описываются на странице, посвященной построению графиков.
Например, на странице атрибутов Plots.jl указано, что толщину линии можно задать посредством аргумента linewidth
. Кроме того, можно задать заголовок с помощью title
. Добавим эти параметры в команду plot, чтобы получить нужный результат, настроим подписи осей и изменим условные обозначения (обратите внимание, что условные обозначения можно отключить с помощью параметра legend=false
), чтобы получить хорошо оформленный график:
plot(sol, linewidth = 5, title = "Solution to the linear ODE with a thick line",
xaxis = "Time (t)", yaxis = "u(t) (in μm)", label = "My Thick Line!") # legend=false
После этого можно добавлять информацию на график с помощью команды plot!
:
plot!(sol.t, t -> 0.5 * exp(1.01t), lw = 3, ls = :dash, label = "True Solution!")
Пример 2. Решение систем уравнений
В этом примере мы будем решать уравнения Лоренца:
Определение функции для решения ОДУ, обновляемой на месте, может быть выгодным с точки зрения производительности. Под этим имеется в виду, что вместо написания функции, возвращающей решение, вы создаете функцию, которая обновляет вектор, содержащий решение. Благодаря этому пакеты решателей DifferentialEquations.jl могут сократить объем памяти, выделяемой для массивов, и обеспечить более высокую производительность.
Для этого мы просто записываем результат в качестве первого элемента входных данных функции. Например, задача решения уравнений Лоренца определяется в виде такой функции:
function lorenz!(du, u, p, t)
du[1] = 10.0 * (u[2] - u[1])
du[2] = u[1] * (28.0 - u[3]) - u[2]
du[3] = u[1] * u[2] - (8 / 3) * u[3]
end
lorenz! (generic function with 1 method)
Затем эту функцию можно использовать для решения задачи:
using DifferentialEquations
u0 = [1.0; 0.0; 0.0]
tspan = (0.0, 100.0)
prob = ODEProblem(lorenz!, u0, tspan)
sol = solve(prob)
retcode: Success
Interpolation: specialized 4th order "free" interpolation, specialized 2nd order "free" stiffness-aware interpolation
t: 1263-element Vector{Float64}:
0.0
3.5678604836301404e-5
0.0003924646531993154
0.0032624077544510573
0.009058075635317072
0.01695646895607931
0.02768995855685593
0.04185635042021763
0.06024041165841079
0.08368541255159562
⋮
99.30760258626904
99.39665422328268
99.49536147459878
99.58822928767293
99.68983993598462
99.77864535713971
99.85744078539504
99.93773320913628
100.0
u: 1263-element Vector{Vector{Float64}}:
[1.0, 0.0, 0.0]
[0.9996434557625105, 0.0009988049817849058, 1.781434788799208e-8]
[0.9961045497425811, 0.010965399721242457, 2.146955365838907e-6]
[0.9693591634199452, 0.08977060667778931, 0.0001438018342266937]
[0.9242043615038835, 0.24228912482984957, 0.0010461623302512404]
[0.8800455868998046, 0.43873645009348244, 0.0034242593451028745]
[0.8483309847495312, 0.6915629321083602, 0.008487624590227805]
[0.8495036669651213, 1.0145426355349096, 0.01821208962127994]
[0.9139069574560097, 1.4425599806525806, 0.03669382197085303]
[1.088863826836895, 2.052326595543049, 0.0740257368585531]
⋮
[4.669609096878053, 3.061564434452441, 25.1424735017959]
[4.188801916573263, 4.617474401440693, 21.09864175382292]
[5.559603854699961, 7.905631612648314, 18.79323210016923]
[8.556629716266505, 12.533041060088328, 20.6623639692711]
[12.280585075547771, 14.505154761545633, 29.332088452699942]
[11.736883151600804, 8.279294641640229, 34.68007510231878]
[8.10973327066804, 3.2495066495235854, 31.97052076740117]
[4.958629886040755, 2.194919965065022, 26.948439650907677]
[3.8020065515435855, 2.787021797920187, 23.420567509786622]
С помощью инструментов для работы с шаблонами графиков, представленных на странице о построении графиков, можно выбрать трехмерный график Далитца для представления различных переменных:
using Plots
plot(sol, idxs = (1, 2, 3))
Обратите внимание, что график по умолчанию для многомерных систем представляет собой наложение каждого временного ряда. Мы можем отрисовать временной ряд только второго компонента, еще раз воспользовавшись интерфейсом выбора переменных:
plot(sol, idxs = (0, 2))
Обратите внимание, что здесь «переменная 0» соответствует независимой переменной («время»).
Определение параметризованных функций
Зачастую необходимо явно связать параметры с дифференциальными уравнениями. Это может потребоваться, например, в подпрограммах оценки параметров. В таком случае используются значения p
посредством следующего синтаксиса:
function parameterized_lorenz!(du, u, p, t)
du[1] = p[1] * (u[2] - u[1])
du[2] = u[1] * (p[2] - u[3]) - u[2]
du[3] = u[1] * u[2] - p[3] * u[3]
end
parameterized_lorenz! (generic function with 1 method)
а затем параметры добавляются в ODEProblem
:
u0 = [1.0, 0.0, 0.0]
tspan = (0.0, 1.0)
p = [10.0, 28.0, 8 / 3]
prob = ODEProblem(parameterized_lorenz!, u0, tspan, p)
ODEProblem with uType Vector{Float64} and tType Float64. In-place: true
timespan: (0.0, 1.0)
u0: 3-element Vector{Float64}:
1.0
0.0
0.0
Есть несколько приемов, позволяющих сделать функции красивее. Например:
function parameterized_lorenz!(du, u, p, t)
x, y, z = u
σ, ρ, β = p
du[1] = dx = σ * (y - x)
du[2] = dy = x * (ρ - z) - y
du[3] = dz = x * y - β * z
end
parameterized_lorenz! (generic function with 1 method)
Обратите внимание, что параметры p
могут быть любого типа: массивы, статические массивы, именованные кортежи и т. д. — любой тип, подходящий для вашей задачи.
Так как параметры существуют внутри функции, определенные таким образом функции также можно использовать для анализа чувствительности, подпрограмм оценки параметров и построения графиков бифуркации. Это делает DifferentialEquations.jl исчерпывающим решением для анализа дифференциальных уравнений, которое также обеспечивает высокую производительность.
Пример 3. Решение неоднородных уравнеий с помощью параметризованных функций
Параметризованные функции также можно использовать для построения неоднородных обыкновенных дифференциальных уравнений (которые также называются ОДУ с ненулевой правой частью). Они часто применяются для моделирования динамических систем с внешними входами (как правило, зависящими от времени). В качестве примера рассмотрим модель маятника, состоящего из тонкого стержня длиной l
и массой m
:
Здесь θ
и ω
— это соответственно угловое отклонение маятника от вертикального (отвесного) положения и угловая скорость, M
— момент внешних сил (возникающий, например, из-за действия ветра или двигателя), а g
— ускорение свободного падения.
using DifferentialEquations
using Plots
l = 1.0 # длина [м]
m = 1.0 # масса [кг]
g = 9.81 # ускорение свободного падения [м/с²]
function pendulum!(du, u, p, t)
du[1] = u[2] # θ'(t) = ω(t)
du[2] = -3g / (2l) * sin(u[1]) + 3 / (m * l^2) * p(t) # ω'(t) = -3g/(2l) sin θ(t) + 3/(ml^2)M(t)
end
θ₀ = 0.01 # начальное угловое отклонение [рад]
ω₀ = 0.0 # начальная угловая скорость [рад/с]
u₀ = [θ₀, ω₀] # начальный вектор состояния
tspan = (0.0, 10.0) # временной интервал
M = t -> 0.1sin(t) # момент внешних сил [Н·м]
prob = ODEProblem(pendulum!, u₀, tspan, M)
sol = solve(prob)
plot(sol, linewidth = 2, xaxis = "t", label = ["θ [rad]" "ω [rad/s]"], layout = (2, 1))
Обратите внимание на то, как момент внешних сил M
, зависящий от времени, вводится в функцию pendulum!
в качестве параметра. Ведь согласно общему правилу параметры могут быть любого типа. Здесь мы указываем, что параметр M
зависит от времени, представляя его в виде функции. Это выражается путем добавления зависимости от времени (t)
к имени параметра.
Кроме того, обратите внимание, что, в отличие от зависящего от времени параметра, вектор переменных состояния u
, который в общем случае также зависит от времени, используется без явного указания на зависимость от времени (t)
.
Пример 4. Использование других типов для систем уравнений
DifferentialEquations.jl поддерживает множество различных типов зависимых переменных (в общем случае должен подойти любой тип с линейным индексом). Поэтому вместо решения векторного уравнения давайте пусть u
будет матрицей! Для этого достаточно, чтобы переменная u0
была матрицей, и нужно определить функцию f
так, чтобы она принимала матрицу и выдавала матрицу. Матрицу линейных ОДУ можно определить следующим образом:
using DifferentialEquations
using Plots
A = [1.0 0 0 -5
4 -2 4 -3
-4 0 0 1
5 -2 2 3]
u0 = rand(4, 2)
tspan = (0.0, 1.0)
f(u, p, t) = A * u
prob = ODEProblem(f, u0, tspan)
ODEProblem with uType Matrix{Float64} and tType Float64. In-place: false
timespan: (0.0, 1.0)
u0: 4×2 Matrix{Float64}:
0.738679 0.732413
0.415883 0.479094
0.153059 0.458707
0.467867 0.34872
Здесь ОДУ представляет собой матрицу 4x2. Это линейная система, определяемая путем умножения на A
. Для решения ОДУ мы выполняем те же действия, что и ранее.
sol = solve(prob)
plot(sol)