Learning variable step solvers¶
In this example we will see how the choice of solver affects the computation of the Foucault pendulum trajectory from a mathematical model. We will demonstrate the performance of DP5 and Tsit5 (ode45), Rodas4 (ode15s), BS3 (ode23), and Trapezoid (ode23t). The names of similar solvers adopted in Simulink are given in brackets.
The solution of the problem will be given in the form of differential equations having the features of a rigid system. For such systems, explicit integration methods (which use data from pairs of neighbouring points to find the derivative) are not well suited. An incomparably better result is usually obtained by implicit methods, which find the solution to the SLAE from the data of a few neighbouring points from the desired one. The definition of rigid systems is rather vague, but one common feature of their behaviour is that they can require a very small integration step, without which the computation is unstable (e.g., the result goes to infinity).
Rigidity is often inherent in systems of equations consisting of two components whose state changes at different rates: very fast and very slow.
Foucault's pendulum¶
The Foucault pendulum model is an example of a rigid system. The pendulum oscillates every few seconds (fast changing component), but it is also affected by the rotation of the Earth (slow changing component), with a periodicity of days, years, etc. The motion of the pendulum is dictated by inertia and the earth's gravitation. And the plane in which the oscillations occur, rotates synchronously with the daily rotation of the Earth.
We will specify the model through two differential equations that describe the motion of the weight at the end of the pendulum in the x-y plane (top view). We will evaluate the performance of different solvers by studying the law of conservation of total energy of the system.
The model does not take into account friction and air resistance, so the equations of motion look very concise:
$$\ddot x = 2 \Omega sin \lambda \dot y - \frac{g}{L}x$$
$$\ddot y = - 2 \Omega sin \lambda \dot x - \frac{g}{L} y$$
where $\Omega$ is the angular velocity of the Earth's rotation around its axis, $\lambda$ is the geographic latitude, $g$ is the acceleration of free fall, $x$ and $y$ are the coordinates of the centre of mass of the weight at the end of the pendulum.
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));
The model that realises these equations is shown below:
Variable pitch solvers¶
We are going to compare the following integrators (solvers):
| Name (Julia) | Name (Simulink) | Type of problem | Accuracy | When to use |
| ----------- | ----------- | ----------- | ----------- | ----------- |
| DP5 but better Tsit5 | ode45 | Not hard | Medium | Default. Try first |
| BS3 | ode23 | Not rigid | Low | If the task is error prone or for medium rigid systems |
| VCABM but better Vern7 | ode113 | Not hard | Low to high | Low error, suitable for computationally intensive tasks |
| QNDF, FBDF, but better Rodas4, KenCarp4, TRBDF2, RadauIIA5 | ode15s | Tough | Low to medium | If ode45 is too slow due to the hardness of the task |
| Rosenbrock23, but better Rodas4 | ode23s | Tough | Low | Low | When the accuracy requirements for hard problems are low and the [mass matrix] is unchanged (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 | Moderately stiff | Low | For moderately stiff problems, if you want a solution without numerical damping |
| TRBDF2 | ode23tb | Stiff | Low | For low accuracy requirements, for stiff problems |
If you don't find the right solver when solving problems and you are more used to Simulink solver names, take a look at the comparison from this article. Engee predominantly uses solver names adopted in Julia packages.
Comparison criteria¶
The quality of the solution of each integration algorithm can be evaluated in several ways. Firstly, if the problem has an analytical solution, we can compare the result of numerical integration with the theoretical result. It is also possible to estimate the time for which this or that algorithm solves the problem and use the time cost as a basis for comparison.
For the Foucault pendulum, alas, there is no exact analytical solution, and approximate ones have their own inherent errors.
Conservation of the total energy of the system¶
We will evaluate the quality of each solver's solution by checking whether the law of conservation of energy in the system is satisfied.
In the absence of friction, the total energy of the pendulum must remain constant. Numerical solutions will not show us perfect conservation of energy due to accumulation of errors.
So, in the example, we will calculate the total energy of the system at each integration step. Ideally, the total energy should not change at all and the difference should be zero. Total energy is the sum of potential and kinetic energy. It is this value that is calculated by the block 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)}$$
where $E(0)$ is the total energy at the initial moment of time.
The change of the total energy from the beginning to the end of the calculations will allow us to evaluate the quality of the numerical solution.
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;
Let us compare two integrators by this parameter:
- BS3 (
ode23) - rather coarse algorithm, unsuitable for rigid systems, - Trapezoid (
ode23t) - also rather coarse algorithm, but suitable for rigid systems.
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 )
As can be easily seen, the "losses" of total energy due to inaccuracies in the integrator are almost imperceptible in the first case, while in the second case they amount to about 20% per day.
As the total energy of the pendulum decreases, the amplitude of the oscillation drops, as we will see by plotting the following graphs. We will also see that the plane of oscillation of the pendulum changes with time (unfolds).
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) )
We have observed how the choice of solver determines the quality of the solution for rigid and non-rigid systems. The plane in which the pendulum oscillates does not make a complete reversal during the day: the angle of rotation depends on the geographical latitude.
We will construct a graph that will allow us to study in detail the speed and quality of different solvers on the task of modelling the Foucault pendulum. We will compare the results for different numerical integration methods and see how the energy of the system is conserved for different values of the relative error.
Running this script may take a few minutes, so to save time when running the script again you have the option of loading the data from a previously saved file 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
Now we can plot and analyse the graph.
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)
)
We see the zone where the non-prediction criterion of the total energy of the system practically ceases to change (at relative_tolerance < 1e-6).
This is due to the behaviour of the particular model. But, in general, it's worth keeping in mind that decreasing the relative accuracy of the solver does not necessarily lead to increased accuracy. When we set a desired limit on the accuracy of the solution and choose other solver settings, we should always think about the allowable simulation time.
We are likely to be satisfied with the accuracy of estimating the position of the weight at the end of the pendulum in a few centimetres. And it is no longer reasonable to calculate with accuracy to microns, because it will be very difficult to conduct the same accurate real experiment to validate the model.
Conclusion¶
We have studied how the choice of the numerical method for solving the equations of the system (the choice of the solver) affects the results and the progress of solving a rigid system. A lot depends on the solver: the algorithm should not only allow us to achieve the desired accuracy of the solution, but also to finish in a reasonable time.
The relative accuracy of a variable step solver (as well as the integration step of a constant step solver) must be small enough to prevent the solution from becoming unstable.
All solvers have their advantages and disadvantages, but some are known to work better when solving hard problems. And variable-pitch solvers are generally more stable than constant-pitch solvers.