Модели классической физики

Если вы не решаетесь окунуться в мир DiffEq, вот несколько специально придуманных небольших задач с дифференциальными уравнениями, которые помогут вам освоить азы.

Линейное обыкновенное дифференциальное уравнение первого порядка

Радиоактивный распад углерода-14

Задача радиоактивного распада — это задача с линейным ОДУ первого порядка, включающим экспоненциальную зависимость с отрицательным коэффициентом, которая представляет период полураспада. Если изменить коэффициент на положительный, получится уравнение роста населения.

using OrdinaryDiffEq, Plots
gr()

#Период полураспада углерода-14 составляет 5730 лет.
C₁ = 5.730

#Настройка
u₀ = 1.0
tspan = (0.0, 1.0)

#Определение задачи
radioactivedecay(u, p, t) = -C₁ * u

#Передача в решатель
prob = ODEProblem(radioactivedecay, u₀, tspan)
sol = solve(prob, Tsit5())

#График
plot(sol, linewidth = 2, title = "Carbon-14 half-life",
     xaxis = "Time in thousands of years", yaxis = "Percentage left",
     label = "Numerical Solution")
plot!(sol.t, t -> exp(-C₁ * t), lw = 3, ls = :dash, label = "Analytical Solution")

Линейное обыкновенное дифференциальное уравнение второго порядка

Простой гармонический осциллятор

Еще один классический пример — гармонический осциллятор, описываемый следующим уравнением:

с известным аналитическим решением:

где

Константы определяются начальными условиями так, что — это начальное положение, а — начальная скорость.

Вместо преобразования этого уравнения в систему ОДУ для решения с помощью ODEProblem мы можем использовать SecondOrderODEProblem следующим образом.

# Задача простого гармонического осциллятора
using OrdinaryDiffEq, Plots

#Параметры
ω = 1

#Начальные условия
x₀ = [0.0]
dx₀ = [π / 2]
tspan = (0.0, 2π)

ϕ = atan((dx₀[1] / ω) / x₀[1])
A = √(x₀[1]^2 + dx₀[1]^2)

#Определение задачи
function harmonicoscillator(ddu, du, u, ω, t)
    ddu .= -ω^2 * u
end

#Передача в решатели
prob = SecondOrderODEProblem(harmonicoscillator, dx₀, x₀, tspan, ω)
sol = solve(prob, DPRKN6())

#График
plot(sol, vars = [2, 1], linewidth = 2, title = "Simple Harmonic Oscillator",
     xaxis = "Time", yaxis = "Elongation", label = ["x" "dx"])
plot!(t -> A * cos(ω * t - ϕ), lw = 3, ls = :dash, label = "Analytical Solution x")
plot!(t -> -A * ω * sin(ω * t - ϕ), lw = 3, ls = :dash, label = "Analytical Solution dx")

Обратите внимание, что переменные (и начальные условия) следуют в порядке dx, x. Таким образом, если мы хотим, чтобы первым рядом была переменная x, необходимо поменять порядок с помощью выражения vars=[2,1].

Нелинейное обыкновенное дифференциальное уравнение второго порядка

Математический маятник

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

Но у нас есть численные решатели ОДУ! Почему бы не решить задачу реального маятника?

Обратите внимание, что теперь у нас ОДУ второго порядка. Чтобы можно было применить тот же метод, что и ранее, необходимо преобразовать это уравнение в систему ОДУ первого порядка с использованием нотации .

# Задача простого маятника
using OrdinaryDiffEq, Plots

#Константы
const g = 9.81
L = 1.0

#Начальные условия
u₀ = [0, π / 2]
tspan = (0.0, 6.3)

#Определение задачи
function simplependulum(du, u, p, t)
    θ = u[1]
    dθ = u[2]
    du[1] = dθ
    du[2] = -(g / L) * sin(θ)
end

#Передача в решатели
prob = ODEProblem(simplependulum, u₀, tspan)
sol = solve(prob, Tsit5())

#График
plot(sol, linewidth = 2, title = "Simple Pendulum Problem", xaxis = "Time",
     yaxis = "Height", label = ["\\theta" "d\\theta"])

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

p = plot(sol, vars = (1, 2), xlims = (-9, 9), title = "Phase Space Plot",
         xaxis = "Velocity", yaxis = "Position", leg = false)
function phase_plot(prob, u0, p, tspan = 2pi)
    _prob = ODEProblem(prob.f, u0, (0.0, tspan))
    sol = solve(_prob, Vern9()) # Используем решатель Vern9 для повышения точности
    plot!(p, sol, vars = (1, 2), xlims = nothing, ylims = nothing)
end
for i in (-4pi):(pi / 2):(4π)
    for j in (-4pi):(pi / 2):(4π)
        phase_plot(prob, [j, i], p)
    end
end
plot(p, xlims = (-9, 9))

Двойной маятник

Более сложным примером может служить двойной маятник. Уравнения, описывающие его движение, выглядят так (они взяты из этого вопроса на Stack Overflow):

#Задача двойного маятника
using OrdinaryDiffEq, Plots

#Константы и настройка
const m₁, m₂, L₁, L₂ = 1, 2, 1, 2
initial = [0, π / 3, 0, 3pi / 5]
tspan = (0.0, 50.0)

#Вспомогательная функция для преобразования полярных координат в декартовы
function polar2cart(sol; dt = 0.02, l1 = L₁, l2 = L₂, vars = (2, 4))
    u = sol.t[1]:dt:sol.t[end]

    p1 = l1 * map(x -> x[vars[1]], sol.(u))
    p2 = l2 * map(y -> y[vars[2]], sol.(u))

    x1 = l1 * sin.(p1)
    y1 = l1 * -cos.(p1)
    (u, (x1 + l2 * sin.(p2),
         y1 - l2 * cos.(p2)))
end

#Определение задачи
function double_pendulum(xdot, x, p, t)
    xdot[1] = x[2]
    xdot[2] = -((g * (2 * m₁ + m₂) * sin(x[1]) +
                 m₂ * (g * sin(x[1] - 2 * x[3]) +
                  2 * (L₂ * x[4]^2 + L₁ * x[2]^2 * cos(x[1] - x[3])) * sin(x[1] - x[3]))) /
                (2 * L₁ * (m₁ + m₂ - m₂ * cos(x[1] - x[3])^2)))
    xdot[3] = x[4]
    xdot[4] = (((m₁ + m₂) * (L₁ * x[2]^2 + g * cos(x[1])) +
                L₂ * m₂ * x[4]^2 * cos(x[1] - x[3])) * sin(x[1] - x[3])) /
              (L₂ * (m₁ + m₂ - m₂ * cos(x[1] - x[3])^2))
end

#Передача в решатели
double_pendulum_problem = ODEProblem(double_pendulum, initial, tspan)
sol = solve(double_pendulum_problem, Vern7(), abstol = 1e-10, dt = 0.05);

retcode: Success
Interpolation: specialized 7th order lazy interpolation
t: 303-element Vector{Float64}:
  0.0
  0.05
  0.12217751559731858
  0.21737046504966892
  0.32557871102475433
  0.4538001831296513
  0.6093195612067883
  0.7728295888824992
  0.9531773177034342
  1.1789364907473971
  ⋮
 48.69306260963273
 48.91021204184887
 49.07673144456275
 49.26950145313383
 49.413714336117856
 49.58639229759933
 49.735034200347606
 49.92709346084078
 50.0
u: 303-element Vector{Vector{Float64}}:
 [0.0, 1.0471975511965976, 0.0, 1.8849555921538759]
 [0.05276815671595484, 1.071438957072351, 0.09384176137401032, 1.8607344940613866]
 [0.13346404949254442, 1.1746261361732084, 0.22461544162151867, 1.7498682935067824]
 [0.25340014853291515, 1.3397062143942808, 0.38068162335792477, 1.5194590635443397]
 [0.40362067267647705, 1.4025944942000763, 0.5297354584756716, 1.240443699864503]
 [0.5722968301269858, 1.171450199082293, 0.6713718147242973, 0.9861902079902164]
 [0.7085423342556784, 0.5399035287557558, 0.8079790490828, 0.7794821214727354]
 [0.7355031463861064, -0.1908205822316285, 0.9166518878844778, 0.5304216598453477]
 [0.6478012338021546, -0.7138455590058932, 0.9733991108628701, 0.057112350467565534]
 [0.4700665413819645, -0.734304300371504, 0.888698401649624, -0.8582767128313853]
 ⋮
 [-0.6790142112405692, 0.1666713249565996, -0.6449607059987224, 1.445707398332814]
 [-0.453711223749631, 1.907772695621923, -0.3662711842230573, 1.0713657099932998]
 [-0.08088755404381312, 2.283085230936697, -0.19537310870031996, 1.1164019590612861]
 [0.25211550091772744, 1.0205834408401855, 0.08848962912628953, 1.8493662815063088]
 [0.3327108424404666, 0.2519694839950239, 0.37935712531815796, 2.073023391326819]
 [0.39200740461082473, 0.5635940086140278, 0.6950083099731763, 1.4813704445243048]
 [0.5030292288143082, 0.8775432091086321, 0.8641350711861532, 0.8027494923114267]
 [0.6687647784296428, 0.7437942081020751, 0.9462230447555693, 0.0961700207091559]
 [0.7159266100692314, 0.5372537705329169, 0.9459035223663264, -0.09804066261876586]

#Получение координат в декартовой системе
ts, ps = polar2cart(sol, l1 = L₁, l2 = L₂, dt = 0.01)
plot(ps...)

Сечение Пуанкаре

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

В данном случае сечение Пуанкаре задается коллекцией , где и .

#Константы и настройка
using OrdinaryDiffEq
initial2 = [0.01, 0.005, 0.01, 0.01]
tspan2 = (0.0, 500.0)

#Определение задачи
function double_pendulum_hamiltonian(udot, u, p, t)
    α = u[1]
    lα = u[2]
    β = u[3]
    lβ = u[4]
    udot .= [2(lα - (1 + cos(β))lβ) / (3 - cos(2β)),
        -2sin(α) - sin(α + β),
        2(-(1 + cos(β))lα + (3 + 2cos(β))lβ) / (3 - cos(2β)),
        -sin(α + β) - 2sin(β) * (((lα - lβ)lβ) / (3 - cos(2β))) +
        2sin(2β) * ((lα^2 - 2(1 + cos(β))lα * lβ + (3 + 2cos(β))lβ^2) / (3 - cos(2β))^2)]
end

# Создание ContiunousCallback
condition(u, t, integrator) = u[1]
affect!(integrator) = nothing
cb = ContinuousCallback(condition, affect!, nothing,
                        save_positions = (true, false))

# Формирование задачи
poincare = ODEProblem(double_pendulum_hamiltonian, initial2, tspan2)
sol2 = solve(poincare, Vern9(), save_everystep = false, save_start = false,
             save_end = false, callback = cb, abstol = 1e-16, reltol = 1e-16)

function poincare_map(prob, u₀, p; callback = cb)
    _prob = ODEProblem(prob.f, u₀, prob.tspan)
    sol = solve(_prob, Vern9(), save_everystep = false, save_start = false,
                save_end = false, callback = cb, abstol = 1e-16, reltol = 1e-16)
    scatter!(p, sol, vars = (3, 4), markersize = 3, msw = 0)
end

poincare_map (generic function with 1 method)

lβrange = -0.02:0.0025:0.02
p = scatter(sol2, vars = (3, 4), leg = false, markersize = 3, msw = 0)
for lβ in lβrange
    poincare_map(poincare, [0.01, 0.01, 0.01, lβ], p)
end
plot(p, xlabel = "\\beta", ylabel = "l_\\beta", ylims = (0, 0.03))

Система Эно — Эйлеса

Потенциал Хенона-Хейлеса возникает при нелинейном движении звезды вокруг центра галактики в одной плоскости.

где

В данном случае мы выбираем , поэтому

Тогда полную энергию системы можно выразить так:

Полная энергия должна сохраняться при эволюции системы.

using OrdinaryDiffEq, Plots

#Настройка
initial = [0.0, 0.1, 0.5, 0]
tspan = (0, 100.0)

#Помните, что V — это потенциал системы, а T — полная кинетическая энергия, поэтому E будет
#полной энергией системы.
V(x, y) = 1 // 2 * (x^2 + y^2 + 2x^2 * y - 2 // 3 * y^3)
E(x, y, dx, dy) = V(x, y) + 1 // 2 * (dx^2 + dy^2);

#Определение функции
function Hénon_Heiles(du, u, p, t)
    x = u[1]
    y = u[2]
    dx = u[3]
    dy = u[4]
    du[1] = dx
    du[2] = dy
    du[3] = -x - 2x * y
    du[4] = y^2 - y - x^2
end

#Передача в решатели
prob = ODEProblem(Hénon_Heiles, initial, tspan)
sol = solve(prob, Vern9(), abstol = 1e-16, reltol = 1e-16);

retcode: Success
Interpolation: specialized 9th order lazy interpolation
t: 1846-element Vector{Float64}:
   0.0
   0.011644644151844885
   0.02264023706311373
   0.03351618364134173
   0.04503735852340375
   0.05892780544806496
   0.07129866971822844
   0.08699618005665914
   0.10247212827681021
   0.1196036739006969
   ⋮
  99.59389021446746
  99.65049695936257
  99.70206554871487
  99.75689255207548
  99.81105766571665
  99.86910986733137
  99.92437958961601
  99.97670729241455
 100.0
u: 1846-element Vector{Vector{Float64}}:
 [0.0, 0.1, 0.5, 0.0]
 [0.005822164178948816, 0.09999389777392866, 0.49995932143722593, -0.0010481306031208865]
 [0.011318958086603432, 0.09997692919997468, 0.49984623673766104, -0.0020384490164109133]
 [0.016754327180830298, 0.09994942744649325, 0.4996630516825552, -0.0030191412361624836]
 [0.022509545576229435, 0.09990865029656407, 0.4993916606128625, -0.0040598744112650765]
 [0.02944344598665176, 0.09984352324358772, 0.4989587515150901, -0.00531808262509163]
 [0.03561310537856846, 0.09977078214677518, 0.4984760158664954, -0.006442691986661458]
 [0.0434322895116093, 0.09965840510082306, 0.4977318621044846, -0.007876507439329966]
 [0.05112855611141907, 0.09952551398935892, 0.49685438693972966, -0.009298984640121525]
 [0.0596309453399941, 0.0993526352773756, 0.49571692895968983, -0.010885809100767009]
 ⋮
 [0.11965128358447895, 0.4263383139792122, 0.31908837067746176, -0.0298319578621376]
 [0.1373412483027802, 0.42423261985963345, 0.3056204514639669, -0.04460527014653252]
 [0.15275165038003788, 0.4215807427072929, 0.2918078825857258, -0.05827819587252601]
 [0.16831341097865576, 0.4179817391236895, 0.27560879868039123, -0.07304800625486678]
 [0.18277600548206052, 0.41362452979806796, 0.25818464707238603, -0.08787741630229802]
 [0.19718748124579646, 0.40805552146512836, 0.23808236201505298, -0.10402807176258946]
 [0.20978857949226154, 0.4018755386827257, 0.21771627624490925, -0.11963832667413636]
 [0.22065520422382953, 0.3952242213625343, 0.19746807322814197, -0.13460853683713261]
 [0.22514698918544104, 0.3920106535852837, 0.1881879739924063, -0.14132567096174772]

# График орбиты
plot(sol, vars = (1, 2), title = "The orbit of the Hénon-Heiles system", xaxis = "x",
     yaxis = "y", leg = false)

#Дополнительная проверка правильности — что, на ваш взгляд, возвращает этот код и почему?
@show sol.retcode

#График
plot(sol, vars = (1, 3), title = "Phase space for the Hénon-Heiles system",
     xaxis = "Position", yaxis = "Velocity")
plot!(sol, vars = (2, 4), leg = false)

#Мы сопоставляем полные энергии на временных интервалах решения (в данном случае sol.u) с новым вектором
#и передаем его в построитель графика, что немного удобнее
energy = map(x -> E(x...), sol.u)

#Здесь мы используем @show, чтобы легко выявлять ошибочное поведение системы исходя из слишком большой потери энергии.
@show ΔE = energy[1] - energy[end]

#График
plot(sol.t, energy .- energy[1], title = "Change in Energy over Time",
     xaxis = "Time in iterations", yaxis = "Change in Energy")

Симплектическая интеграция

Чтобы избежать дрейфа энергии, можно вместо этого воспользоваться симплектическим интегратором. Можно напрямую определить и решить SecondOrderODEProblem:

function HH_acceleration!(dv, v, u, p, t)
    x, y = u
    dx, dy = dv
    dv[1] = -x - 2x * y
    dv[2] = y^2 - y - x^2
end
initial_positions = [0.0, 0.1]
initial_velocities = [0.5, 0.0]
prob = SecondOrderODEProblem(HH_acceleration!, initial_velocities, initial_positions, tspan)
sol2 = solve(prob, KahanLi8(), dt = 1 / 10);

retcode: Success
Interpolation: 3rd order Hermite
t: 1001-element Vector{Float64}:
   0.0
   0.1
   0.2
   0.30000000000000004
   0.4
   0.5
   0.6
   0.7
   0.7999999999999999
   0.8999999999999999
   ⋮
  99.19999999999864
  99.29999999999863
  99.39999999999863
  99.49999999999862
  99.59999999999862
  99.69999999999861
  99.7999999999986
  99.8999999999986
 100.0
u: 1001-element Vector{ArrayPartition{Float64, Tuple{Vector{Float64}, Vector{Float64}}}}:
 ([0.5, 0.0], [0.0, 0.1])
 ([0.497004124813899, -0.009071101031878595], [0.049900082497367014, 0.09954822053953173])
 ([0.488065986503409, -0.01856325999777532], [0.09920263962777168, 0.0981717138550656])
 ([0.4733339415930383, -0.028870094978230895], [0.1473200401467506, 0.09580835287880228])
 ([0.45305474121099776, -0.040331734000937175], [0.1936844146808317, 0.09235911550101675])
 ([0.4275722362076315, -0.0532114683862374], [0.2377574398051348, 0.08769462857458447])
 ([0.39732459499168415, -0.06767627591286327], [0.279039924450448, 0.08166386202131776])
 ([0.36283926518715554, -0.08378216859669584], [0.3170810117622551, 0.07410453631898412])
 ([0.3247249463611115, -0.10146505675930295], [0.35148673325356056, 0.06485472395637552])
 ([0.2836600192614762, -0.12053752272622162], [0.3819275865822665, 0.05376507097805092])
 ⋮
 ([0.35754893663581877, 0.06817150704505637], [-0.01689940464435974, 0.41858730953689555])
 ([0.3573608346071073, 0.043775799566172786], [0.018901094264629065, 0.42418546702794857])
 ([0.35056546788922516, 0.01917894635581023], [0.054352325047751296, 0.4273357603171658])
 ([0.3372619546371046, -0.00582469608876845], [0.08879705585509527, 0.4280076678341491])
 ([0.317723101892544, -0.031415069954177], [0.12159668542307678, 0.42615121224761354])
 ([0.2923883191107272, -0.057726485993246215], [0.15214830714764738, 0.42170054930538214])
 ([0.2618505981512214, -0.0848309469334548], [0.17990076750513168, 0.4145793965132324])
 ([0.22683770400643935, -0.11272570532234003], [0.2043691308538718, 0.40470792450133325])
 ([0.1881879739856503, -0.1413256709667938], [0.22514698919068493, 0.39201065357989734])

Обратите внимание, что результаты будут теми же:

# Строим график орбиты
plot(sol2, vars = (3, 4), title = "The orbit of the Hénon-Heiles system", xaxis = "x",
     yaxis = "y", leg = false)

plot(sol2, vars = (3, 1), title = "Phase space for the Hénon-Heiles system",
     xaxis = "Position", yaxis = "Velocity")
plot!(sol2, vars = (4, 2), leg = false)

Однако в данном случае изменение энергии практически равно нулю:

energy = map(x -> E(x[3], x[4], x[1], x[2]), sol2.u)
#Здесь мы используем @show, чтобы легко выявлять ошибочное поведение системы исходя из слишком большой потери энергии.
@show ΔE = energy[1] - energy[end]

#График
plot(sol2.t, energy .- energy[1], title = "Change in Energy over Time",
     xaxis = "Time in iterations", yaxis = "Change in Energy")

И давайте еще попробуем решить эту задачу с помощью решателя Рунге-Кутты-Нюстрёма. Обратите внимание, что решатель Рунге-Кутты-Нюстрёма не симплектический.

sol3 = solve(prob, DPRKN6());
energy = map(x -> E(x[3], x[4], x[1], x[2]), sol3.u)
@show ΔE = energy[1] - energy[end]
gr()
plot(sol3.t, energy .- energy[1], title = "Change in Energy over Time",
     xaxis = "Time in iterations", yaxis = "Change in Energy")

Обратите внимание, что мы используем решатель DPRKN6 с погрешностью reltol=1e-3 (значение по умолчанию), однако изменение энергии меньше, чем при использовании Vern9 с погрешностями abstol=1e-16, reltol=1e-16. Таким образом, специализированные решатели очень эффективны для решения соответствующих задач.