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

Изучаем решатели с переменным шагом

В этом примере мы разберем, как выбор решателя сказывается на вычислении траектории маятника Фуко по математической модели. Мы продемонстрируем работу решателей DP5 и Tsit5 (ode45), Rodas4 (ode15s), BS3 (ode23) а также Trapezoid (ode23t). В скобках приведены названия аналогичных решателей, принятые в Simulink.

Решение задачи будет задано в виде дифференциальных уравнений, имеющих черты жесткой системы. Для таких систем плохо подходят явные методы интегрирования (использующие данные пары соседних точек чтобы найти производную). Несравненно лучший результат обычно дают неявные методы, которые находят решение СЛАУ от данных нескольких соседних точек от искомой. Определение жестких систем довольно расплывчато, но есть одна общая черта их поведения – они могут требовать очень маленького шага интегрирования, без чего вычисления идут неустойчиво (например, результат уходит в бесконечность).

Жесткость часто присуща системам уравнений, состоящим из двух компонентов, состояние которых меняется с разной скоростью: очень быстро и очень медленно.

Маятник Фуко

Модель маятника Фуко – пример жесткой системы. Маятник колеблется раз в несколько секунд (быстро изменяющийся компонент), но на него также влияет вращение Земли (медленно изменяющийся компонент), периодичность которого исчисляется сутками, годами и т.д. Движение маятника продиктованы инерцией и земным тяготением. А плоскость, в которой происходят колебания, поворачивается синхронно с суточным вращением Земли.

Зададим модель через два дифференциальных уравнения, которые описывают движение груза на конце маятника в плоскости x-y (вид сверху). Оценку работы разных решателей мы будем производить, изучая закон сохранения полной энергии системы.

Модель не учитывает трение и сопротивление воздуха, поэтому уравнения движения выглядят очень лаконично:

$$\ddot x = 2 \Omega sin \lambda \dot y - \frac{g}{L}x$$

$$\ddot y = - 2 \Omega sin \lambda \dot x - \frac{g}{L} y$$

где $\Omega$ – угловая скорость вращения Земли вокруг оси, $\lambda$ – географическая широта, $g$ – ускорение свободного падения, $x$ и $y$ – координаты центра масс груза на конце маятника.

In [ ]:
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));

Модель, которая реализует эти уравнения, показана ниже:

image.png

Решатели с переменным шагом

Мы собираемся сравнить следующие интеграторы (решатели):

Название (Julia) Название (Simulink) Вид задач Точность Когда пользоваться
DP5, но лучше Tsit5 ode45 Не жесткие Средняя По умолчанию. Попробовать первым
BS3 ode23 Не жесткие Низкая Если в задаче допустимы большие ошибки или для средне жестких систем
VCABM, но лучше Vern7 ode113 Не жесткие От низкой до высокой Низкая ошибка, подходит для задач с большим объемом вычислений
QNDF, FBDF, но лучше Rodas4, KenCarp4, TRBDF2, RadauIIA5 ode15s Жесткие От низкой до средней Если ode45 работает слишком медленно из-за жесткости задачи
Rosenbrock23, но лучше Rodas4 ode23s Жесткие Низкая При низких требованиях к точности решения жестких задач и при неизменной матрице масс
Trapezoid ode23t Умеренно жесткие Низкая Для умеренно жестких задач, если нужно решение без численного демпфирования
TRBDF2 ode23tb Жесткие Низкая При низких требованиях к точности решения, для решения жестких задач

Если при решении задач вы не находите нужного решателя, при этом вам более привычны названия решателей Simulink – обратите внимание на сравнение из этой статьи. Engee преимущественно использует названия решателей, принятые в пакетах Julia.

Критерии сравнения

Качество решения каждого алгоритма интегрирования можно оценить несколькими способами. Во-первых, если у задачи есть аналитическое решение, то можно сравнить результат численного интегрирования с теоретическим результатом. Еще можно оценить время, за которое тот или иной алгоритм решает задачу и использовать временные затраты как базис для сравнения.

Для маятника Фуко, увы, нет точного аналитического решения, а приблизительные имеют собственные, свойственные им погрешности.

Сохранение полной энергии системы

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

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

Итак, в примере мы рассчитаем полную энергию системы на каждом шаге интегрирования. В идеале полная энергия вообще не должна меняться и разница должна быть равна нулю. Полная энергия – сумма потенциальной и кинетической энергии. Именно это значение и рассчитывает блок NormalizeEnergy:

$$E = \frac{v_x^2 + v_y^2}{2} + g(L - \sqrt{L^2 - x^2 - y^2})$$

$$E_{norm} (t) = \frac{E(t)}{E(0)}$$

где $E(0)$ – полная энергия в начальный момент времени.

Изменение полной энергии от момента начала до окончания расчетов и позволит нам оценить качество численного решения.

In [ ]:
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) – тоже довольно грубый алгоритм, но зато подходящий для жестких систем.
In [ ]:
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 )
Out[0]:
No description has been provided for this image

Как легко заметить, "потери" полной энергии, связанные с неточностями в интеграторе, в первом случае почти незаметны, во втором – составляют около 20% за сутки.

При снижении полной энергии у маятника падает амплитуда колебаний, в чем мы убедимся, построив следующие графики. Также мы увидим, что плоскость колебаний маятника меняется со временем (разворачивается).

In [ ]:
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) )
Out[0]:
No description has been provided for this image

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

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

Выполнение этого скрипта может занять несколько минут, поэтому чтобы сэкономить время при повторном запуске скрипта у вас есть возможность загрузки данных из ранее сохраненного файла results.bin.

In [ ]:
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
Out[0]:
(["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])

Теперь мы можем построить и проанализировать график.

In [ ]:
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)
    )
Out[0]: