Модели классической физики
Линейное обыкновенное дифференциальное уравнение первого порядка
Радиоактивный распад углерода-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))