Изучаем решатели с переменным шагом
В этом примере мы разберем, как выбор решателя сказывается на вычислении траектории маятника Фуко по математической модели. Мы продемонстрируем работу решателей DP5 и Tsit5 (ode45), Rodas4 (ode15s), BS3 (ode23) а также Trapezoid (ode23t). В скобках приведены названия аналогичных решателей, принятые в Simulink.
Решение задачи будет задано в виде дифференциальных уравнений, имеющих черты жесткой системы. Для таких систем плохо подходят явные методы интегрирования (использующие данные пары соседних точек чтобы найти производную). Несравненно лучший результат обычно дают неявные методы, которые находят решение СЛАУ от данных нескольких соседних точек от искомой. Определение жестких систем довольно расплывчато, но есть одна общая черта их поведения – они могут требовать очень маленького шага интегрирования, без чего вычисления идут неустойчиво (например, результат уходит в бесконечность).
Жесткость часто присуща системам уравнений, состоящим из двух компонентов, состояние которых меняется с разной скоростью: очень быстро и очень медленно.
Маятник Фуко
Модель маятника Фуко – пример жесткой системы. Маятник колеблется раз в несколько секунд (быстро изменяющийся компонент), но на него также влияет вращение Земли (медленно изменяющийся компонент), периодичность которого исчисляется сутками, годами и т.д. Движение маятника продиктованы инерцией и земным тяготением. А плоскость, в которой происходят колебания, поворачивается синхронно с суточным вращением Земли.
Зададим модель через два дифференциальных уравнения, которые описывают движение груза на конце маятника в плоскости x-y (вид сверху). Оценку работы разных решателей мы будем производить, изучая закон сохранения полной энергии системы.
Модель не учитывает трение и сопротивление воздуха, поэтому уравнения движения выглядят очень лаконично:
где – угловая скорость вращения Земли вокруг оси, – географическая широта, – ускорение свободного падения, и – координаты центра масс груза на конце маятника.
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));Модель, которая реализует эти уравнения, показана ниже:
Решатели с переменным шагом
Мы собираемся сравнить следующие интеграторы (решатели):
| Название (Julia) | Название (Simulink) | Вид задач | Точность | Когда пользоваться |
|---|---|---|---|---|
DP5, но лучше Tsit5 |
|
Не жесткие |
Средняя |
По умолчанию. Попробовать первым |
BS3 |
|
Не жесткие |
Низкая |
Если в задаче допустимы большие ошибки или для средне жестких систем |
VCABM, но лучше Vern7 |
|
Не жесткие |
От низкой до высокой |
Низкая ошибка, подходит для задач с большим объемом вычислений |
QNDF, FBDF, но лучше Rodas4, KenCarp4, TRBDF2, RadauIIA5 |
|
Жесткие |
От низкой до средней |
Если |
Rosenbrock23, но лучше Rodas4 |
|
Жесткие |
Низкая |
При низких требованиях к точности решения жестких задач и при неизменной матрице масс |
Trapezoid |
|
Умеренно жесткие |
Низкая |
Для умеренно жестких задач, если нужно решение без численного демпфирования |
TRBDF2 |
|
Жесткие |
Низкая |
При низких требованиях к точности решения, для решения жестких задач |
Если при решении задач вы не находите нужного решателя, при этом вам более привычны названия решателей Simulink – обратите внимание на сравнение из этой статьи. Engee преимущественно использует названия решателей, принятые в пакетах Julia.
Критерии сравнения
Качество решения каждого алгоритма интегрирования можно оценить несколькими способами. Во-первых, если у задачи есть аналитическое решение, то можно сравнить результат численного интегрирования с теоретическим результатом. Еще можно оценить время, за которое тот или иной алгоритм решает задачу и использовать временные затраты как базис для сравнения.
Для маятника Фуко, увы, нет точного аналитического решения, а приблизительные имеют собственные, свойственные им погрешности.
Сохранение полной энергии системы
Мы будем оценивать качество решения каждого решателя, проверяя соблюдение закона сохранения энергии в системе.
При отсутствии трения полная энергия маятника должна оставаться постоянной. Численные решения не продемонстрируют нам идеального сохранения энергии из-за накопления ошибок.
Итак, в примере мы рассчитаем полную энергию системы на каждом шаге интегрирования. В идеале полная энергия вообще не должна меняться и разница должна быть равна нулю. Полная энергия – сумма потенциальной и кинетической энергии. Именно это значение и рассчитывает блок NormalizeEnergy:
где – полная энергия в начальный момент времени.
Изменение полной энергии от момента начала до окончания расчетов и позволит нам оценить качество численного решения.
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;Сравним по этому параметру два интегратора:
-
BS3 (
ode23) – довольно грубый алгоритм, неподходящий для жестких систем, -
Trapezoid (
ode23t) – тоже довольно грубый алгоритм, но зато подходящий для жестких систем.
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 )
Как легко заметить, "потери" полной энергии, связанные с неточностями в интеграторе, в первом случае почти незаметны, во втором – составляют около 20% за сутки.
При снижении полной энергии у маятника падает амплитуда колебаний, в чем мы убедимся, построив следующие графики. Также мы увидим, что плоскость колебаний маятника меняется со временем (разворачивается).
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) )
Мы пронаблюдали, как от выбора решателя зависит качество решения для жестких и нежестких систем. Плоскость, в которой происходят колебания маятника, не совершает полного разворота в ходе дня: угол поворота зависит от географической широты.
Построим график, который позволит нам в деталях изучить скорость и качество работы разных решателей на задаче моделирования маятника Фуко. Мы сравним результаты для разных методов численного интегрирования и посмотрим, как сохраняется энергия системы при разных значениях относительной ошибки.
Выполнение этого скрипта может занять несколько минут, поэтому чтобы сэкономить время при повторном запуске скрипта у вас есть возможность загрузки данных из ранее сохраненного файла 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(["Tsit5" "DP5" … "Trapezoid" "BS3"], [1.0000000000000001e-7, 3.162277660168379e-7, 1.0e-6, 3.162277660168379e-6, 1.0e-5, 3.1622776601683795e-5, 0.0001, 0.00031622776601683794, 0.001], [7.929647437576737e-5 7.993461600175333e-5 … 0.0008342080580365785 0.0010495257229566901; 0.006522783592120329 0.006803880028364184 … 0.0432470795415133 0.04324716090734526; … ; 0.00020914013837636247 0.0002237451074058594 … 0.008214673902898095 0.01263204365201287; 0.07957109906963433 0.08391745090197478 … 0.8201637378884246 0.957540689453349], [25.203170908 25.070636927 … 24.543856697 24.405600441; 24.731999073 24.538390433 … 23.893544084 24.32454816; … ; 30.82383809 31.923476917 … 24.9366935 24.482693973; 28.296808969 28.368053586 … 25.189954714 24.714920212])Теперь мы можем построить и проанализировать график.
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)
)
Мы видим зону, где критерий недопрогноза полной энергии системы практически перестает меняться (при relative_tolerance < 1e-6).
Это связано с поведением конкретной модели. Но, в целом, стоит учитывать, что уменьшение относительной точности решателя не обязательно приводит к повышению точности. Когда мы задаемся желаемым пределом точности решения и выбираем другие настройки решателя, всегда следует думать о допустимых временных затрат на симуляцию.
Нас скорее всего устроит точность оценки положения груза на конце маятника в несколько сантиметров. А выполнять расчет с точностью до микрон уже нецелесообразно, поскольку провести столь же точный реальный эксперимент для валидации модели будет очень затруднительно.
Заключение
Мы изучили, как выбор численного метода решения уравнений системы (выбор решателя) влияет на результаты и ход решения жесткой системы. От решателя зависит очень многое: ведь алгоритм не только долен позволить нам добиться нужной точности решения, но и завершиться за допустимое время.
Относительная точность решателя с переменным шагом (равно как и же шаг интегрирования решателя с постоянным шагом) должны быть достаточно малы, чтобы решение не становилось нестабильным.
У всех решателей есть свои преимущества и недостатки, но некоторые заведомо лучше работают при решении жестких задач. А решатели с переменным шагом в целом более стабильны, чем решатели с постоянным шагом.