Документация Engee
Notebook

Модель маятника с неопределенностью по параметрам

Открыть пример в Engee

В этой демонстрации мы изучим процесс запуска модели математического маятника, значения параметров которого заданы с разбросом.

Заодно мы проверим, какие особенности есть у библиотеки DifferentialEquations (для решения ОДУ) и Plots (вывод графиков) в плане реализации ключевой для Julia технологии мультиметодов (multiple dispatch). Той технологии, которая позволяет пакетам в Julia подбирать правильные методы для работы с самыми неожиданными типами данных.

Подготовка окружения

Для этой демонстрации нам потребуются библиотеки Measurements и DifferentialEquations. Если одна из них не установлена в вашей системе, наберите ]add Measurements (и ]add DifferentialEquations) в отдельной новой ячейке или в командной строке и выполните команду.

In [ ]: 
using Measurements, DifferentialEquations

Первым делом зададим параметры окружения и самого маятника.

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

Символ допуска ± можно ввести через команду LaTeX в командной строке: \pm. После нажатия табуляции вы получите символ, который можно использовать в вычислениях и который интерпретируется в рамках некоторых библиотек, в частности Measurements.

Роль библиотеки Measurements состоит в том, чтобы правильно обрабатывать разбросы (допуски). Например, чтобы вычитание измерения из него же самого не приводило к удваиванию разброса ошибки.

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

В этом плане сложение и вычитание работают неодинаково.

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

Постановка задачи о решении дифференциального уравнения

В этом примере мы реализуем упрощенную модель математического маятника:

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

где $g$ – постоянная гравитации, $L$ – длина маятника, $\theta$ – угол отклонения маятника в радианах. Применяется аппроксимация малых углов $\theta \approx 0$, поэтому считается что $sin(\theta) \approx \theta$.

Зададим изначальную скорость и угол отклонения, а также временные границы симуляции.

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

Ноль в нижнем регистре можно ввести как команду LaTeX в командной строке: \_0 и нажать клавишу табуляции. После преобразования вы получаете новый символ, который можно скопировать в буфер обмена и использовать в названиях переменных.

Зададим уравнения качающегося маятника.

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)

Каждый вектор состоит из двух значений. Вектор состояния u состоит из значения отклонения маятника θ и его скорости dθ. В свою очередь, вектор производных включает в себя скорость движения маятника dθ и силу, которая воздействует на груз.

Формирование вычислительного задания и запуск

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

Мы видим, что аргумент u этой функции имеет тип Measurement, аргумент t имеет тип Float64, и они отлично работают в связке. Теперь запустим симуляцию.

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]:

Как мы видим, функция plot автоматически проставляет доверительные границы, по сути, дополняя графики аргументом yerr.

Аналитическое решение

Первым практическим применением этого подхода может быть сравнение ошибки численного и аналитического решений этой задачи.

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]:

Мы почти не видим разницы ни в значениях функции, ни в разбросе вокруг каждой точки обоих вариантов решения.

Замечания

На время публикации, текст справки по этой библиотеке еще не был включен в документацию Engee, поэтому приводим ссылку на страницу разработчиков: https://github.com/JuliaPhysics/Measurements.jl (которую на первое время можно перевести встроенными средствами браузера).

Особенности работы DifferentialEquations с пакетом Measurements

Мы реализовали функцию pendulum(u, p, t), которая использует return. Пакет DifferentialEquations может также реализовывать функции in-place, метод вызова которых будет выглядеть как pendulum(du, u, p, t), но из-за того, что мы не контролируем напрямую тип переменой du, распространение типа данных Measurement прерывается и этот вариант описания маятника не будет работать. Поэтому мы пользуемся return.

Особенности работы Plots с пакетом Measurements

Как видно, вывод первого из двух наших графиков тоже разделен на две команды. Функция plot, хотя и принимает на вход данные типа Measurement, не имеет конкретного рецепта для работы с несколькими векторами таких данных. Мы не можем успешно вызвать plot( sol.t, sol.u ), происходит ошибка преобразования типов.

Заключение

Мы успешно нашли численное решение дифференциального уравнения и вывели график с допусками для задачи моделирования математического маятника, входные значения которой заданы в виде чисел с допусками по точности.