Engee documentation
Notebook

A model of a pendulum with uncertainty in parameters

In this demonstration, we will study the process of starting a mathematical pendulum model, the values of which are set with a spread.

At the same time, we will check what features the library has. DifferentialEquations (to solve the ODE) and Plots (graph output) in terms of implementing Julia's key technology multimethods (multiple dispatch). The technology that allows Julia packages to select the right methods for working with the most unexpected types of data.

Environment preparation

We will need libraries for this demonstration. Measurements and DifferentialEquations. If one of them is not installed on your system, type ]add Measurements (and ]add DifferentialEquations) in a separate new cell or in the command prompt and run the command.

In [ ]: 
Pkg.add(["DifferentialEquations", "Measurements"])
In [ ]: 
using Measurements, DifferentialEquations

First of all, we will set the parameters of the environment and the pendulum itself.

In [ ]: 
g = 9.81 ± 0.01; # Гравитационная постоянная 
L = 1.00 ± 0.02; # Длина подвеса маятника

Access symbol ± you can enter it using the LaTeX command on the command line: \pm. After pressing the tab, you will receive a symbol that can be used in calculations and which is interpreted within some libraries, in particular Measurements.

The role of the library Measurements The goal is to handle the spreads (tolerances) correctly. For example, so that subtracting a measurement from the same measurement does not double the error spread.

In [ ]: 
x = 1 ± 0.01;
x - x
Out[0]:
$0.0 \pm 0.0$

In this regard, addition and subtraction work differently.

In [ ]: 
x + x
Out[0]:
$2.0 \pm 0.02$

Statement of the problem of solving a differential equation

In this example, we implement a simplified model of a mathematical pendulum.:

where – gravity constant, – the length of the pendulum, – the angle of deflection of the pendulum in radians. Small angle approximation is used , therefore , it is assumed that .

We will set the initial velocity and deflection angle, as well as the time limits of the simulation.

In [ ]: 
u₀ = [0 ± 0.0, π/6 ± 0.01]
tspan = (0.0, 6.3)
Out[0]:
(0.0, 6.3)

Zero in lowercase can be entered as a LaTeX command on the command line.: \_0 and press the tab key. After the conversion, you get a new character that you can copy to the clipboard and use in variable names.

Let's define the equations of a swinging pendulum.

In [ ]: 
function pendulum( u, p, t )
    θ = u[1]
    dθ = u[2]
    du = [dθ, (-g/L) * θ]
    return du
end
Out[0]:
pendulum (generic function with 1 method)

Each vector consists of two values. The state vector u It consists of the value of the deflection of the pendulum θ and its velocity dθ. In turn, the vector of derivatives includes the velocity of the pendulum Dθ and the force that acts on the load.

Creating a computational task and launching it

In [ ]: 
prob = ODEProblem( pendulum, u₀, tspan )
Out[0]:
ODEProblem with uType Vector{Measurement{Float64}} and tType Float64. In-place: false
timespan: (0.0, 6.3)
u0: 2-element Vector{Measurement{Float64}}:
   0.0 ± 0.0
 0.524 ± 0.01

We see that the argument u this function has the type Measurement, the argument t has a type Float64 and they work great together. Now let's run the simulation.

In [ ]: 
sol = solve( prob, Tsit5(), reltol=1e-6 );
In [ ]: 
plot( sol.t, first.(sol.u), label="положение" )
plot!( sol.t, last.(sol.u), label="скорость" )
Out[0]:

As we can see, the function plot automatically sets confidence bounds, in fact, supplementing the graphs with an argument yerr.

Analytical solution

The first practical application of this approach may be to compare the error of numerical and analytical solutions to this problem.

In [ ]: 
u = u₀[2] .* cos.( sqrt( g/L ) .* sol.t );
In [ ]: 
plot( sol.t, last.(sol.u), label="Численным методом" )
plot!( sol.t, u, label="Аналитическим", title="Сравнение решений задачи моделирования маятника" )
Out[0]:

We see almost no difference either in the values of the function or in the spread around each point of both solutions.

Remarks

At the time of publication, the help text for this library had not yet been included in the Engee documentation, so here is a link to the developers' page: https://github.com/JuliaPhysics/Measurements.jl (which, for the first time, can be translated using the browser's built-in tools).

Features of DifferentialEquations with the Measurements package

We have implemented the function pendulum(u, p, t), which uses return. Package DifferentialEquations It can also implement functions in-place, the method of calling which will look like pendulum(du, u, p, t) but due to the fact that we do not directly control the type of change du, data type propagation Measurement it is interrupted and this version of the description of the pendulum will not work. That's why we use return.

Features of how Plots work with the Measurements package

As you can see, the output of the first of our two graphs is also divided into two commands. Function plot, although it accepts input data of the type Measurement, does not have a specific recipe for working with multiple vectors of such data. We are unable to successfully call plot( sol.t, sol.u ), a type conversion error occurs.

Conclusion

We have successfully found a numerical solution to the differential equation and derived a graph with tolerances for the mathematical pendulum modeling problem, the input values of which are set as numbers with precision tolerances.