Engee documentation
Notebook

Pendulum model with uncertainty in parameters

In this demonstration we will study the process of running a model of a mathematical pendulum whose parameter values are specified with scatter.

At the same time, we will check what features the library DifferentialEquations (for solving ODEs) and Plots (plotting) have in terms of implementing the multimethods (multiple dispatch) technology that is key to Julia. The technology that allows packages in Julia to select the right methods to work with the most unexpected data types.

Preparing the environment

For this demonstration, we will need the libraries 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 at the command line and run the command.

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

The first thing to do is to set the parameters of the environment and the pendulum itself.

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

The tolerance symbol ± can be entered via the LaTeX command on the command line: \pm. After pressing tab, you will get a symbol that can be used in calculations and is interpreted within some libraries, in particular Measurements.

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

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

In this respect, addition and subtraction do not work in the same way.

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:

$$\ddot \theta = -\frac{g}{L} \theta$$

where $g$ is the gravitational constant, $L$ is the length of the pendulum, $\theta$ is the angle of deflection of the pendulum in radians. We use the approximation of small angles $\theta \approx 0$, so it is assumed that $sin(\theta) \approx \theta$.

Let us set the initial velocity and angle of deflection, 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 lower case 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 can be copied to the clipboard and used in variable names.

Let's set the equations of the 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 consists of the value of the pendulum's deflection θ and its velocity dθ. In turn, the vector of derivatives includes the velocity of the pendulum dθ and the force that acts on the weight.

Formation of the computational task and start

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 of this function has the type Measurement, the argument t has the type Float64, and they work well 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, essentially augmenting the graphs with the argument yerr.

Analytical solution

The first practical application of this approach may be to compare the error of the 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 function values or in the scatter around each point of both solutions.

Remarks

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

Features of DifferentialEquations with the Measurements package

We have implemented the function pendulum(u, p, t), which uses return. The DifferentialEquations package can also implement the function in-place, whose call method will look like pendulum(du, u, p, t), but because we do not directly control the variable type du, the propagation of the data type Measurement is interrupted and this variant of the pendulum description will not work. So we use return.

Features of Plots with the Measurements package

As you can see, the output of the first of our two plots is also divided into two commands. The function plot, although it accepts as input data of type Measurement, has no specific recipe for working with several vectors of such data. We cannot 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 tolerance plot for a mathematical pendulum modelling problem whose input values are given as numbers with tolerances in accuracy.