Численное моделирование физических систем

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

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

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

try
    using DifferentialEquations
catch
    Pkg.add(["DifferentialEquations"])
    using DifferentialEquations
end

Пакет DifferentialEquations.jl – мощная библиотека для решения дифференциальных уравнений, разработанная в рамках экосистемы SciML (Scientific Machine Learning). Эта библиотека предоставляет:

Широкий выбор решателей: от классических методов Рунге–Кутты (Tsit5(), Vern7()) до специализированных алгоритмов для жёстких систем, дифференциально-алгебраических уравнений, стохастических и дифференциальных уравнений с запаздыванием.
Гибкое управление шагом: решатели автоматически подбирают шаг интегрирования для достижения заданной точности, что освобождает пользователя от ручной настройки.
Механизм событий (callbacks): позволяет останавливать интегрирование при выполнении определённых условий (например, касание земли в баллистике) или изменять параметры в процессе расчёта.
Удобное сохранение результатов: параметр saveat позволяет зафиксировать решение в заранее заданные моменты времени, что критически важно для построения плавных анимаций.

Работа с библиотекой строится по единой схеме:

Описание системы ОДУ в виде функции, возвращающей производные (du).
Задание начальных условий и временно́го интервала.
Создание задачи (ODEProblem) и её решение с помощью выбранного решателя.
Обработка и визуализация полученного решения.

Баллистика

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

Движение описывается системой из 4 ОДУ первого порядка для координат и скоростей:

Начальные условия задаются модулем скорости и углом . Интегрирование ведется до момента касания земли ().

Используется пакет DifferentialEquations.jl:

Решатель: Tsit5() (адаптивный метод Рунге–Кутты 5-го порядка).
События (Callbacks): Применяется ContinuousCallback для точного определения момента пересечения . Это останавливает интегрирование (terminate!) ровно в момент падения, исключая ошибки дискретной сетки.
Визуализация: Строится график траектории с отмеченными и подписанными ключевыми точками: старт, максимум и место падения.

g = 9.81
v0 = 20.0
angle = deg2rad(45)
t_flight = 2 * v0 * sin(angle) / g  
tspan = (0.0, t_flight)
u0 = [0.0, 0.0, v0*cos(angle), v0*sin(angle)]

function ballistic!(du, u, p, t)
    du[1] = u[3]
    du[2] = u[4]
    du[3] = 0.0
    du[4] = -g
end

prob = ODEProblem(ballistic!, u0, tspan)
sol = solve(prob, Tsit5())

y_vals = sol[2, :]
idx_max = argmax(y_vals)

x_max = sol[1, idx_max]
y_max = y_vals[idx_max]
x_start, y_start = sol[1, 1], sol[2, 1]
x_end,   y_end   = sol[1, end], sol[2, end]

plot(sol[1, :], sol[2, :], label="Траектория", lw=2,
     xlabel="x (м)", ylabel="y (м)", title="Баллистика",
     legend=:topright)

scatter!([x_start], [y_start], color=:blue,  markersize=8, label="Старт")
scatter!([x_max],   [y_max],   color=:green, markersize=8, label="Макс. высота")
scatter!([x_end],   [y_end],   color=:red,   markersize=8, label="Падение")

annotate!(x_start, y_start - 0.5, text("Старт", :blue, 10))
annotate!(x_max,   y_max   + 0.5, text("Макс. высота", :green, 10))
annotate!(x_end,   y_end   - 0.5, text("Падение", :red, 10))

Гармонические колебания пружинного маятника

Далее моделируются затухающие колебания пружинного маятника. Движение описывается системой из 2 ОДУ первого порядка для смещения и скорости :

где — собственная частота, — коэффициент затухания. Начальные условия: , .

Этот пример расширяет понимание осцилляторов и показывает, как простые изменения в модели (добавление затухания) качественно меняют поведение системы.

ω = 2.0
β = 0.3
A = 1.0
tspan = (0.0, 15.0)

u0 = [A, 0.0]

function harmonic_damped!(du, u, p, t)
    du[1] = u[2]
    du[2] = -ω^2 * u[1] - β * u[2]
end

prob = ODEProblem(harmonic_damped!, u0, tspan)
sol = solve(prob, Tsit5(), saveat=0.02)

@gif for i in 1:2:length(sol.t)
    y = sol[1,i]
    plot([0, 0], [0, -y], color=:black, lw=2, label="Пружина",
         xlims=(-0.5, 0.5), ylims=(-1.5, 1.5),
         aspect_ratio=:equal, title="Затухающие колебания (β=$β)",
         xlabel="", ylabel="Смещение y")
    scatter!([0], [-y], color=:red, markersize=10, label="Груз")
    hline!([0], color=:gray, linestyle=:dash, label="")
end

Движение частицы в постоянном магнитном поле.

Теперь рассмотрим движение заряженной частицы в однородном постоянном магнитном поле, направленном вдоль оси . Под действием силы Лоренца частица движется по окружности (циклотронное движение). Система описывается 4 ОДУ первого порядка для координат и скоростей:

где — циклотронная частота. Начальные условия: частица стартует из начала координат с начальной скоростью вдоль оси .

ω = 1.0
v0 = 2.0
tspan = (0.0, 2π/ω * 4)

u0 = [0.0, 0.0, v0, 0.0]

function lorentz!(du, u, p, t)
    du[1] = u[3]
    du[2] = u[4]
    du[3] = ω * u[4]
    du[4] = -ω * u[3]
end

prob = ODEProblem(lorentz!, u0, tspan)
sol = solve(prob, Tsit5(), saveat=0.02)

@gif for i in 1:length(sol.t)
    plot(sol[1, 1:i], sol[2, 1:i], 
         label="Траектория", lw=2, color=:blue,
         xlabel="x", ylabel="y", 
         title="Движение заряженной частицы в магнитном поле",
         aspect_ratio=:equal, 
         xlims=(-3, 3), ylims=(-5, 1))
    scatter!([sol[1,i]], [sol[2,i]], 
             color=:red, markersize=8, label="Частица")
end

Нелинейный маятник

Перейдём к модели колебания математического маятника без приближения малых углов. В классической задаче часто используют линеаризацию: при малых отклонениях , что приводит к простому гармоническому осциллятору. Однако это приближение работает только при рад (примерно до 10–15°). При больших амплитудах (в данном примере ) нелинейность существенна, и период колебаний зависит от амплитуды.

Система описывается 2 ОДУ первого порядка для угла и угловой скорости :

где — коэффициент затухания. Начальные условия: , .

Этот пример демонстрирует, как численное решение позволяет исследовать системы, не поддающиеся аналитическому решению, и показывает качественное отличие нелинейных колебаний от линейных.

g = 9.81
L = 1.0
β = 0.5
θ0 = deg2rad(120)
ω0 = 0.0
tspan = (0.0, 15.0)

u0 = [θ0, ω0]

function pendulum_damped!(du, u, p, t)
    du[1] = u[2]
    du[2] = -(g/L) * sin(u[1]) - β * u[2]
end

prob = ODEProblem(pendulum_damped!, u0, tspan)
sol = solve(prob, Tsit5(), saveat=0.02)

x_vals = L * sin.(sol[1,:])
y_vals = -L * cos.(sol[1,:])

@gif for i in 1:length(sol.t)
    plot(x_vals[1:i], y_vals[1:i],
         label="Траектория", lw=2, color=:blue,
         xlabel="x", ylabel="y",
         title="Затухающий маятник (β=$β)",
         aspect_ratio=:equal,
         xlims=(-1.5, 1.5), ylims=(-1.5, 1.5))
    plot!([0, x_vals[i]], [0, y_vals[i]],
          color=:black, lw=2, label="Стержень")
    scatter!([x_vals[i]], [y_vals[i]],
             color=:red, markersize=10, label="Груз")
    plot!([0, 0], [-1.5, 0.1],
          color=:gray, linestyle=:dash, label="")
end

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

И в конце рассмотрим двойной математический маятник — два груза, соединённых невесомыми стержнями, колеблющихся в вертикальной плоскости под действием силы тяжести. Это классический пример хаотической системы: движение крайне чувствительно к начальным условиям, а траектория второго груза имеет сложную, непредсказуемую на больших временах форму. В отличие от одинарного маятника, здесь нет аналитического решения даже в консервативном случае.

Система описывается 4 ОДУ первого порядка для углов и угловых скоростей . Уравнения получаются из лагранжиана системы и имеют громоздкий вид из-за нелинейной связи между звеньями:

где . Добавлены коэффициенты затухания в шарнирах. Начальные условия: , .

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

L1 = L2 = 1.0
m1 = m2 = 1.0
β1 = 0.2
β2 = 0.2

θ1_0 = deg2rad(90)
θ2_0 = deg2rad(90)
ω1_0 = 0.0
ω2_0 = 0.0
u0 = [θ1_0, ω1_0, θ2_0, ω2_0]
tspan = (0.0, 15.0)

function double_pendulum_damped!(du, u, p, t)
    θ1, ω1, θ2, ω2 = u
    Δ = θ1 - θ2
    den1 = (m1 + m2) * L1 - m2 * L1 * cos(Δ)^2
    den2 = (L2 / L1) * den1

    du[1] = ω1
    du[2] = ( -g * (m1 + m2) * sin(θ1) - m2 * L1 * ω1^2 * sin(Δ) * cos(Δ)
              - m2 * g * sin(θ2) * cos(Δ) - β1 * ω1 ) / den1
    du[3] = ω2
    du[4] = ( -g * (m1 + m2) * sin(θ1) * cos(Δ)
              - (m1 + m2) * L1 * ω1^2 * sin(Δ)
              + g * (m1 + m2) * sin(θ2) - β2 * ω2 ) / den2
end

prob = ODEProblem(double_pendulum_damped!, u0, tspan)
sol = solve(prob, Tsit5(), saveat=0.02, abstol=1e-8, reltol=1e-8)

x1_vals = L1 * sin.(sol[1,:])
y1_vals = -L1 * cos.(sol[1,:])
x2_vals = x1_vals .+ L2 * sin.(sol[3,:])
y2_vals = y1_vals .- L2 * cos.(sol[3,:])

@gif for i in 1:length(sol.t)
    plot(x2_vals[1:i], y2_vals[1:i],
         label="Траектория второго груза", lw=2, color=:blue,
         xlabel="x", ylabel="y",
         title="Затухающий двойной маятник",
         aspect_ratio=:equal,
         xlims=(-2.5, 2.5), ylims=(-2.5, 2.5))
    plot!([0, x1_vals[i], x2_vals[i]], [0, y1_vals[i], y2_vals[i]],
          color=:black, lw=2, label="Стержни")
    scatter!([x1_vals[i], x2_vals[i]], [y1_vals[i], y2_vals[i]],
             color=:red, markersize=8, label="Грузы")
    plot!([0, 0], [-2.5, 0.1], color=:gray, linestyle=:dash, label="")
end

Вывод

В данном материале мы рассмотрели пять задач классической механики, демонстрирующих возможности численного моделирования в Engee с использованием пакета DifferentialEquations.jl.

Основные результаты:

Баллистика — показала применение событийного механизма (ContinuousCallback) для точного определения момента выполнения условия (касание земли), что критично во многих инженерных задачах.
Гармонический осциллятор с затуханием — продемонстрировал, как диссипативные силы качественно меняют поведение системы, и как с помощью анимации можно наглядно наблюдать экспоненциальное затухание амплитуды.
Движение частицы в магнитном поле — проиллюстрировал циклотронное движение и показал, как постепенная отрисовка траектории в анимации помогает визуализировать эволюцию системы во времени.
Нелинейный маятник — продемонстрировал важность численных методов для систем, не поддающихся аналитическому решению, и показал качественное отличие нелинейных колебаний от линейных при больших амплитудах.
Двойной маятник — представил классический пример хаотической системы, где траектория крайне чувствительна к начальным условиям, и показал необходимость высоких допусков точности (abstol=1e-8, reltol=1e-8) для корректного моделирования таких систем.