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

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