Engee 文档
Notebook

学习变步求解器

在本例中，我们将看到求解器的选择如何影响数学模型的傅科摆轨迹计算。我们将演示 DP5 和 Tsit5 (ode45)、Rodas4 (ode15s)、BS3 (ode23) 和 Trapezoid (ode23t) 的性能。括号中给出了 Simulink 中采用的类似求解器的名称。

问题的解将以具有刚性系统特征的微分方程的形式给出。对于此类系统，显式积分法（利用相邻点对的数据求导数）并不十分适用。通常采用_隐式方法_会得到好得多的结果，这种方法从与所需值相邻的几个点的数据中找到 SLAE 的解。刚性系统的定义相当模糊，但其行为的一个共同特征是，它们可能需要一个非常小的积分步长，否则计算将不稳定（例如，结果将变为无穷大）。

刚性通常是由状态变化速度不同的两个部分组成的方程系统的固有特性：非常快和非常慢。

傅科摆

傅科摆模型是刚性系统的一个例子。摆每几秒钟摆动一次（快速变化部分），但它也受到地球自转的影响（慢速变化部分），周期为天、年等。摆的运动是由惯性和地球引力决定的。发生摆动的平面与地球的日自转同步旋转。

我们将通过两个描述摆锤末端砝码在 x-y 平面（俯视图）上运动的微分方程来指定模型。我们将通过研究系统总能量守恒定律来评估不同求解器的性能。

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

$$\ddot x = 2 \Omega sin \lambda \dot y - \frac{g}{L}x$$

$$\ddot y = - 2 \Omega sin \lambda \dot x - \frac{g}{L} y$$

其中$\Omega$ 是地球绕地轴自转的角速度，$\lambda$ 是地理纬度，$g$ 是自由落体的加速度，$x$ 和$y$ 是摆锤末端重物质心的坐标。

In [ ]: 
Pkg.add(["Measures", "Serialization"])
In [ ]: 
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));

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

image.png

可变螺距解算器

我们将对下列积分器（求解器）进行比较：

| 名称（Julia） | 名称（Simulink） | 问题类型 | 精度 | 使用时间 | ----------- | ----------- | ----------- | ----------- | ----------- | | DP5 但更好的是 Tsit5 |ode45 | 不难 | 中等 | 默认。先试试 BS3 |ode23 | 不硬 | 低 | 如果任务容易出错或适用于中等硬度的系统 VCABM 但更好 Vern7 |ode113 | 不难 | 从低到高 | 低误差，适合计算密集型任务 | QNDF 但更好 Vern7 | | 不难 | 从低到高 | 低误差，适合计算密集型任务 | QNDF，FBDF，但更好 Rodas4，KenCarp4，TRBDF2，RadauIIA5 |ode15s | 艰难 | 低到中等 | 如果ode45 由于任务的艰难性而速度太慢 | | VCABM，但更好 Vern7。 | Rosenbrock23，但更好的是Rodas4 |ode23s | 艰难 | 低 | 低 | 当难题的精度要求不高，且[质量矩阵]不变时（https://ru.wikipedia.org/wiki/%D0%9C%D0%B0%D1%82%D1%80%D0%B8%D1%86%D0%B0_%D0%BC%D0%B0%D1%81%D1%81） | | | | Trapezoid |ode23t | 中等刚度 | 低 | 对于中等刚度问题，如果你想要一个没有数值阻尼的解 | TRBDF2 | **TRBDF2 |ode23tb | 刚性 | 低 | 对于精度要求较低的刚性问题

如果您在解决问题时没有找到正确的求解器，而您又更习惯于 Simulink 求解器名称，请参阅 本文 中的比较。Engee 主要使用 Julia 软件包中采用的求解器名称。

比较标准

每种积分算法的求解质量都可以通过几种方式进行评估。首先，如果问题有分析解，我们可以将数值积分的结果与理论结果进行比较。还可以估算这种或那种算法解决问题所需的时间，并将时间成本作为比较的基础。

可惜，傅科摆没有精确的解析解，而近似解析解也有其固有误差。

系统总能量守恒

我们将通过检查系统中的能量守恒定律是否得到满足来评估每个求解器求解的质量。

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

因此，在示例中，我们将计算每个积分步的系统总能量。理想情况下，总能量应该没有任何变化，差值应该为零。总能量是势能和动能之和。块NormalizeEnergy 计算的就是这个值 :

$$E = \frac{v_x^2 + v_y^2}{2} + g(L - \sqrt{L^2 - x^2 - y^2})$$

$$E_{norm} (t) = \frac{E(t)}{E(0)}$$

其中$E(0)$ 是初始时刻的总能量。

通过计算开始到结束时总能量的变化，我们可以评估数值解法的质量。

In [ ]: 
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) - 算法相当粗糙，不适合刚性系统、 Trapezoid (ode23t) - 算法也相当粗糙，但适用于刚性系统。

In [ ]: 
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 )
Out[0]:
No description has been provided for this image

不难看出，在第一种情况下，由于积分器的不准确而造成的总能量 "损失 "几乎微乎其微，而在第二种情况下，每天的 "损失 "约为 20%。

随着摆锤总能量的减少，振荡的振幅也随之减小，我们可以通过绘制下面的图表看到这一点。我们还将看到，摆的振荡平面随时间变化（展开）。

In [ ]: 
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 )
plot!( x_ode23["x"].value[Int32(ceil(end/2)-bar_w):Int32(ceil(end/2))], x_ode23["y"].value[Int32(ceil(end/2)-bar_w):Int32(ceil(end/2))], 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 )
plot!( x_ode23t["x"].value[Int32(ceil(end/2)-bar_w):Int32(ceil(end/2))], x_ode23t["y"].value[Int32(ceil(end/2)-bar_w):Int32(ceil(end/2))], lc=:yellow, lw=3 )
p2 = plot!( [initial_x], [0], mc=:black, st=:scatter, ms=8 )

plot( p1, p2, size = (1000,450) )
Out[0]:
No description has been provided for this image

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

我们将绘制一张图表，以便详细研究不同求解器对傅科摆建模的速度和质量。我们将比较不同数值积分方法的结果，看看在不同的相对误差值下系统的能量是如何守恒的。

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

In [ ]: 
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
Out[0]:
(["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])

现在我们可以绘制和分析图表了。

In [ ]: 
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)
    )
Out[0]:
No description has been provided for this image

我们可以看到系统总能量_非预测标准_几乎停止变化的区域（relative_tolerance < 1e-6 ）。



.

这是由特定模型的行为造成的。但总的来说，值得注意的是，降低求解器的相对精度并不一定会提高精度。当我们对求解的精度设定一个期望的限制并选择其他求解器设置时，我们应该始终考虑到允许的模拟时间。

我们可能会对估计摆锤末端砝码位置几厘米的精度感到满意。而精确到微米的计算不再合理，因为很难进行同样精确的实际实验来验证模型。

结论

我们研究了求解系统方程的数值方法选择（求解器的选择）如何影响刚性系统的求解结果和进度。这在很大程度上取决于求解器：算法不仅要让我们达到所需的求解精度，还要在合理的时间内完成求解。

变步长求解器的相对精度（以及恒步长求解器的积分步长）必须足够小，以防止解变得不稳定。

所有求解器都各有优缺点，但有些求解器在求解难题时效果更好。可变步长求解器通常比恒定步长求解器更稳定。

示例中使用的块