Engee 文档
Notebook

参数不确定的摆式模型

在本演示中，我们将研究一个数学摆模型的运行过程，该模型的参数值是以散点表示的。

与此同时，我们还将检查库DifferentialEquations （用于求解 ODEs）和Plots （绘图）在实现multimethods（多重调度）技术方面有哪些特点，该技术是 Julia 的关键所在。这项技术允许 Julia 中的软件包选择正确的方法来处理最意想不到的数据类型。

准备环境

本演示需要使用Measurements 和DifferentialEquations 库。如果您的系统中没有安装这两个库，请在单独的新单元格或命令行中键入]add Measurements （和]add DifferentialEquations ）并运行命令。

In [ ]: 
Pkg.add(["DifferentialEquations", "Measurements"])
In [ ]: 
using Measurements, DifferentialEquations

首先要做的是设置环境和摆锤本身的参数。

In [ ]: 
g = 9.81 ± 0.01; # Гравитационная постоянная 
L = 1.00 ± 0.02; # Длина подвеса маятника

公差符号± 可以通过命令行上的 LaTeX 命令输入：\pm 。按下 tab 键后，您将得到一个可用于计算的符号，并可在某些库中进行解释，尤其是Measurements 。

Measurements 库的作用是正确处理差值（公差）。例如，从自身减去一个测量值不会使误差差加倍。

In [ ]: 
x = 1 ± 0.01;
x - x
Out[0]:
$0.0 \pm 0.0$

在这方面，加法和减法的作用并不相同。

In [ ]: 
x + x
Out[0]:
$2.0 \pm 0.02$

求解微分方程问题的说明

在本例中，我们实现了一个数学摆的简化模型：

$$\ddot \theta = -\frac{g}{L} \theta$$

其中$g$ 是重力常数，$L$ 是摆的长度，$\theta$ 是摆的偏转角度（以弧度为单位）。我们使用小角度近似$\theta \approx 0$ ，因此假定$sin(\theta) \approx \theta$ 。

让我们设定初始速度和偏转角，以及模拟的时间限制。

In [ ]: 
u₀ = [0 ± 0.0, π/6 ± 0.01]
tspan = (0.0, 6.3)
Out[0]:
(0.0, 6.3)

小写的 "0 "可以作为 LaTeX 命令在命令行中输入：\_0 ，然后按制表符键。转换后，您将得到一个新字符，可将其复制到剪贴板并用于变量名中。

让我们来设置摆锤的方程。

In [ ]: 
function pendulum( u, p, t )
    θ = u[1]
    dθ = u[2]
    du = [dθ, (-g/L) * θ]
    return du
end
Out[0]:
pendulum (generic function with 1 method)

每个向量由两个值组成。状态向量u 由摆的偏转值 θ 和速度 dθ 组成。而导数向量则包括摆的速度 dθ 和作用在砝码上的力。

计算任务的形成和启动

In [ ]: 
prob = ODEProblem( pendulum, u₀, tspan )
Out[0]:
ODEProblem with uType Vector{Measurement{Float64}} and tType Float64. In-place: false
timespan: (0.0, 6.3)
u0: 2-element Vector{Measurement{Float64}}:
   0.0 ± 0.0
 0.524 ± 0.01

我们看到，这个函数的参数u 的类型是Measurement ，参数t 的类型是Float64 ，它们配合得很好。现在让我们运行模拟。

In [ ]: 
sol = solve( prob, Tsit5(), reltol=1e-6 );
In [ ]: 
plot( sol.t, first.(sol.u), label="положение" )
plot!( sol.t, last.(sol.u), label="скорость" )
Out[0]:

我们可以看到，函数plot 自动设置了置信区间，实质上是用参数yerr 增强了图形。

分析解决方案

这种方法的第一个实际应用可能是比较这个问题的数值解和分析解的误差。

In [ ]: 
u = u₀[2] .* cos.( sqrt( g/L ) .* sol.t );
In [ ]: 
plot( sol.t, last.(sol.u), label="Численным методом" )
plot!( sol.t, u, label="Аналитическим", title="Сравнение решений задачи моделирования маятника" )
Out[0]:

我们发现，两种解法的函数值或各点周围的散度几乎没有差别。

备注

本文发布时，Engee文档中尚未包含该库的帮助文本，因此此处提供了开发者页面的链接：https://github.com/JuliaPhysics/Measurements.jl（可通过浏览器内置工具首次翻译）。

使用测量软件包的微分方程的特点

我们已经实现了函数pendulum(u, p, t) ，它使用了return 。DifferentialEquations 软件包也可以实现函数in-place ，其调用方法与pendulum(du, u, p, t) 类似，但由于我们不直接控制变量类型du ，数据类型Measurement 的传播会中断，这种变体的摆式描述将不起作用。因此我们使用return 。

使用测量软件包绘制图表的特点

如图所示，两幅图中第一幅的输出也分为两条命令。函数plot 虽然接受Measurement 类型的数据作为输入，但并没有处理多个此类数据向量的具体方法。我们无法成功调用plot( sol.t, sol.u ) ，出现了类型转换错误。

结论

我们成功地找到了微分方程的数值解，并推导出了一个数学摆模型问题的公差图，该问题的输入值为具有精度公差的数字。