Mass-spring-damper model¶
In this example, we will study two approaches to modelling a classical dynamic model involving a single point object with a given mass, connected to a support by a spring and a damper. We will run the models and compare the results when choosing constant and variable integration step solvers.
Pkg.add(["Interpolations", "DelimitedFiles", "PrettyTables"])
using Plots
gr()
Description of the models¶
A dynamic system including a mass, a spring and a damper is very often given to study the nature of oscillatory processes, kinematics and effects of chaotic motion (if there are several objects) or numerical integration methods.
We consider a directed model msd_causal.engee, which performs numerical integration of the equation $m \ddot x + k x + b \dot x = 0$.
It can be clearly seen that if we express from this formula the equation for the acceleration of the load, this is what comes to the input of the integrator $\frac{1}{s^2}$.
$$\ddot x = \frac{-k x - b \dot x}{m}$$
The second model msd_acausal.engee is created from non-directional (acausal) blocks of the physical modelling library Engee. In the upper part there is a sensor that measures the displacement of the mass relative to the fulcrum, as well as a block SolverConfiguration, which indicates to the global model solver that the signals and blocks connected to it will have their own integrator with its own settings.
In the lower part of the model there is a group of elements resembling standard educational illustrations - the mass $m$ is connected to the support by means of a spring with stiffness $k$ and a damper with absorption coefficient $b$ (variable parameters are specified in the properties of these blocks).
.
Due to the graphical nature of the model, these elements can be easily arranged in many different ways, connected in parallel or in series, or repeatedly duplicated.
Setting model parameters¶
Both models use identical variables: weight mass m, spring stiffness k and damper absorption coefficient b.
b,k,m = 1,2,3;
Once the desired variables have been created, you can manipulate the models by changing the values of the variables in the Variables Window. With the values b=1, k=2, m=3 we should observe a rather slow, inertial process with noticeable damping.
Loading and exploring the models¶
With the one-line command, we will either load (load) or open an already loaded model (open). Here we use the syntax of one-line conditional statements x = условие ? если_истинно : если_ложно.
# Однострочная команда для открытия/загрузки модели
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" );
Let's see which integrator parameters are prescribed in both models.
# Отсеим свойства модели, которые нас интересуют
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]
)
We see that both models have the same integrator settings at the top level. The integrator of the physical model is set in the block Solver Configuration, at the top level the model is polled with a constant step 0.01 s.
Execution of models with variable step¶
Let's run the models and store the results into the variables rc and ra for the directional and physical models, respectively.
rc = engee.run( mc, verbose=false );
ra = engee.run( ma, verbose=false );
Comparison of the execution results:
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) )
In general, there will be a different number of points on the plots, so we apply interpolation to compare the results of different models. By doing so, we introduce an additional error, but hopefully not a significant one.
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) )
The difference $\varepsilon$ between the calculation results obtained with a physical model with variable integration step and a directed model with constant step reaches values of the order of $10^{-3}$.
Let us switch to the constant step solver¶
Now we will obtain the solution of the directional model using the method with variable integration step and see how the error changes with respect to the physical model.
# Настроим направленную модель так, чтобы она решалась с переменным шагом интегрирования
engee.set_param!( "msd_causal", "SolverType"=>"variable-step" )
# Запустим модели
rc = engee.run( mc, verbose=false );
ra = engee.run( ma, verbose=false );
To compare plots with different resolution we need interpolation.
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) )
The divergence of the models is also oscillatory. But, this time, the difference between the physical model and the directed model is much smaller. The divergence reaches values of the order of $10^{-12}$.
The Engee has a Data Inspector built in, which will allow you to compare the different results of the runs with each other.
Closing models without saving¶
If we later want to run this example again from the beginning, we must close the models without saving the changes made.
#engee.close( "msd_causal", force=true );
#engee.close( "msd_acausal", force=true );
Conclusion¶
In contrast to directed models, physical models generally cannot be solved using solvers with constant integration step. Difficulties in substituting implicit initial conditions can lead to numerical instability - for example, infinite gradients.
The possibility to investigate models of different nature and to apply different integration methods in an automated way allows us to choose the optimal solution in both respects: the most comprehensible modelling formalism and the most stable solver.