Fig. 7.31 - MMSB

Hasty synthetic oscillator model

using OrdinaryDiffEq
using ComponentArrays: ComponentArray
using SimpleUnPack
using CairoMakie

function model731!(du, u, p, t)
    @unpack alpha, sigma, gammax, gammay, ay, D = p
    @unpack x1, y1, x2, y2 = u
    den1 = (1 + x1^2 + sigma * x1^4) * (1 + y1^4)
    den2 = (1 + x2^2 + sigma * x2^4) * (1 + y2^4)
    du.x1 = (1 + x1^2 + alpha * sigma * x1^4) / den1 - gammax * x1 + D * (x2 - x1)
    du.y1 = ay * (1 + x1^2 + alpha * sigma * x1^4) / den1 - gammay * y1
    du.x2 = (1 + x2^2 + alpha * sigma * x2^4) / den2 - gammax * x2 + D * (x1 - x2)
    du.y2 = ay * (1 + x2^2 + alpha * sigma * x2^4) / den2 - gammay * y2
    nothing
end

ps731 = ComponentArray(
    alpha=11,
    sigma=2,
    gammax=0.2,
    gammay=0.012,
    ay=0.2,
    D=0.015,
)

u0731 = ComponentArray(
    x1=0.3963,
    y1=2.3346,
    x2=0.5578,
    y2=1.9317,
)

tend = 500.0
prob731 = ODEProblem(model731!, u0731, (0.0, tend), ps731)
@time sol731 = solve(prob731, KenCarp47())

  2.382448 seconds (9.39 M allocations: 631.638 MiB, 7.74% gc time, 99.94% compilation time)

retcode: Success
Interpolation: 3rd order Hermite
t: 57-element Vector{Float64}:
   0.0
   0.11740367751851252
   1.2914404527036376
   4.511907461040367
   7.7323744693770955
  12.902005480058696
  18.281725612831522
  25.964102328817496
  35.37767386279916
  44.79124539678083
   ⋮
 440.82866211510526
 448.7663348469433
 453.2647770327574
 459.7155510928191
 464.72930124851143
 473.9750360356151
 483.2207708227188
 494.47524985828164
 500.0
u: 57-element Vector{ComponentArrays.ComponentVector{Float64, Vector{Float64}, Tuple{ComponentArrays.Axis{(x1 = 1, y1 = 2, x2 = 3, y2 = 4)}}}}:
 ComponentVector{Float64}(x1 = 0.3963, y1 = 2.3346, x2 = 0.5578, y2 = 1.9317)
 ComponentVector{Float64}(x1 = 0.392700105296416, y1 = 2.3323867602242525, x2 = 0.5624596801401064, y2 = 1.9325985904186846)
 ComponentVector{Float64}(x1 = 0.36048669011930534, y1 = 2.3101526487371995, x2 = 0.6157022192530186, y2 = 1.944303374291833)
 ComponentVector{Float64}(x1 = 0.30633703482186236, y1 = 2.2498619616642235, x2 = 0.8526128340898065, y2 = 2.0121405680489555)
 ComponentVector{Float64}(x1 = 0.2945818938260461, y1 = 2.193069032689244, x2 = 1.1689193875728283, y2 = 2.1331175922275176)
 ComponentVector{Float64}(x1 = 0.321919289906436, y1 = 2.1130104722389222, x2 = 1.3259655230067424, y2 = 2.308056546511433)
 ComponentVector{Float64}(x1 = 0.3671580012546514, y1 = 2.047103971983285, x2 = 1.166571170312196, y2 = 2.4092074186399257)
 ComponentVector{Float64}(x1 = 0.45245387780726537, y1 = 1.9897408512629327, x2 = 0.8305956994535477, y2 = 2.438956766773974)
 ComponentVector{Float64}(x1 = 0.7074438091524423, y1 = 2.021692292452505, x2 = 0.4329553774603207, y2 = 2.3252281436364735)
 ComponentVector{Float64}(x1 = 1.2148531958524662, y1 = 2.2779845345881142, x2 = 0.320691587147811, y2 = 2.1630507293442416)
 ⋮
 ComponentVector{Float64}(x1 = 0.36392421685705995, y1 = 1.9415348028467294, x2 = 0.37447514558298495, y2 = 1.9358100134206084)
 ComponentVector{Float64}(x1 = 0.5273708932464454, y1 = 1.9274169498739582, x2 = 0.553967302215615, y2 = 1.930951114134332)
 ComponentVector{Float64}(x1 = 0.809981593787924, y1 = 1.9937694275980937, x2 = 0.8755726816250095, y2 = 2.0136106002739864)
 ComponentVector{Float64}(x1 = 1.4047029365556867, y1 = 2.246788308065704, x2 = 1.4089431544995152, y2 = 2.266038277716898)
 ComponentVector{Float64}(x1 = 1.3703508712604588, y1 = 2.383745275264092, x2 = 1.3482558747129998, y2 = 2.393946908080564)
 ComponentVector{Float64}(x1 = 0.9866015800648426, y1 = 2.4754117031286023, x2 = 0.9657174373582494, y2 = 2.4755875710438926)
 ComponentVector{Float64}(x1 = 0.5633995688954301, y1 = 2.4063171030773507, x2 = 0.5476702671231246, y2 = 2.3999991379056844)
 ComponentVector{Float64}(x1 = 0.2592079151213729, y1 = 2.201596068957154, x2 = 0.2581560204327481, y2 = 2.195143760499976)
 ComponentVector{Float64}(x1 = 0.24693735366300354, y1 = 2.111217825368062, x2 = 0.24832158893536715, y2 = 2.10574018655429)

fig731 = Figure()
ax731 = Axis(fig731[1, 1];
    xlabel="Time (a.u.)",
    ylabel="Concentration (a.u.)",
    title="Fig. 7.31\nHasty Synthetic Oscillator Model",
)
lines!(ax731, 0..tend, t->sol731(t).x1, label="x1")
lines!(ax731, 0..tend, t->sol731(t).y1, label="y1")
lines!(ax731, 0..tend, t->sol731(t).x2, label="x2")
lines!(ax731, 0..tend, t->sol731(t).y2, label="y2")
axislegend(ax731; position=:rc)
fig731

This notebook was generated using Literate.jl.