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

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

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

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

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

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

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

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

Осциллятор

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

и независимо от начальных условий (х0, у0), приходит в состояние равновесия, точку (0,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) ) # слепим горизонтально и вертикально

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

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

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 г. Лефевром и Пригожиным

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

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

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) )

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