我们研究可变间距的求解器

在这个例子中，我们将分析求解器的选择如何影响使用数学模型计算福柯摆的轨迹。我们将演示求解器的工作DP5和Tsit5（ode45），Rodas4（ode15s），BS3（ode23）和梯形（ode23t). Simulink中采用的类似求解器的名称在括号中给出。

该问题的解决方案将以具有刚性系统特征的微分方程的形式给出。显式积分方法（使用来自相邻点对的数据来查找导数）不适用于此类系统。一个无与伦比的更好的结果通常由隐式方法给出，这些方法根据所需的几个相邻点的数据找到一个解决方案。刚性系统的定义相当模糊，但它们的行为有一个共同特征–它们可能需要一个非常小的积分步骤，没有这个步骤计算是不稳定的（例如，结果走向无穷大）。

刚性通常是由两个分量组成的方程组所固有的，其状态以不同的速率变化：非常快和非常慢。

福柯的钟摆

福柯的钟摆模型是刚性系统的一个例子。钟摆每隔几秒钟振荡一次（一个快速变化的分量），但它也受到地球自转的影响（一个缓慢变化的分量），其频率以天，年等计算。钟摆的运动是由惯性和重力决定的. 并且发生振荡的平面与地球的每日旋转同步旋转。

让我们通过两个微分方程来定义模型，这些微分方程描述了摆锤末端负载在x-y平面（俯视图）中的运动。我们将通过研究系统总能量守恒定律来评估不同求解器的工作。

该模型没有考虑摩擦和空气阻力，因此运动方程看起来非常简洁。:

哪里 -地球绕轴旋转的角速度, -地理纬度, -自由落体加速, 和 -摆锤末端负载质量中心的坐标。

Pkg.add(["Measures", "Serialization"])

using Plots
gr( fmt=:png )

# Параметры модели и начальные условия
Ω = 2*pi/86400;    # Угловая скорость вращения Земли вокруг оси (рад/с)
λ = 59.93/180*pi;  # Широта (рад)
g = 9.8195;        # Ускорение свободного падения (м/с^2)
L = 97;            # Длина маятника (м)

initial_x = L/100; # Начальное значение координаты x (м)
initial_y = 0;     # Начальное значение координаты y (м)
initial_xdot = 0;  # Начальное значение скорости по оси x (м/с)
initial_ydot = 0;  # Начальное значение скорости по оси y (м/с)

TotalEnergy = 1/2*(initial_xdot^2 + initial_ydot^2) + g*(L-sqrt( L^2-initial_x^2-initial_y^2));

实现这些方程的模型如下所示。:

变间距求解器

我们将比较以下积分器（求解器）:

/Name(Julia)/Name(Simulink)/任务类型/精度/何时使用|
| ----------- | ----------- | ----------- | ----------- | ----------- |
/DP5，但更好Tsit5| ode45 /不是刚性/平均/默认。先试试吧|
|BS3| ode23 /非刚性/低/如果在任务中允许大的错误或中等硬系统|
/VCABM，但更好Vern7*|/ ode113 /不死板/低到高/低误差，适用于大容量计算任务|
|QNDF，FBDF，但更好Rodas4**，KenCarp4，TRBDF2，RadauIIA5| ode15s /硬/从低到中/如果 ode45 由于任务的刚性，它的工作速度太慢|
/Rosenbrock23，但更好**Rodas4*|/ ode23s /硬/低/精度要求低，用于解决困难的任务，并具有恒定的[矩阵масс]（https://ru.wikipedia.org/wiki/Матрица_масс )|
|梯形/ ode23t /中等硬/低/对于中等苛刻的任务，如果需要没有数值阻尼的解决方案|
|TRBDF2| ode23tb /硬/低/对解决方案精度要求低，用于解决艰巨的任务|

如果您在解决问题时没有找到合适的求解器，并且您更熟悉Simulink求解器的名称，请注意[本文]中的比较（https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve /)。 Engee主要使用Julia包中接受的求解器的名称。

比较标准

每个积分算法的解的质量可以通过几种方式来评估。首先，如果问题具有解析解，那么您可以将数值积分的结果与理论结果进行比较。您还可以估计算法解决问题所需的时间，并将时间成本用作比较的基础。

不幸的是，福柯钟摆没有确切的解析解，近似的有自己固有的误差。

节约了系统的总能量

我们将通过检查系统中能量守恒定律的遵守情况来评估每个求解器的解决方案的质量。

在没有摩擦的情况下，摆锤的总能量必须保持恒定。由于误差的积累，数值解不会向我们展示理想的能量守恒。

因此，在示例中，我们将在每个积分步骤中计算系统的总能量。理想情况下，总能量不应改变，差值应为零。总能量是势能和动能的总和。这正是块计算的值。 NormalizeEnergy:

哪里 -在初始时刻的总能量。

从计算开始到结束的时刻总能量的变化将使我们能够评估数值解的质量。

modelName = "foucault_pendulum_model"
model = modelName in [m.name for m in engee.get_all_models()] ? engee.open( modelName ) : engee.load( "$(@__DIR__)/$modelName.engee" );

engee.set_param!( modelName, "SolverName"=>"BS3", "RelTol"=>1e-5 ) # Установим решатель, аналогичный ode23
dt_ode23 = @elapsed begin x_ode23 = engee.run( modelName ); end;

engee.set_param!( modelName, "SolverName"=>"Trapezoid", "RelTol"=>1e-5 ) # Установим решатель, аналогичный ode23t
dt_ode23t = @elapsed begin x_ode23t = engee.run( modelName ); end;

让我们通过这个参数比较两个积分器。:

*BS3（ode23）是一个相当粗糙的算法，不适合刚性系统,
*梯形(ode23t）也是一个相当粗糙的算法，但它适用于刚性系统。

using Plots

plot( x_ode23["nE"].time, x_ode23["nE"].value, label="ode23", title="Сравнение качества решений при RelTol = 1e-5", lw=2 )
plot!( x_ode23t["nE"].time, x_ode23t["nE"].value, ylimits=(0,1.1), label="ode23t", leg=:bottomleft, lw=2 )
plot!( xlabel="Время (с)", ylabel="Относительная полная энергия",
       titlefont=11, guidefontsize=10, tickfontsize=9, legendfontsize=10 )

很容易看出，与积分器中的不准确性相关的总能量的"损失"在第一种情况下几乎察觉不到，并且在第二种情况下每天约为20％。

随着摆锤的总能量减小，振荡幅度减小，如下图所示。我们还会看到摆锤的振荡平面随时间变化（展开）。

using Measures

bar_w = 150 # "Толщина" линий фиксирующих плоскость колебания маятника в разные моменты

lastPlane_y = maximum( x_ode23["y"].value[end-bar_w:end] )
lastPlane_x = maximum( x_ode23["x"].value[end-bar_w:end] )
plot( x_ode23["x"].value[1:2:end], x_ode23["y"].value[1:2:end],
      st=:scatter, ms=.8, markerstrokecolor=:blue, markercolor=:blue, 
      title="Положение груза маятника Фуко (вид сверху)\nРешатель: ode23; RelTol = 1e-5; Время выполнения=$(round(dt_ode23,digits=2)) с",
      titlefont=9, guidefontsize=8, tickfontsize=8, leg=:none, xlabel="Координата x (м)", ylabel="Координата y (м)",
      left_margin=13mm, bottom_margin=13mm, top_margin=5mm, aspect_ratio = :equal )
scatter!( [initial_x], [0], mc=:black, ms=8 )                                                # Черный шар
plot!( x_ode23["x"].value[1:bar_w], x_ode23["y"].value[1:bar_w], lc=:black, lw=3 )           # Черная линия
plot!( x_ode23["x"].value[end-bar_w:end], x_ode23["y"].value[end-bar_w:end], lc=:red, lw=3 )

vlen = length(x_ode23["x"].value)
half_vec = ceil(Int32, vlen/2);
plot!( x_ode23["x"].value[(half_vec-bar_w):half_vec], x_ode23["y"].value[(half_vec-bar_w):half_vec], lc=:yellow, lw=3 )

annotate!( initial_x, 0.05, text("Начальное положение", :black, rotation = 90, 7, :left) )
plot!( [lastPlane_x-0.1, lastPlane_x-0.8], [lastPlane_y, lastPlane_y], lw=13, lc=:white ) # Белая линия
p1 = annotate!( lastPlane_x-0.1, lastPlane_y, text("Конечная плоскость", :red, 7, :right) )

plot( x_ode23t["x"].value[1:2:end], x_ode23t["y"].value[1:2:end], st=:scatter, ms=.8, markerstrokecolor=:blue, markercolor=:blue,
          title="Положение груза маятника Фуко (вид сверху)\nРешатель: ode23t; RelTol = 1e-5; Время выполнения=$(round(dt_ode23t,digits=2)) с", 
          titlefont=9, guidefontsize=8, tickfontsize=8, leg=:none, xlabel="Координата x (м)", ylabel="Координата y (м)",
          left_margin=13mm, bottom_margin=13mm, top_margin=5mm, aspect_ratio = :equal )
plot!( x_ode23t["x"].value[1:bar_w], x_ode23t["y"].value[1:bar_w], lc=:black, lw=3 )           # Черная линия
plot!( x_ode23t["x"].value[end-bar_w:end], x_ode23t["y"].value[end-bar_w:end], lc=:red, lw=3 )

vlen = length(x_ode23t["x"].value)
half_vec = ceil(Int32, vlen/2);
plot!( x_ode23t["x"].value[(half_vec-bar_w):half_vec], x_ode23t["y"].value[(half_vec-bar_w):half_vec], lc=:yellow, lw=3 )

p2 = plot!( [initial_x], [0], mc=:black, st=:scatter, ms=8 )

plot( p1, p2, size = (1000,450) )

我们观察到刚性和非刚性系统的解决方案的质量如何取决于求解器的选择。钟摆摆动的平面在白天不会完全转动：转动的角度取决于地理纬度。

让我们建立一个图表，使我们能够详细研究不同求解器在福柯钟摆模拟问题上的速度和工作质量。我们将比较不同的数值积分方法的结果，看看系统的能量在不同的相对误差值下是如何守恒的。

执行此脚本可能需要几分钟，因此为了节省再次运行脚本时的时间，您可以选择从以前保存的文件下载数据。 results.bin.

using Serialization

# Чтобы обновить результаты, удалите файл "results"
if !isfile( "$(@__DIR__)/results.bin" )

    solvers = ["Tsit5" "DP5" "Rosenbrock23" "Trapezoid" "BS3"]
    reltols = 10 .^( collect(range(-7,-3,length=9)) )

    nE_results = zeros( length(solvers), length(reltols) )
    dt_results = zeros( length(solvers), length(reltols) )

    # Изучим все комбинации (решатель; относительная точность)
    for (i,solver) in enumerate(solvers)
        print( solver, ": ")
        for (j,reltol) in enumerate(reltols)
            print( round(reltol, sigdigits=3), ", " )
            # Установим для модели новый решатель и относительную точность
            engee.set_param!( modelName, "SolverName"=>solver, "RelTol"=>reltol )
            # Вычисляем модель и заодно измеряем время
            dt_results[i,j] = @elapsed begin sim_data = engee.run( modelName ); end;
            # Вычисляем соотношение 
            nE_results[i,j] = abs( sim_data["nE"].value[1] - sim_data["nE"].value[end] )
        end
        println()
    end
    
    serialize( "$(@__DIR__)/results.bin", (solvers, reltols, nE_results, dt_results) )
else
    (solvers, reltols, nE_results, dt_results) = deserialize( "$(@__DIR__)/results.bin" );
end

(["Tsit5" "DP5" … "Trapezoid" "BS3"], [1.0000000000000001e-7, 3.162277660168379e-7, 1.0e-6, 3.162277660168379e-6, 1.0e-5, 3.1622776601683795e-5, 0.0001, 0.00031622776601683794, 0.001], [7.929647437576737e-5 7.993461600175333e-5 … 0.0008342080580365785 0.0010495257229566901; 0.006522783592120329 0.006803880028364184 … 0.0432470795415133 0.04324716090734526; … ; 0.00020914013837636247 0.0002237451074058594 … 0.008214673902898095 0.01263204365201287; 0.07957109906963433 0.08391745090197478 … 0.8201637378884246 0.957540689453349], [25.203170908 25.070636927 … 24.543856697 24.405600441; 24.731999073 24.538390433 … 23.893544084 24.32454816; … ; 30.82383809 31.923476917 … 24.9366935 24.482693973; 28.296808969 28.368053586 … 25.189954714 24.714920212])

现在我们可以构建和分析图表了。

markerShapes_ = [ :circle :rect :star5 :diamond :hexagon ]
markerStrokeColors_ = [1 2 3 4 5]

plot(
      plot( reltols, nE_results', xscale=:log10, yscale=:log10,
            label=solvers, ylabel="Макс. относительная\nошибка полной энергии",
            titlefont=11, guidefontsize=8, tickfontsize=9, legendfontsize=7, left_margin=13mm, legend=:bottomleft,
            marker=markerShapes_, markerstrokecolor=markerStrokeColors_, markercolor=RGBA(1,1,1,0) ),
      plot( reltols, dt_results', xscale=:log10,
            label=solvers, ylabel="Время выполнения\nмодели (с)", xlabel="Относительная точность алгоритма численного интегрирования",
            titlefont=11, guidefontsize=8, tickfontsize=9, legendfontsize=7, left_margin=13mm, bottom_margin=13mm,
            marker=markerShapes_, markerstrokecolor=markerStrokeColors_, markercolor=RGBA(1,1,1,0) ),
      size=(800,500),
      layout=(2,1)
    )

我们看到一个区域，预测系统总能量的标准实际上不再改变（在 relative_tolerance < 1e-6).

这是由于特定模型的行为。但是，总的来说，值得考虑的是，求解器相对精度的降低不一定会导致精度的提高。当我们为解设置所需的精度限制并选择其他求解器设置时，我们应该始终考虑在仿真上花费的可接受时间。

我们很可能对在钟摆末端估计负载位置几厘米的精度感到满意。并且不再建议以微米的精度执行计算，因为进行同样准确的真实实验来验证模型将非常困难。

结论

我们研究了求解系统方程的数值方法的选择（求解器的选择）如何影响求解刚性系统的结果和过程。很大程度上取决于求解器：毕竟，算法不仅应该使我们能够实现解的所需精度，而且还应该在合理的时间内完成。

可变步长求解器（以及恒定步长求解器的积分步长）的相对精度应该足够小，以便解不会变得不稳定。

所有的求解器都有其优点和缺点，但有些在解决棘手问题时显然工作得更好。而变步求解器通常比恒步求解器更稳定。

我们研究可变间距的求解器

福柯的钟摆

变间距求解器

比较标准

节约了系统的总能量

结论

示例中使用的块¶