物理系统的数值模拟

该材料是一组演示示例，旨在说明AnyMath及其生态系统在解决古典和现代物理学问题方面的能力。主要目标是展示如何使用简单而富有表现力的代码来模拟机械和电动系统的动力学，观察它们随时间的行为，并以图形和动画的形式可视化结果。

我们将考虑五个关键任务：弹道学，带衰减的谐振子，带电粒子在恒定磁场中的运动，非线性摆，双摆。

每个例子都伴随着常微分方程（ODES）和动画（或静态图）系统的数值解，这有助于可视化物理过程。

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

P**阿克特 DifferentialEquations.jl**是在SciML（科学机器学习）生态系统框架内开发的用于求解微分方程的强大库。这个图书馆提供:

-广泛的求解器：从经典的Runge-Kutta方法（Tsit5(), Vern7()）到刚性系统，微分代数方程，随机和延迟微分方程的专门算法。
-灵活的步骤控制：求解器自动选择积分步骤以达到指定的精度，用户无需手动配置。
-事件机制（回调）：允许您在满足某些条件时停止集成（例如，在弹道学中接触地面）或在计算过程中更改参数。
-方便保存结果：参数 saveat 它允许您在预先设定的时间点捕获决策，这对于构建流畅的动画至关重要。

与图书馆的工作是基于一个单一的方案:

ODE系统的描述以函数返回导数的形式（du).
设置初始条件和时间间隔。
创建任务(ODEProblem)及其使用所选求解器的解。
处理和可视化所得到的解决方案。

弹道学

这个例子模拟了在重力场中以初始速度与地平线成一定角度抛出的材料点的运动。这是学校和大学物理学的经典问题，它清楚地说明了运动学和动力学的基本原理。在我们的例子中，我们考虑一个没有空气阻力的理想化情况，这使我们能够获得抛物线轨迹并使用事件机制来准确确定落到地面的时刻。

运动由4个一阶ODE系统描述，用于坐标和速度:

初始条件由速度模数设定和角度 . 整合进行到接触地面的那一刻（).

包正在使用中 DifferentialEquations.jl:

*求解器: Tsit5() （第5阶的自适应Runge–Kutta方法）。
*事件（回调）：应用 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个一阶ODE系统描述，用于位移和速度 :

哪里 -固有频率, -衰减系数。初始条件: , .

此示例扩展了对振荡器的理解，并显示了模型中的简单更改（添加衰减）如何定性地改变系统的行为。

ω = 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个一阶ODE描述:

哪里 -回旋加速器频率。初始条件：粒子以初始速度从原点开始沿轴线 .

ω = 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

非线性钟摆

让我们继续讨论数学摆的振荡模型，而不近似小角度。在经典问题中，经常使用线性化：对于小偏差，从而产生简谐振动器。但是，此近似值仅适用于以下情况 rad(高达约10-15°)。在大振幅（在本例中非线性是显着的，振荡周期取决于振幅。

该系统由角度的2个一阶ODE描述和角速度 :

哪里 -衰减系数。初始条件: , .

这个例子演示了数值解如何使研究无法解析求解的系统成为可能，并显示了非线性和线性振荡之间的定性差异。

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个角度的一阶ODE描述和角速度 . 这些方程是从系统的拉格朗日导出的，由于链接之间的非线性关系，具有繁琐的外观。:

哪里 . 增加衰减系数在铰链里。初始条件: , .

这个例子清楚地表明了数值方法如何使研究根本不可能解析解的系统成为可能，并显示了混沌特征相空间轨迹的混合。

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

结论

在本文中，我们回顾了经典力学的五个问题，这些问题展示了使用该软件包在AnyMath中进行数值建模的可能性 DifferentialEquations.jl.

主要结果:

弹道学-显示使用事件机制（ContinuousCallback）以准确地确定何时满足条件（接触地面），这在许多工程任务中至关重要。
带衰减的谐振子-演示了耗散力如何定性地改变系统的行为，以及如何在动画的帮助下，可以直观地观察振幅的指数衰减。
粒子在磁场中的运动--展示了回旋加速器运动，并展示了动画中轨迹的逐步绘制如何帮助可视化系统随时间的演变。
非线性摆锤-证明了数值方法对于无法解析求解的系统的重要性，并显示了大振幅下非线性和线性振荡之间的定性差异。
双摆-提出了一个混沌系统的经典例子，其中轨迹对初始条件非常敏感，并表明需要高精度公差（abstol=1e-8, reltol=1e-8），以便对这类系统进行正确的建模。