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

Если вы не решаетесь окунуться в мир 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...)