Julia и фазовые портреты динамических систем

  • Tutorial


Продолжаем осваивать молодой и перспективный язык общего назначения Julia. Но для начала нам нужна как раз таки сугубо узкая возможность применения — для решения типичных задач физики. Это самая удобная тактика освоения инструмента: чтоб набить руку, будем решать насущные проблемы, постепенно увеличивая сложность и находя способы облегчения своей жизни. Короче, будем решать дифуры и строить графики.


Для начала качаем графический пакет Gadfly, да сразу полную версию разработчика, чтоб хорошо работала с нашей Julia 1.0.1. Вбиваем в нашем REPL, JUNO или Jupyter notebook:


using Pkg
pkg"add Compose#master"
pkg"add Gadfly#master"

Не самый удобный вариант, но пока ждём обновления PlotlyJS можно ради сравнения и попробовать.


Также нужен могущественнейший инструмент для решения дифференциальных уравнений DifferentialEquations


]add DifferentialEquations # альтернативный способ добавления пакета

Это самый обширный и поддерживаемый пакет, предоставляющий уйму методов решения дифуров. Теперь приступим к фазовым портретам.


Решения обыкновенных дифференциальных уравнений часто удобнее изображать не в привычном виде $y^1(t),...,y^L(t)$, а в фазовом пространстве, по осям которого откладываются значения каждой из найденных функций. При этом аргумент t входит в графики лишь параметрически. В случае двух ОДУ такой график – фазовый портрет системы, является кривой на фазовой плоскости и поэтому особенно нагляден.


Осциллятор


Динамика гармонического осциллятора с затуханием $\gamma$ описывается следующей системой уравнений:


$ \dot{x}=y\\ \dot{y}=-\omega^2x+\gamma y $


и независимо от начальных условий (х0, у0), приходит в состояние равновесия, точку (0,0) с нулевым углом отклонения и нулевой скоростью.


При $\gamma=0$ решение принимает вид свойственный для классического осциллятора.


using DifferentialEquations, Gadfly # включаем пакеты в проект - страшно долгий процесс

function garmosc(ω = pi, γ = 0, x0 = 0.1, y0 = 0.1)

    A = [0 1
        -ω^2 γ]

    syst(u,p,t) = A * u # ODE system

    u0 = [x0, y0] # start cond-ns
    tspan = (0.0, 4pi) # time period

    prob = ODEProblem(syst, u0, tspan) # problem to solve
    sol = solve(prob, RK4(),reltol=1e-6, timeseries_steps = 4)

    N = length(sol.u)
    J = length(sol.u[1])

    U = zeros(N, J)

    for i in 1:N, j in 1:J
        U[i,j] = sol.u[i][j] # перепишем ответы в матрицу поудобней
    end
    U
end

Разберем код. Функция принимает значения частоты и параметра затухания, а также начальные условия. Вложенная функция syst() задает систему. Для того чтоб вышло однострочно воспользовались матричным умножением. Функция solve() принимая огромное количество параметров весьма гибко настраивается для решения проблемы, но мы лишь указали метод решения — Рунге-Кутты 4 (есть ещё много других), относительную погрешность, и то, что надо сохранять не все точки, а только каждую четвертую. В переменную sol сохраняется матрица-ответ, причем, sol.t содержит вектор значений времен, а sol.u решения системы на этих временах. Все это спокойно рисуется в Plots, а для Gadfly придется укомплектовывать ответ в матрицу поудобней. Времена нам не нужны, так что возвращаем только решения.


Построим фазовый портрет:


Ans0 = garmosc()
plot(x = Ans0[:,1], y = Ans0[:,2])


Среди функций высокого порядка для нас особенно приметна broadcast(f, [1, 2, 3]), которая в функцию f поочередно подставляет значения массива [1, 2, 3]. Краткая запись f.([1, 2, 3]). Это пригодится чтобы поварьировать частоту, параметр затухания, начальную координату и начальную скорость.


Упрятано под спойлер
L = Array{Any}(undef, 3)# создали массив чего угодно (для слоев графиков)
clr = ["red" "green" "yellow"]

function ploter(Answ)
    for i = 1:3
        L[i] = layer(x = Answ[i][:,1], y = Answ[i][:,2], Geom.path, Theme(default_color=clr[i]) )
    end
    p = plot(L[1], L[2], L[3])
end

Ans1 = garmosc.( [pi 2pi 3pi] )# частоту
p1 = ploter(Ans1)

Ans2 = garmosc.( pi, [0.1 0.3 0.5] )# затухание
p2 = ploter(Ans2)

Ans3 = garmosc.( pi, 0.1, [0.2 0.5 0.8] ) 
p3 = ploter(Ans3)

Ans4 = garmosc.( pi, 0.1, 0.1, [0.2 0.5 0.8] )
p4 = ploter(Ans4)

set_default_plot_size(16cm, 16cm)

vstack( hstack(p1,p2), hstack(p3,p4) ) # слепим горизонтально и вертикально


Теперь рассмотрим не малые колебания, а колебания произвольной амплитуды:


$ \dot{x}=y\\ \dot{y}=-sin(\omega x)+\gamma y $


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


Упрятано под спойлер
function ungarmosc(ω = pi, γ = 0, x0 = 0.1, y0 = 0.1)

    function syst(du,u,p,t)
        du[1] = u[2]
        du[2] = -sin( ω*u[1] )+γ*u[2]
    end

    u0 = [x0, y0] # start cond-ns
    tspan = (0.0, 2pi) # time period

    prob = ODEProblem(syst, u0, tspan) # problem to solve
    sol = solve(prob, RK4(),reltol=1e-6, timeseries_steps = 4)

    N = length(sol.u)
    J = length(sol.u[1])

    U = zeros(N, J)

    for i in 1:N, j in 1:J
        U[i,j] = sol.u[i][j] # перепишем ответы в матрицу поудобней
    end
    U
end

Ans1 = ungarmosc.( [pi 2pi 3pi] )
p1 = ploter(Ans1)

Ans2 = ungarmosc.( pi, [0.1 0.4 0.7] )
p2 = ploter(Ans2)

Ans3 = ungarmosc.( pi, 0.1, [0.1 0.5 0.9] )
p3 = ploter(Ans3)

Ans4 = ungarmosc.( pi, 0.1, 0.1, [0.1 0.5 0.9] )
p4 = ploter(Ans4)

vstack( hstack(p1,p2), hstack(p3,p4) )


Брюсселятор


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


$ \dot{u}_1=-(\mu+1)u_1+u^2_1u_2+1 \\ \dot{u}_2=\mu u_1-u^2_1u_2 $


Неизвестные функции $u_1(t),\,u_2(t)$ отражают динамику концентрации промежуточных продуктов химической реакции. Параметр модели $\mu$ имеет смысл исходной концентрации катализатора (третьего вещества).


Более детально, эволюцию фазового портрета брюсселятора можно наблюдать, проводя расчеты с различным параметром $\mu$. При его увеличении узел будет сначала постепенно смещаться в точку с координатами (1, $\mu$), пока не достигнет бифуркационного значения $\mu$=2. В этой точке происходит качественная перестройка портрета, выражающаяся в рождении предельного цикла. При дальнейшем увеличении $\mu$ происходит лишь количественное изменение параметров этого цикла.


Упрятано под спойлер
function brusselator(μ = 0.1, u0 = [0.1; 0.1])

    function syst(du,u,p,t)
        du[1] = -(μ+1)*u[1] + u[1]*u[1]*u[2] + 1
        du[2] =      μ*u[1] - u[1]*u[1]*u[2]
    end

    tspan = (0.0, 10) # time period

    prob = ODEProblem(syst, u0, tspan) # problem to solve
    sol = solve(prob, RK4(),reltol=1e-6, timeseries_steps = 1)

    N = length(sol.u)
    J = length(sol.u[1])

    U = zeros(N, J)

    for i in 1:N, j in 1:J
        U[i,j] = sol.u[i][j] # перепишем ответы в матрицу поудобней
    end
    U
end

L = Array{Any}(undef, 10)

function ploter(Answ)
    for i = 1:length(Answ)
        L[i] = layer(x = Answ[i][:,1], y = Answ[i][:,2], Geom.path )
    end
    plot(L[1], L[2], L[3], L[4], L[5], L[6], L[7], L[8], L[9], L[10])
end

SC = [ [0 0.5], [0 1.5], [2.5 0], [1.5 0], [0.5 1], [1 0], [1 1], [1.5 2], [0.1 0.1], [0.5 0.2] ]

Ans1 = brusselator.( 0.1, SC )
Ans2 = brusselator.( 0.8, SC )
Ans3 = brusselator.( 1.6, SC )
Ans4 = brusselator.( 2.5, SC )
p1 = ploter(Ans1)
p2 = ploter(Ans2)
p3 = ploter(Ans3)
p4 = ploter(Ans4)

set_default_plot_size(16cm, 16cm)

vstack( hstack(p1,p2), hstack(p3,p4) )


На сегодня достаточно. В следующий раз попробуем научиться пользоваться другим графическим пакетом на новых задачах, попутно набивая руку на синтаксис Юлии. В процессе решения задач потихоньку начинают прослеживаться не то чтобы прям плюсы и минусы… скорее приятности и неудобства — этому следует посвятить отдельный разговор. И конечно хотелось бы на хабре увидеть что-то посерьезней моих игр с функциями — посему призываю джулистов поопытней рассказать о более менее серьезных проектах, это помогло бы многим и мне в частности, в изучении этого замечательного языка.

Поделиться публикацией

Комментарии 9

    0
    Хорошим тоном будет указывать версии всего своего хозяйства:
    Статус
    julia>]
    (v1.0) pkg> status
        Status `C:\Users\Игорь\.julia\environments\v1.0\Project.toml`
      [537997a7] AbstractPlotting v0.9.0
      [ad839575] Blink v0.8.1
      [159f3aea] Cairo v0.5.6
      [5ae59095] Colors v0.9.5
      [8f4d0f93] Conda v1.1.1
      [0c46a032] DifferentialEquations v5.3.1
      [a1bb12fb] Electron v0.2.0
      [5789e2e9] FileIO v1.0.2
      [5752ebe1] GMT v0.5.0
      [28b8d3ca] GR v0.35.0
      [c91e804a] Gadfly v1.0.0+ #master (https://github.com/GiovineItalia/Gadfly.jl.
    git)
      [4c0ca9eb] Gtk v0.16.4
      [a1b4810d] Hexagons v0.2.0
      [7073ff75] IJulia v1.13.0+ [`C:\Users\Игорь\.julia\dev\IJulia`]
      [6218d12a] ImageMagick v0.7.1
      [c601a237] Interact v0.9.0
      [b964fa9f] LaTeXStrings v1.0.3
      [ee78f7c6] Makie v0.9.0+ #master (https://github.com/JuliaPlots/Makie.jl.git)
      [7269a6da] MeshIO v0.3.1
      [47be7bcc] ORCA v0.1.1
      [58dd65bb] Plotly v0.2.0
      [f0f68f2c] PlotlyJS v0.12.0+ #master (https://github.com/sglyon/PlotlyJS.jl.gi
    t)
      [91a5bcdd] Plots v0.21.0
      [438e738f] PyCall v1.18.5
      [d330b81b] PyPlot v2.6.3
      [c4c386cf] Rsvg v0.2.2
      [60ddc479] StatPlots v0.8.1
      [b8865327] UnicodePlots v0.3.1
      [0f1e0344] WebIO v0.4.0
      [c2297ded] ZMQ v1.0.0
    

      +1
      стрелок нет, на фазовых портретах это важно
        0

        Объективности ради, в первых двух случаях не нужно (правило часовой стрелки). А в "Брюсселяторе" вообще оси неверно подписаны.

        0
        Уже стоит переходить с Fortran на Julia или подождать еще пару (десяти-)лет(-ий)?
        Будет эта Julia параллелится на 128 процессоров? а на 128000?
          0
          docs.julialang.org/en/v1/manual/parallel-computing
          Всё можно, если осторожно, а еще можно использовать CUDA, но лучше годок подождать, пока окружение подтянется, и документация расширится
            0
            Пробовал переходить с фортрана на Джулию, но как то все не получается.
            Вроде и синтаксис приятнее, и работает тоже быстро. Пробовал тестировать на примитивных задачках типа табуляции всяких спецфункций — вылезли какие то баги (я так понял конкретно в пакете спецфункций). В то же время реализация тех же спецфункций в Fortran IMSL — отлично сработала. Общий вывод — пока окружение сыроватое. Если не хотите взять участие в его тестировании и отладке, а может и исправлении — стоит подождать.
            +1
            using DifferentialEquations, Gadfly # включаем пакеты в проект — страшно долгий процесс

            «страшно долгий процесс» — это сарказм, или реальность?
              0
              Реальность (Jupyter задумывается на минуту на две). В оправдание можно было бы сказать, что DE очень сложный и обширный пакет, а Гадфлай мастер пока бета и вместе они притормаживают. В принципе, подключил пакеты и на протяжении сессии о них не думаешь — в юпитерском ноутбуке можно насобирать кучу однотипных программ и параллельно разбираться с ними, а не подключать пакеты для каждой. Но все равно это портит впечатление: тот же Scilab запускается не дольше 30 сек и уже сходу можно и дифур решить и нарисовать всё что нужно.
                +1
                спасибо.
                все как-то руки не доходят поставить какой-нибудь из этих пакетов да попробовать. По роду деятельности это не надо, а любопытство гложет. ибо «в наши времена такого не было»©

            Только полноправные пользователи могут оставлять комментарии. Войдите, пожалуйста.

            Самое читаемое