Модель маятника с неопределенностью по параметрам
В этой демонстрации мы изучим процесс запуска модели математического маятника, значения параметров которого заданы с разбросом.
Заодно мы проверим, какие особенности есть у библиотеки DifferentialEquations
(для решения ОДУ) и Plots
(вывод графиков) в плане реализации ключевой для Julia технологии мультиметодов (multiple dispatch). Той технологии, которая позволяет пакетам в Julia подбирать правильные методы для работы с самыми неожиданными типами данных.
Подготовка окружения
Для этой демонстрации нам потребуются библиотеки Measurements
и DifferentialEquations
. Если одна из них не установлена в вашей системе, наберите ]add Measurements
(и ]add DifferentialEquations
) в отдельной новой ячейке или в командной строке и выполните команду.
Pkg.add(["DifferentialEquations", "Measurements"])
using Measurements, DifferentialEquations
Первым делом зададим параметры окружения и самого маятника.
g = 9.81 ± 0.01; # Гравитационная постоянная
L = 1.00 ± 0.02; # Длина подвеса маятника
Символ допуска ±
можно ввести через команду LaTeX в командной строке: \pm
. После нажатия табуляции вы получите символ, который можно использовать в вычислениях и который интерпретируется в рамках некоторых библиотек, в частности Measurements
.
Роль библиотеки Measurements
состоит в том, чтобы правильно обрабатывать разбросы (допуски). Например, чтобы вычитание измерения из него же самого не приводило к удваиванию разброса ошибки.
x = 1 ± 0.01;
x - x
В этом плане сложение и вычитание работают неодинаково.
x + x
Постановка задачи о решении дифференциального уравнения
В этом примере мы реализуем упрощенную модель математического маятника:
где – постоянная гравитации, – длина маятника, – угол отклонения маятника в радианах. Применяется аппроксимация малых углов , поэтому считается что .
Зададим изначальную скорость и угол отклонения, а также временные границы симуляции.
u₀ = [0 ± 0.0, π/6 ± 0.01]
tspan = (0.0, 6.3)
Ноль в нижнем регистре можно ввести как команду LaTeX в командной строке: \_0
и нажать клавишу табуляции. После преобразования вы получаете новый символ, который можно скопировать в буфер обмена и использовать в названиях переменных.
Зададим уравнения качающегося маятника.
function pendulum( u, p, t )
θ = u[1]
dθ = u[2]
du = [dθ, (-g/L) * θ]
return du
end
Каждый вектор состоит из двух значений. Вектор состояния u
состоит из значения отклонения маятника θ и его скорости dθ. В свою очередь, вектор производных включает в себя скорость движения маятника dθ и силу, которая воздействует на груз.
Формирование вычислительного задания и запуск
prob = ODEProblem( pendulum, u₀, tspan )
Мы видим, что аргумент u
этой функции имеет тип Measurement
, аргумент t
имеет тип Float64
, и они отлично работают в связке. Теперь запустим симуляцию.
sol = solve( prob, Tsit5(), reltol=1e-6 );
plot( sol.t, first.(sol.u), label="положение" )
plot!( sol.t, last.(sol.u), label="скорость" )
Как мы видим, функция plot
автоматически проставляет доверительные границы, по сути, дополняя графики аргументом yerr
.
Аналитическое решение
Первым практическим применением этого подхода может быть сравнение ошибки численного и аналитического решений этой задачи.
u = u₀[2] .* cos.( sqrt( g/L ) .* sol.t );
plot( sol.t, last.(sol.u), label="Численным методом" )
plot!( sol.t, u, label="Аналитическим", title="Сравнение решений задачи моделирования маятника" )
Мы почти не видим разницы ни в значениях функции, ни в разбросе вокруг каждой точки обоих вариантов решения.
Замечания
На время публикации, текст справки по этой библиотеке еще не был включен в документацию 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 )
, происходит ошибка преобразования типов.
Заключение
Мы успешно нашли численное решение дифференциального уравнения и вывели график с допусками для задачи моделирования математического маятника, входные значения которой заданы в виде чисел с допусками по точности.