Задачи ODE

# SciMLBase.ODEProblem — Type

Определяет задачу ODE. Страница документации: https://docs.sciml.ai/DiffEqDocs/stable/types/ode_types/

Математическая спецификация задачи ОДУ

Для определения задачи ОДУ достаточно задать функцию и начальное условие , которые определяют ОДУ:

Задать f можно сделать двумя разными способами.

f(du,u,p,t): на месте. Эффективное использование памяти при отказе от выделения. Лучший вариант для большинства случаев, если изменение запрещено.
f(u,p,t): возврат du. Менее эффективный способ использования памяти, подходящий при запрете изменений (например, с некоторыми пакетами автоматической дифференциации, такими как Zygote).

Условие u₀ должно представлять собой массив AbstractArray (или число), геометрия которого соответствует требуемой геометрии u. Обратите внимание, что тип значения u₀ не ограничен числами или векторами; значением u₀ также могут быть произвольные матрицы или тензоры высокой разрядности.

Сведения о матрице масс см. в документации по ODEFunction.

Тип задачи

Конструкторы

ODEProblem можно построить, предварительно создав ODEFunction, или просто передав конструктору правую часть ОДУ. Доступные конструкторы:

ODEProblem(f::ODEFunction,u0,tspan,p=NullParameters();kwargs...)
ODEProblem{isinplace,specialize}(f,u0,tspan,p=NullParameters();kwargs...): Определяет ОДУ с указанными функциями. isinplace дополнительно задает, является ли функция функцией на месте или нет. Это определяется автоматически, а не выводится. specialize дополнительно определяет уровень специализации. Дополнительные сведения см. в разделе об уровнях специализации в документации по SciMLBase. По умолчанию используется AutoSpecialize.

Дополнительные сведения о параметрах на месте и параметрах специализации см. в документации по ODEFunction.

Параметры являются необязательными, и если они не заданы, будет использоваться одинарный элемент NullParameters(), который при попытке индексирования несуществующих параметров будет выдавать понятные ошибки. Все дополнительные именованные аргументы передаются в решатели. Например, если в задаче задать обратный вызов (callback), этот callback будет добавляться при каждом вызове функции solve.

Сведения о задании якобианов и матриц масс см. в документации по ODEFunction.

Поля

f: функция в ОДУ.
u0: начальное условие.
tspan: временной интервал для задачи.
p: параметры.
kwargs: именованные аргументы, передаваемые в решатели.

Пример задачи

using SciMLBase
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
u0 = [1.0;0.0;0.0]
tspan = (0.0,100.0)
prob = ODEProblem(lorenz!,u0,tspan)

# Проверка, что все работает
using OrdinaryDiffEq
sol = solve(prob,Tsit5())
using Plots; plot(sol,vars=(1,2,3))

Дополнительные примеры задач

Примеры задач можно найти на странице DiffEqProblemLibrary.jl.

Использовать пример задачи, например prob_ode_linear, можно так:

#] add ODEProblemLibrary
using ODEProblemLibrary
prob = ODEProblemLibrary.prob_ode_linear
sol = solve(prob)

# SciMLBase.ODEFunction — Type

ODEFunction{iip,F,TMM,Ta,Tt,TJ,JVP,VJP,JP,SP,TW,TWt,TPJ,S,S2,S3,O,TCV} <: AbstractODEFunction{iip,specialize}

Представление функции ОДУ f, определяемое следующим образом:

и все связанные функции, такие как якобиан f, градиент относительно времени и др. Для всех случаев u0 — это начальное условие, p — это параметры, а t — это независимая переменная.

Конструктор

ODEFunction{iip,specialize}(f;
                           mass_matrix = __has_mass_matrix(f) ? f.mass_matrix : I,
                           analytic = __has_analytic(f) ? f.analytic : nothing,
                           tgrad= __has_tgrad(f) ? f.tgrad : nothing,
                           jac = __has_jac(f) ? f.jac : nothing,
                           jvp = __has_jvp(f) ? f.jvp : nothing,
                           vjp = __has_vjp(f) ? f.vjp : nothing,
                           jac_prototype = __has_jac_prototype(f) ? f.jac_prototype : nothing,
                           sparsity = __has_sparsity(f) ? f.sparsity : jac_prototype,
                           paramjac = __has_paramjac(f) ? f.paramjac : nothing,
                           syms = __has_syms(f) ? f.syms : nothing,
                           indepsym= __has_indepsym(f) ? f.indepsym : nothing,
                           paramsyms = __has_paramsyms(f) ? f.paramsyms : nothing,
                           colorvec = __has_colorvec(f) ? f.colorvec : nothing,
                           sys = __has_sys(f) ? f.sys : nothing)

Обратите внимание, что требуется только сама функция f. Эта функция должна быть задана как f!(du,u,p,t) или du = f(u,p,t). Более подробные сведения об обработке на месте и не на месте см. в разделе об iip.

Все остальные функции являются необязательными для улучшения или ускорения использования f. К ним относятся следующие.

mass_matrix: матрица масс M, представленная в функции ОДУ. С ее помощью можно определить, что уравнение является дифференциально-алгебраическим уравнением (ДАУ), если M сингулярна. Обратите внимание, что в этом случае требуются специальные решатели. Подробнее см. на странице о решателях ДАУ: https://docs.sciml.ai/DiffEqDocs/stable/solvers/dae_solve/. Должен быть AbstractArray или AbstractSciMLOperator.
analytic(u0,p,t): используется для передачи функции аналитического решения для аналитического решения уравнения ОДУ. Обычно применяется только для тестирования и разработки решателей.
tgrad(dT,u,p,t) или dT=tgrad(u,p,t): возвращает
jac(J,u,p,t) или J=jac(u,p,t): возвращает
jvp(Jv,v,u,p,t) или Jv=jvp(v,u,p,t): возвращает направленную производную
vjp(Jv,v,u,p,t) или Jv=vjp(v,u,p,t): возвращает сопряженную производную
jac_prototype: матрица-прототип, соответствующая типу, который соответствует якобиану. Например, если якобиан является трехдиагональным, в качестве прототипа может использоваться матрица Tridiagonal соответствующего размера, и по возможности интеграторы будут специализировать эту структуру. В неструктурированных шаблонах разреженности для якобиана следует использовать SparseMatrixCSC с правильным шаблоном разреженности. По умолчанию используется nothing, что означает плотный якобиан.
paramjac(pJ,u,p,t): возвращает якобиан параметра .
syms: имена символов для элементов уравнения. По размеру они должны соответствовать u0. Например, если u0 = [0.0,1.0] и syms = [:x, :y], к значениям будет применено каноническое именование, что позволит использовать sol[:x] в решении и автоматически присваивать имена значениям на графиках.
indepsym: каноническое именование для независимой переменной. Значение по умолчанию — nothing, которое внутренним образом использует t в качестве представления в любых графиках.
paramsyms: имена символов для параметров уравнения. По размеру они должны соответствовать p. Например, если p = [0.0, 1.0] и paramsyms = [:a, :b], к значениям будет применено каноническое именование, что позволит использовать sol[:a] в решении.
colorvec: цветовой вектор в соответствии с определением SparseDiffTools.jl для шаблона разреженности матрицы-прототипа jac_prototype. Это позволяет ускорить построение якобиана при использовании конечных разностей и автоматического дифференцирования на основе шаблона разреженности. Значение по умолчанию — nothing, что означает, что цветовой вектор будет внутренним образом вычисляться по требованию, когда это необходимо. Стоимость этой операции существенно зависит от шаблона разреженности.

iip: на месте и не на месте

iip является необязательным логическим значением, определяющим, написана ли данная функция для использования на месте или не на месте. Функциями на месте являются f!(du,u,p,t), в которых возврат игнорируется, и ожидается, что результат будет изменен на значение du. Функциями не на месте являются du=f(u,p,t).

Обычно это определяется автоматически путем просмотра таблицы методов для f и определения максимального количества аргументов в доступных диспетчеризациях. По этой причине конструктор ODEFunction(f) в целом работает (но нестабилен по типу). Однако для обеспечения стабильности типа или корректности этот параметр передается через ODEFunction{true}(f).

specialize: управление компиляцией и специализацией

Параметр specialize управляет уровнем специализации ODEFunction для функции f. Это позволяет найти компромисс между производительностью во время компиляции и выполнения. Доступны следующие уровни специализации:

SciMLBase.AutoSpecialize: эта форма выполняет отложенное заключение функций ODE в оболочку, чтобы остановить перекомпиляцию решателя ODE, но разрешить prob.f оставаться без оболочки для обычного использования. Этот уровень специализации используется по умолчанию и обеспечивает баланс между производительностью во время компиляции и во время выполнения.
SciMLBase.FullSpecialize: эта форма полностью специализирует ODEFunction по составляющим ее поля функциям. Таким образом, каждая функция ODEFunction в этой форме уникально типизирована, требуя переспециализации и компиляции для каждого нового определения ODE. Эта форма имеет наибольшее время компиляции за счет того, что является наиболее оптимальной во время выполнения. Ее рекомендуется использовать для длительных вычислений (например, внутри циклов оптимизации) и для тестов производительности.
SciMLBase.NoSpecialize: эта форма полностью отменяет специализацию типов функций в определении ODEFunction с помощью объявления типов Any. В итоге это может привести к снижению производительности во время выполнения, однако такая форма обеспечивает минимальное время компиляции.
SciMLBase.FunctionWrapperSpecialize: безотложная форма заключения функции в оболочку. Она небезопасна для многих решателей и поэтому используется в основном для тестирования при разработке.

Дополнительные сведения см. в разделе об уровнях специализации в документации по SciMLBase.

Поля

Поля типа ODEFunction напрямую соответствуют именам входных данных.

Дополнительные сведения о якобианах

В следующем примере создается функция ODEFunction не месте, якобианом которой является Diagonal.

using LinearAlgebra
f = (du,u,p,t) -> du .= t .* u
jac = (J,u,p,t) -> (J[1,1] = t; J[2,2] = t; J)
jp = Diagonal(zeros(2))
fun = ODEFunction(f; jac=jac, jac_prototype=jp)

Обратите внимание, что интеграторы всегда будут создавать глубокую копию fun.jac_prototype, поэтому можно не опасаться назначения псевдонимов.

В общем случае прототипом якобиана может быть все, для чего определена функция mul!. В частности, разреженные матрицы или пользовательские отложенные типы, поддерживающие mul!. Особый случай: когда jac_prototype является оператором AbstractSciMLOperator, не нужно указывать jac, так как для него автоматически задается update_coefficients!. Более подробную информацию о настройке операторов, зависящих от времени или параметров, см. в документации по AbstractSciMLOperator.

Примеры

Объявление явных якобианов для ОДУ

В наиболее стандартном случае объявление функции для якобиана осуществляется путем перегрузки функции f(du,u,p,t) обновляемой на месте функцией для якобиана f_jac(J,u,p,t), где для диспетчеризации используется тип значения. Например, возьмем модель Лотки-Вольтерры:

function f(du,u,p,t)
  du[1] = 2.0 * u[1] - 1.2 * u[1]*u[2]
  du[2] = -3 * u[2] + u[1]*u[2]
end

Для объявления якобиана просто добавим диспетчеризацию:

function f_jac(J,u,p,t)
  J[1,1] = 2.0 - 1.2 * u[2]
  J[1,2] = -1.2 * u[1]
  J[2,1] = 1 * u[2]
  J[2,2] = -3 + u[1]
  nothing
end

Тогда можно указать якобиан для ОДУ в виде:

ff = ODEFunction(f;jac=f_jac)

и использовать это в ODEProblem:

prob = ODEProblem(ff,ones(2),(0.0,10.0))

Символическое построение функций

См. описание функции modelingtoolkitize из ModelingToolkit.jl для автоматического символического построения якобиана и других компонентов из численно-определенных функций.

Тип решения

# SciMLBase.ODESolution — Type

struct ODESolution{T, N, uType, uType2, DType, tType, rateType, P, A, IType, S, AC<:Union{Nothing, Vector{Int64}}} <: SciMLBase.AbstractODESolution{T, N, uType}

Представление решения для обычного дифференциального уравнения, определяемого задачей ODEProblem.

Интерфейс DESolution

Более подробную информацию о взаимодействии с типами DESolution можно найти на странице о работе с решениями в документации по DifferentialEquations.jl.

https://docs.sciml.ai/DiffEqDocs/stable/basics/solution/

Поля

u: представление решения ОДУ. Задается в виде массива решений, где u[i] соответствует решению в момент времени t[i]. В большинстве случаев рекомендуется не обращаться к sol.u напрямую, а использовать интерфейс массива, описанный на странице о работе с решениями в документации по DifferentialEquations.jl.
t: временные точки, соответствующие сохраненным значениям решения ОДУ.
prob: исходная задача ODEProblem, которая была решена.
alg: тип алгоритма, используемого решателем.
stats: статистика по решателю, например количество необходимых вычислений функции, количество вычисленных якобианов и др.
retcode: код возврата от решателя. Используется для определения того, успешно ли решен решатель, завершился ли он досрочно из-за определенного пользователем обратного вызова или завершился ли он из-за ошибки. Дополнительные сведения см. в документации по кодам возврата.

Примеры задач

Примеры задач можно найти на странице DiffEqProblemLibrary.jl.

Использовать пример задачи, например prob_ode_linear, можно так:

#] add DiffEqProblemLibrary
using DiffEqProblemLibrary.ODEProblemLibrary
# загружаем задачи
ODEProblemLibrary.importodeproblems()
prob = ODEProblemLibrary.prob_ode_linear
sol = solve(prob)

# ODEProblemLibrary.prob_ode_linear — Constant

Линейное ODE

с начальным условием , и решением

с Float64s. Параметром является

# ODEProblemLibrary.prob_ode_2Dlinear — Constant

4x2 версия линейного ODE

с начальным условием как все равномерно распределенные случайные числа, и решением

с Float64s

# ODEProblemLibrary.prob_ode_bigfloatlinear — Constant

Линейное ODE

с начальным условием , и решением

с BigFloats

# ODEProblemLibrary.prob_ode_bigfloat2Dlinear — Constant

4x2 версия линейного ODE

с начальным условием как все равномерно распределенные случайные числа, и решением

с BigFloats

# ODEProblemLibrary.prob_ode_large2Dlinear — Constant

100x100 версия линейного ODE

с начальным условием как все равномерно распределенные случайные числа, и решением

с Float64s

# ODEProblemLibrary.prob_ode_2Dlinear_notinplace — Constant

4x2 версия линейного ODE

с начальным условием как все равномерно распределенные случайные числа, и решением

с Float64. Преднамеренно отсутствует в качестве теста.

# ODEProblemLibrary.prob_ode_lotkavolterra — Constant

Уравнения Lotka-Volterra (нежесткие)

с начальным условием

# ODEProblemLibrary.prob_ode_fitzhughnagumo — Constant

Fitzhugh-Nagumo (нежесткие)

с начальным условием

# ODEProblemLibrary.prob_ode_threebody — Constant

Задача трех тел, описанная Хайером (нежесткая)

Источник: From Hairer Norsett Wanner Solving Ordinary Differential Equations I - Nonstiff Problems Page 129

Обычно при такой настройке решается при и периодичном .

# ODEProblemLibrary.prob_ode_pleiades — Constant

Задача Pleiades (нежесткая)

где

и начальными условиями являются

и с , за исключением

Источник: From Hairer Norsett Wanner Solving Ordinary Differential Equations I - Nonstiff Problems Page 244

Обычно решается от 0 до 3.

# ODEProblemLibrary.prob_ode_vanderpol — Constant

Уравнения Van der Pol

с и ]

Нежесткие параметры

# ODEProblemLibrary.prob_ode_vanderpol_stiff — Constant

Уравнения Van der Pol

с и ]

Жесткие параметры

# ODEProblemLibrary.prob_ode_rober — Constant

Биохимические реакции Робертсона (жесткие)

где , , . Подробные сведения см. в документе:

Hairer Norsett Wanner Solving Ordinary Differential Equations I - Nonstiff Problems Page 129

Обычно решается при ]

# ODEProblemLibrary.prob_ode_rigidbody — Constant

Уравнения твердого тела (нежесткие)

с , и .

Начальное условие — ].

Источник: Solving Differential Equations in R by Karline Soetaert

или Hairer Norsett Wanner Solving Ordinary Differential Equations I - Nonstiff Problems Page 244

Обычно решается от 0 до 20.

# ODEProblemLibrary.prob_ode_hires — Constant

Задача Hires (жесткая)

Она представлена в виде:

где определяется следующим образом:

Справка: demohires.pdf Интерактивный скрипт: Hires.ipynb

# ODEProblemLibrary.prob_ode_orego — Constant

Задача Orego (жесткая)

Она представлен в виде с

где определяется следующим образом:

где , и .

Справка: demoorego.pdf Интерактивный скрипт: Orego.ipynb

# ODEProblemLibrary.prob_ode_pollution — Constant

Задача загрязнения (жесткая)

Эта IVP представляет собой жесткую систему из 20 нелинейных обыкновенных дифференциальных уравнений. Она представлена в виде:

где определяется следующим образом:

с начальным условием

Включен аналитический якобиан.

Справка: pollu.pdf Интерактивный скрипт: Pollution.ipynb

# ODEProblemLibrary.prob_ode_nonlinchem — Constant

Нелинейная система реакций с аналитическим решением

с начальным условием ] во временном интервале

Источник:

Liu, L. C., Tian, B., Xue, Y. S., Wang, M., & Liu, W. J. (2012). Analytic solution for a nonlinear chemistry system of ordinary differential equations. Nonlinear Dynamics, 68(1-2), 17-21.

Реализовано аналитическое решение, позволяющее легко тестировать решатели ОДУ.

# ODEProblemLibrary.prob_ode_brusselator_1d — Constant

1D Брюсселатор

и начальными условиями являются

с краевым условием

Источник: Hairer Norsett Wanner Solving Ordinary Differential Equations II - Stiff and Differential-Algebraic Problems Page 6

# ODEProblemLibrary.prob_ode_brusselator_2d — Constant

2D Брюсселатор

где

и начальными условиями являются

с периодическим краевым условием

Источник: Hairer Norsett Wanner Solving Ordinary Differential Equations II - Stiff and Differential-Algebraic Problems Page 152

# ODEProblemLibrary.prob_ode_filament — Constant

Дискретизация PDE накала

Интерактивный скрипт: Filament.ipynb

В этой задаче использована реальная биологическая модель из документа «Magnetic dipole with a flexible tail as a self-propelling microdevice». Это система PDE, представляющая модель Кирхгофа эластичного стержня, в которой уравнения движения задаются аппроксимацией Роуза со свободными краевыми условиями.

# ODEProblemLibrary.prob_ode_thomas — Constant

Циклически симметричные уравнения аттрактора Томаса