Julia and phase portraits of dynamic systems

Julia and phase portraits of dynamic systems We continue to learn a young and promising general-purpose language Julia . But first, we need just the very narrow possibility of application - for solving typical problems of physics. This is the most convenient tactic for mastering an instrument: in order to fill our hands, we will solve urgent problems, gradually increasing complexity and finding ways to make our lives easier. In short, we will solve diffs and build graphs. 3r33333. Gadfly , but immediately the full version of the developer to work well with our Julia ???. We drive in our REPL, JUNO or Jupyter notebook:
3r33385.  3r33394.
3r33333. using Pkg
pkg "add Compose # master"
pkg "add Gadfly # master"
3r33385.  3r33394.
Not the most convenient option, but while waiting for the update 3r33333. PlotlyJS
can be for the sake of comparison and try. 3r33333. 3r33385.  3r33394.
You also need a powerful tool for solving differential equations 3r-338. DifferentialEquations
3r33333. 3r33385.  3r33394.
3r33333. ]add DifferentialEquations # An alternative way to add a package.
3r33385.  3r33394.
This is the most extensive and supported package, providing a lot of methods for solving difures. Now we proceed to the phase portraits. 3r33333. 3r33385.  3r33394.
The solutions of ordinary differential equations are often more conveniently portrayed in a different way than 3r-3319. 3r3355.
, and in
phase space
, along whose axes the values ​​of each of the functions found are deposited. The argument t is included in the graphs only parametrically. In the case of two ODEs, such a graph, the phase portrait of the system, is a curve on the phase plane and is therefore particularly visual. 3r33333. 3r33385.  3r33394. 3r362. Oscillator
3r33385.  3r33394.
Dynamics of a harmonic oscillator with attenuation

described by the following system of equations:
3r33385.  3r33394.
3r33333.
3r33333.
3r33333. 3r33385.  3r33394.
and regardless of the initial conditions (x? y0), it comes to equilibrium, a point (?0) with a zero angle of deviation and zero speed. 3r33333. 3r33385.  3r33394.
At
3r388.
The decision takes the form peculiar to the classical oscillator. 3r33333. 3r33385.  3r33394.
3r33333. using DifferentialEquations, Gadfly # we include packages in the project - a terribly long process
3r33394. function garmosc (ω = pi, γ = ? x0 = 0.? y0 = 0.1)
3r33394. A =[0 1
-ω^2 γ]3r33394. 3r33394. syst (u, p, t) = A * u # ODE system
3r33394. u0 =[x0, y0]# start cond-ns
tspan = (0.? 4pi) # time period
3r33394. prob = ODEProblem (syst, u? tspan) # problem to solve
sol = solve (prob, RK4 (), reltol = 1e-? timeseries_steps = 4)
3r33394. N = length (sol.u)
J = length (sol. U[1])
3r33394. U = zeros (N, J)
3r33394. for i in 1: N, j in 1: J
U[i,j]= sol.u[i] [j]# rewrite the answers in a matrix
more comfortably. end
U
end
3r33385.  3r33394.
Let's sort the code. The function takes the frequency and attenuation parameter values ​​as well as the initial conditions. Nested function
syst ()
sets the system. In order to come out single-line, they used matrix multiplication. Function
solve ()
taking a huge number of parameters is very flexible to adjust the problem, but we only indicated the method of solution - Runge-Kutta 4 (3r3128. there are many other 3r-3129.), relative error, and the fact that not all points need to be saved, but only every fourth . In the variable
sol
the answer matrix is ​​preserved, and
sol.t
contains a vector of time values, and
sol.u
system solutions at these times. All this is quietly drawn in
Plots
, and for
Gadfly
will have to complete the answer in a matrix more comfortably. We do not need times, so we only return solutions. 3r33333. 3r33385.  3r33394.
Construct the phase portrait:
3r33385.  3r33394.
3r33333. Ans0 = garmosc ()
plot (x = Ans0[:,1], y = Ans0[:,2])
3r33385.  3r33394.
Among the functions of a high order for us is particularly noticeable 3r3333330. broadcast (f,[1, 2, 3]) which is in function
f
alternately substitutes the values ​​of the array[1, 2, 3]. Short record f. ([1, 2, 3]) 3r33376. . This is useful to vary the frequency, attenuation parameter, initial coordinate, and initial velocity. 3r33333. 3r33385.  3r33394.
Hid under the spoiler [/b]
  3r33333. L = Array {Any} (undef, 3) # created an array of anything (for graph layers) 3r3r94. clr =["red" "green" "yellow"]3r33394. 3r33394. function ploter (Answers)
for i = 1: 3
L[i]= layer (x = Answ[i] [:,1], y = Answers[i] [:,2], Geom.path, Theme (default_color = clr[i]))
end
p = plot (L[1], L[2], L[3])
end
3r33394. Ans1 = garmosc. (W2w2w235.) # Frequency
p1 = ploter (Ans1)
3r33394. Ans2 = garmosc. (Pi,[0.1 0.3 0.5]) # Attenuation
p2 = ploter (Ans2)
3r33394. Ans3 = garmosc. (Pi, 0.?[0.2 0.5 0.8])
p3 = ploter (Ans3)
3r33394. Ans4 = garmosc. (Pi, 0.? 0.?[0.2 0.5 0.8])
p4 = ploter (Ans4)
3r33394. set_default_plot_size (16cm, 16cm)
3r33394. vstack (hstack (p? p2), hstack (p? p4)) # blind horizontally and vertically
3r33385.  3r33394.

3r3204. 3r33333. 3r33395. 3r33395. 3r33385.  3r33394.

Now we will consider not small oscillations, but oscillations of an arbitrary amplitude: 3r-3388. 3r33385.  3r33394.

3r33333.

3r33333. 3r33333. 3r33385.  3r33394.

The graph of the solution on the phase plane is not an ellipse (which indicates the nonlinearity of the oscillator). The less we choose the amplitude of oscillations (i.e., the initial conditions), the less non-linearity will appear (therefore, the small oscillations of the pendulum can be considered harmonic). 3r33333. 3r33385.  3r33394.

Hid under the spoiler [/b]
  3r33333. function ungarmosc (ω = pi, γ = ? x0 = 0.? y0 = 0.1)
3r33394. function syst (du, u, p, t)
du[1]= u[2]3r33394. du[2]= -sin (ω * u[1]) + γ * u[2]3r33394. end
3r33394. u0 =[x0, y0]# start cond-ns
tspan = (0.? 2pi) # time period
3r33394. prob = ODEProblem (syst, u? tspan) # problem to solve
sol = solve (prob, RK4 (), reltol = 1e-? timeseries_steps = 4)
3r33394. N = length (sol.u)
J = length (sol. U[1])
3r33394. U = zeros (N, J)
3r33394. for i in 1: N, j in 1: J
U[i,j]= sol.u[i] [j]# rewrite the answers in a matrix
more comfortably. end
U
end
3r33394. Ans1 = ungarmosc. (W2w2w235.)
p1 = ploter (Ans1)
3r33394. Ans2 = ungarmosc. (Pi,[0.1 0.4 0.7])
p2 = ploter (Ans2)
3r33394. Ans3 = ungarmosc. (Pi, 0.?[0.1 0.5 0.9])
p3 = ploter (Ans3)
3r33394. Ans4 = ungarmosc. (Pi, 0.? 0.?[0.1 0.5 0.9])
p4 = ploter (Ans4)
3r33394. vstack (hstack (p? p2), hstack (p? p4))
3r33385.  3r33394.

Brusselator 3r33385.  3r33394.

This model describes some autocatalytic chemical reaction in which diffusion plays a certain role. The model was proposed in 1968 by Lefebvre and Prigogine 3r-3388. 3r33385.  3r33394.

3r33333.

3r33333. 3r33333. 3r33385.  3r33394.

Unknown features 3r33232. reflect the dynamics of the concentration of intermediate products of a chemical reaction. Parameter model 3r33320. makes sense of the initial concentration of the catalyst (the third substance). 3r33333. 3r33385.  3r33394.

In more detail, the evolution of the phase portrait of the Brusselator can be observed by performing calculations with a different parameter 3r-3319. 3r33320. . When it is increased, the node will first gradually shift to a point with coordinates (? 3r33333. ) Until it reaches 3r3333314. bifurcation r3r3315. values. 3r33320. = 2. At this point, a qualitative reorganization of the portrait takes place, which is expressed in the birth of the limiting cycle. With a further increase of 3r33333. 3r33320. there is only a quantitative change in the parameters of this cycle. 3r33333. 3r33385.  3r33394.

Hid under the spoiler [/b]
  3r33333. function brusselator (μ = 0.? u0 =[0.1; 0.1])
3r33394. 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]3r33394. end
3r33394. tspan = (0.? 10) # time period
3r33394. prob = ODEProblem (syst, u? tspan) # problem to solve
sol = solve (prob, RK4 (), reltol = 1e-? timeseries_steps = 1)
3r33394. N = length (sol.u)
J = length (sol. U[1])
3r33394. U = zeros (N, J)
3r33394. for i in 1: N, j in 1: J
U[i,j]= sol.u[i] [j]# rewrite the answers in a matrix
comfortably. end
U
end
3r33394. L = Array {Any} (undef, 10)
3r33394. function ploter (Answers)
for i = 1: length (Answers)
L[i]= layer (x = Answ[i] [:,1], y = Answers[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 w222. L end
3r33394. 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]]
3r33394. Ans1 = brusselator. (0.? SC)
Ans2 = brusselator. (0.? SC)
Ans3 = brusselator. (1.? SC)
Ans4 = brusselator. (2.? SC)
p1 = ploter (Ans1)
p2 = ploter (Ans2)
p3 = ploter (Ans3)
p4 = ploter (Ans4)
3r33394. set_default_plot_size (16cm, 16cm)
3r33394. vstack (hstack (p? p2), hstack (p? p4))
3r33385.  3r33394.
3r33381. 3r33333. 3r33395. 3r33395. 3r33385.  3r33394.
It's enough for today. Next time, we’ll try to learn how to use another graphic package on new tasks, at the same time getting a hand on Julia’s syntax. In the process of solving problems, they slowly begin to be traced not so directly to the pros and cons rather, to pleasantness and inconvenience — a separate conversation should be devoted to this. And of course, I would like to see something more serious in my games with functions - so I urge julists to tell more about serious projects, this would help many people and me in particular in learning this wonderful language. 3r33333. 3r33395. 3r33394. 3r33394.
! function (e) {function t (t, n) {if (! (n in e)) {for (var r, a = e.document, i = a.scripts, o = i.length; o-- ;) if (-1! == i[o].src.indexOf (t)) {r = i[o]; break} if (! r) {r = a.createElement ("script"), r.type = "text /jаvascript", r.async =! ? r.defer =! ? r.src = t, r.charset = "UTF-8"; var d = function () {var e = a.getElementsByTagName ("script")[0]; e.parentNode.insertBefore (r, e)}; "[object Opera]" == e.opera? a.addEventListener? a.addEventListener ("DOMContentLoaded", d,! 1): e.attachEvent ("onload", d ): d ()}}} t ("//mediator.mail.ru/script/2820404/"""_mediator") () (); 3r33393. 3r33394. 3r33395.
+ 0 -

Add comment