Engee 文档
Notebook

质量-弹簧-阻尼模型

在本例中，我们将研究一种经典动态模型的两种建模方法，该模型涉及一个具有给定质量的单点物体，该物体通过弹簧和阻尼器与支撑物相连。我们将运行模型，并比较选择恒定和可变积分步骤求解器时的结果。

In [ ]: 
Pkg.add(["Interpolations", "DelimitedFiles", "PrettyTables"])
In [ ]: 
using Plots
gr()
Out[0]:
Plots.GRBackend()

模型描述

为了研究振荡过程的性质、运动学和混沌运动的影响（如果有多个物体）或数值积分方法，通常会给出一个包括质量、弹簧和阻尼器的动态系统。

image.png

我们考虑一个定向模型msd_causal.engee ，它对方程$m \ddot x + k x + b \dot x = 0$ 进行数值积分。

image.png

可以清楚地看到，如果我们从这个公式中表达出负载的加速度方程，这就是积分器$\frac{1}{s^2}$ 的输入。

$$\ddot x = \frac{-k x - b \dot x}{m}$$

第二个模型msd_acausal.engee 是由物理建模库Engee的非定向（无因果关系）块创建的。上部有一个传感器，用于测量质量相对于支点的位移，还有一个块SolverConfiguration ，它向全局模型求解器表明，与之相连的信号和块将有自己的积分器，并有自己的设置。

image.png

在模型的下部有一组类似于标准教学插图的元素--质量$m$ 通过刚度为$k$ 的弹簧和吸收系数为$b$ 的阻尼器与支架相连（可变参数在这些块的属性中指定）。

.

由于模型的图形性质，这些元素可以很容易地以多种不同的方式排列、并联或串联，或重复复制。

设置模型参数

两个模型使用相同的变量：重量质量m 、弹簧刚度k 和减震器吸收系数b 。

In [ ]: 
b,k,m = 1,2,3;

一旦创建了所需的变量，您就可以通过更改变量窗口中的变量值来操作模型。使用b=1,k=2,m=3 的值，我们可以观察到一个相当缓慢的惯性过程，并带有明显的阻尼。

加载和探索模型

通过单行命令，我们可以加载 (load) 或打开已经加载的模型 (open) 。这里我们使用单行条件语句的语法x = условие ? если_истинно : если_ложно 。

In [ ]: 
# Однострочная команда для открытия/загрузки модели
mc = "msd_causal" in [m.name for m in engee.get_all_models()] ? engee.open( "msd_causal" ) : engee.load( "$(@__DIR__)/msd_causal.engee" );
ma = "msd_acausal" in [m.name for m in engee.get_all_models()] ? engee.open( "msd_acausal" ) : engee.load( "$(@__DIR__)/msd_acausal.engee" );

让我们来看看这两个模型都规定了哪些积分器参数。

In [ ]: 
# Отсеим свойства модели, которые нас интересуют
mc_s = split( string( engee.get_param( "msd_causal" ) ), "#Параметры интегратора:\n")[2];
ma_s = split( string( engee.get_param( "msd_acausal" ) ), "#Параметры интегратора:\n")[2];
# Приведем их в соответствие стандартному форматированию CSV файлов
mc_s = replace( strip(mc_s, [' ', '\n', ')'] ), "=>"=>";" )
ma_s = replace( strip(ma_s, [' ', '\n', ')'] ), "=>"=>";" )

# Таблица для сравнения параметров моделей
using DelimitedFiles, DataFrames, PrettyTables
pretty_table(
    coalesce.(outerjoin(
        DataFrame( readdlm( IOBuffer(mc_s), ';' ), ["Свойство", "Направленная модель"] ),
        DataFrame( readdlm( IOBuffer(ma_s), ';' ), ["Свойство", "Физическая модель"] ),
        on="Свойство" ), "-"),
    backend=:html, nosubheader=true, alignment=[:r, :c, :c]
    )
Свойство Направленная модель Физическая модель
:StartTime 0.0 0.0
:StopTime 10.0 10.0
:SolverType fixed-step fixed-step
:SolverName Euler Euler
:FixedStep 0.01 0.01

我们可以看到，两个模型在顶层都有相同的积分器设置。物理模型的积分器设置在程序块Solver Configuration 中，在顶层，模型以0.01 s 的恒定步长进行轮询。

执行步长可变的模型

让我们运行模型，并将结果分别存储到定向模型和物理模型的变量rc 和ra 中。

In [ ]: 
rc = engee.run( mc, verbose=false );
ra = engee.run( ma, verbose=false );

比较执行结果：

In [ ]: 
plot( rc["p"].time, rc["p"].value, title="Положение" )
p_pos = plot!( ra["p"].time, ra["p"].value )

plot( rc["v"].time, rc["v"].value, title="Скорость" )
p_vel = plot!( ra["v"].time, ra["v"].value )

plot( p_pos, p_vel, size=(600,200) )
Out[0]:

一般来说，图上的点数会有所不同，因此我们采用插值法来比较不同模型的结果。这样做会带来额外的误差，但希望不会很大。

In [ ]: 
using Interpolations
rc_p_interp_linear = linear_interpolation(rc["p"].time, rc["p"].value)
rc_v_interp_linear = linear_interpolation(rc["v"].time, rc["v"].value)

using LaTeXStrings
d_pos = plot( ra["p"].value - rc_p_interp_linear.(ra["p"].time), title=L"\varepsilon_{положение}", legend=false)
d_vel = plot( ra["v"].value - rc_v_interp_linear.(ra["v"].time), title=L"\varepsilon_{скорость}", legend=false )

plot( d_pos, d_vel, size=(600,200) )
Out[0]:

使用积分步长可变的物理模型和积分步长恒定的定向模型得出的计算结果之间的差值$\varepsilon$ 达到$10^{-3}$ 的数量级。

让我们切换到恒定步长求解器

现在，我们将使用积分步长可变的方法求得方向模型的解，看看误差与物理模型的关系如何变化。

In [ ]: 
# Настроим направленную модель так, чтобы она решалась с переменным шагом интегрирования
engee.set_param!( "msd_causal", "SolverType"=>"variable-step" )

# Запустим модели
rc = engee.run( mc, verbose=false );
ra = engee.run( ma, verbose=false );

要比较不同分辨率的曲线图，我们需要插值法。

In [ ]: 
rc_p_interp_linear = linear_interpolation(rc["p"].time, rc["p"].value)
rc_v_interp_linear = linear_interpolation(rc["v"].time, rc["v"].value)

e_pos = plot( ra["p"].value - rc_p_interp_linear.(ra["p"].time), title=L"\varepsilon_{положение}", legend=false )
e_vel = plot( ra["v"].value - rc_v_interp_linear.(ra["v"].time), title=L"\varepsilon_{скорость}", legend=false )

plot( e_pos, e_vel, size=(600,200) )
Out[0]:

模型的分歧也是振荡的。不过，这次物理模型和定向模型之间的差异要小得多。发散值达到了$10^{-12}$ 的数量级。

Engee**内置了一个 "数据检查器"，可以比较不同的运行结果。

image.png

关闭模型而不保存

如果我们以后想从头开始重新运行这个示例，就必须在不保存更改的情况下关闭模型。

In [ ]: 
#engee.close( "msd_causal", force=true );
#engee.close( "msd_acausal", force=true );

结论

与定向模型相比，物理模型通常无法使用恒定积分步长的求解器求解。替换隐含初始条件的困难会导致数值不稳定--例如，无限梯度。

研究不同性质的模型和自动应用不同积分方法的可能性，使我们能够从两方面选择最佳解决方案：最易理解的建模形式和最稳定的求解器。

示例中使用的块