using PyPlot
matplotlib[:rc]("text", usetex = true)
matplotlib[:rc]("figure", figsize = (4, 3), dpi = 300)
INFO: Recompiling stale cache file /home/jgraw/.julia/lib/v0.6/MacroTools.ji for module MacroTools. INFO: Recompiling stale cache file /home/jgraw/.julia/lib/v0.6/PyPlot.ji for module PyPlot.
function rhs(x, h, t, τs, τc, σ = tanh)
return -x - σ(h(t - τs) - h(t - τs - τc))
end
rhs (generic function with 2 methods)
using ApproxFun
using NLsolve
INFO: Recompiling stale cache file /home/jgraw/.julia/lib/v0.6/ApproxFun.ji for module ApproxFun. INFO: Recompiling stale cache file /home/jgraw/.julia/lib/v0.6/NLsolve.ji for module NLsolve.
function residual!(res, xs, basis, τs, τc, σ)
f = Fun(basis, xs)
df = f'
resid = t -> df(t) - rhs(f(t), f, t, τs, τc, σ)
res .= coefficients(Fun(resid, basis, length(xs)))
return nothing
end
residual! (generic function with 1 method)
τs = 0.2
τc = 1.0
κ = 5.0
σ = x -> tanh(κ * x)
(::#3) (generic function with 1 method)
residuals = Array{Float64}(0)
amplitudes = similar(residuals)
n = 20
periods = linspace(0.6, 1.6, 2000)
for period in periods
s = Fourier(0..period)
initial_f = Fun(x -> 0.25sin(2pi * x / period), s, n)
res = nlsolve(
(res, xs) -> residual!(res, xs, s, τs, τc, σ),
coefficients(initial_f),
iterations = 10)
push!(residuals, res.residual_norm)
push!(amplitudes, maximum(abs.(res.zero)))
end
plot(periods, residuals, label="residual")
plot(periods, amplitudes, label="amplitude")
yscale("log")
ylim(ymin=1e-8)
legend()
grid()
xlabel("periodicity")
PyObject Text(0.5,72,'periodicity')
# signal-to-noise ratio (significance)
plot(periods, amplitudes ./ residuals, label="amplitude / residual", "-")
yscale("log")
#ylim(1e-2, 1e2)
#xlim(0.4, 1.6)
legend()
grid()
xlabel("periodicity")
axhline(1, color="black", ls="--")
PyObject <matplotlib.lines.Line2D object at 0x7febfe612ef0>
$\Rightarrow$ the two limit cycles coexist at $\tau_c=1$
using Optim
INFO: Recompiling stale cache file /home/jgraw/.julia/lib/v0.6/Optim.ji for module Optim.
function prominence(period, τs, τc, σ, n=20)
s = Fourier(0..period)
initial_f = Fun(x -> 0.25sin(2pi * x / period), s, n)
res = nlsolve(
(res, xs) -> residual!(res, xs, s, τs, τc, σ),
coefficients(initial_f),
iterations = 10)
return maximum(abs.(res.zero)) / res.residual_norm
end
prominence (generic function with 2 methods)
# example
res = optimize(0.5, 0.9, GoldenSection()) do period
return -prominence(period, τs, τc, σ, 20)
end
period = Optim.minimizer(res)
display(res)
Results of Optimization Algorithm * Algorithm: Golden Section Search * Search Interval: [0.500000, 0.900000] * Minimizer: 6.880456e-01 * Minimum: -2.581134e+05 * Iterations: 34 * Convergence: max(|x - x_upper|, |x - x_lower|) <= 2*(1.5e-08*|x|+2.2e-16): true * Objective Function Calls: 35
residuals = Array{Float64}(0)
amplitudes = similar(residuals)
n = 20
periods = linspace(1.3, 1.35, 1_000)
for period in periods
s = Fourier(0..period)
initial_f = Fun(x -> 0.25sin(2pi * x / period), s, n)
res = nlsolve(
(res, xs) -> residual!(res, xs, s, τs, τc, σ),
coefficients(initial_f),
iterations = 10)
push!(residuals, res.residual_norm)
push!(amplitudes, maximum(abs.(res.zero)))
end
plot(periods, residuals, label="residual", ".")
plot(periods, amplitudes, label="amplitude", ".")
yscale("log")
legend()
grid()
xlabel("periodicity")
periods[indmin(residuals)] # approximate location of limit cycle
1.3176176176176175
# example
res = optimize(1.2, 1.4, GoldenSection()) do period
return -prominence(period, τs, τc, σ, 20)
end
period = Optim.minimizer(res)
display(res)
Results of Optimization Algorithm * Algorithm: Golden Section Search * Search Interval: [1.200000, 1.400000] * Minimizer: 1.317599e+00 * Minimum: -4.455105e+06 * Iterations: 32 * Convergence: max(|x - x_upper|, |x - x_lower|) <= 2*(1.5e-08*|x|+2.2e-16): true * Objective Function Calls: 33
n = 20
#period = 1.3176176176176175
s = Fourier(0..period)
initial_f = Fun(x -> sin(2pi * x / period), s, n)
res = nlsolve(
(res, xs) -> residual!(res, xs, s, τs, τc, σ),
coefficients(initial_f),
ftol = 1e-6)
u1 = Fun(s, res.zero)
display(res)
Results of Nonlinear Solver Algorithm * Algorithm: Trust-region with dogleg and autoscaling * Starting Point: [-3.55101e-17, 1.0, -5.02205e-17, 6.19944e-17, -1.38075e-16, -4.00619e-17, -1.34016e-16, 1.16442e-16, -3.83097e-17, 2.22045e-17, -3.44509e-17, 1.00309e-16, 5.43213e-17, 1.2888e-16, -3.50143e-18, 1.30333e-16, 1.04436e-16, 8.88178e-17, 1.60956e-16, -1.1437e-16] * Zero: [3.52262e-16, 0.201403, -0.00520023, -2.99654e-16, -9.14944e-17, -0.00678723, 0.00654776, 8.85457e-17, -1.34121e-17, 0.000283827, -0.000964783, 1.94363e-17, 1.17696e-17, 2.3924e-5, 0.000124478, -1.44738e-17, 1.10596e-17, -8.02232e-6, -1.52482e-5, 1.22766e-18] * Inf-norm of residuals: 0.000000 * Iterations: 4 * Convergence: true * |x - x'| < 0.0e+00: false * |f(x)| < 1.0e-06: true * Function Calls (f): 5 * Jacobian Calls (df/dx): 5
t2 = linspace(0, 5, 1000)
plot(t2, u1.(t2))
1-element Array{PyCall.PyObject,1}: PyObject <matplotlib.lines.Line2D object at 0x7febfd45cf28>
residuals = Array{Float64}(0)
amplitudes = similar(residuals)
n = 20
periods = linspace(0.6875, 0.6885, 1_000)
for period in periods
s = Fourier(0..period)
initial_f = Fun(x -> 0.25sin(2pi * x / period), s, n)
res = nlsolve(
(res, xs) -> residual!(res, xs, s, τs, τc, σ),
coefficients(initial_f),
iterations = 10)
push!(residuals, res.residual_norm)
push!(amplitudes, maximum(abs.(res.zero)))
end
plot(periods, residuals, label="residual", ".")
plot(periods, amplitudes, label="amplitude", ".")
yscale("log")
legend()
grid()
xlabel("periodicity")
periods[indmin(residuals)]
0.6880455455455455
# example
res = optimize(0.6, 0.8, GoldenSection()) do period
return -prominence(period, τs, τc, σ, 20)
end
period = Optim.minimizer(res)
display(res)
Results of Optimization Algorithm * Algorithm: Golden Section Search * Search Interval: [0.600000, 0.800000] * Minimizer: 6.880456e-01 * Minimum: -1.616248e+06 * Iterations: 33 * Convergence: max(|x - x_upper|, |x - x_lower|) <= 2*(1.5e-08*|x|+2.2e-16): true * Objective Function Calls: 34
n = 20
#period = 0.6880455455455455
s = Fourier(0..period)
initial_f = Fun(x -> 0.25sin(2pi * x / period), s, n)
res = nlsolve(
(res, xs) -> residual!(res, xs, s, τs, τc, σ),
coefficients(initial_f))
display(res)
t2 = linspace(0, 5, 1000)
u2 = Fun(s, res.zero)
Results of Nonlinear Solver Algorithm * Algorithm: Trust-region with dogleg and autoscaling * Starting Point: [-1.24475e-18, 0.25, -8.08673e-17, 4.92372e-17, -1.93931e-17, 3.33067e-17, 1.02233e-17, 3.41692e-17, -1.18889e-19, 1.66533e-17, 3.57962e-17, 2.06417e-17, 1.16386e-17, 3.33067e-17, 4.78372e-17, -8.06624e-18, -1.87379e-17, 6.66134e-17, 2.72138e-17, 1.2347e-17] * Zero: [-5.4121e-17, 0.0566217, 0.000303099, 2.41751e-18, -2.29399e-18, -0.000352873, 7.54974e-5, -3.24691e-18, -1.74099e-19, 4.10706e-6, -1.25653e-6, 9.28576e-19, 4.07465e-19, -5.69598e-8, 2.22776e-8, -2.80564e-18, -1.21155e-18, 8.35589e-10, -4.17137e-10, 1.23465e-18] * Inf-norm of residuals: 0.000000 * Iterations: 23 * Convergence: true * |x - x'| < 0.0e+00: false * |f(x)| < 1.0e-08: true * Function Calls (f): 24 * Jacobian Calls (df/dx): 17
Fun(Fourier(【0.0,0.6880456129134734❫),[-5.4121e-17, 0.0566217, 0.000303099, 2.41751e-18, -2.29399e-18, -0.000352873, 7.54974e-5, -3.24691e-18, -1.74099e-19, 4.10706e-6, -1.25653e-6, 9.28576e-19, 4.07465e-19, -5.69598e-8, 2.22776e-8, -2.80564e-18, -1.21155e-18, 8.35589e-10, -4.17137e-10, 1.23465e-18])
plot(t2, u1.(t2))
plot(t2, u2.(t2), "--")
#ylim(-0.22, 0.22)
xlabel(L"time $t$")
ylabel(L"process $A(t)$")
tight_layout()
using ApproxFun
using NLsolve
using Optim
using OrdinaryDiffEq
using DelayDiffEq
INFO: Recompiling stale cache file /home/jgraw/.julia/lib/v0.6/RecursiveArrayTools.ji for module RecursiveArrayTools. INFO: Recompiling stale cache file /home/jgraw/.julia/lib/v0.6/OrdinaryDiffEq.ji for module OrdinaryDiffEq. INFO: Recompiling stale cache file /home/jgraw/.julia/lib/v0.6/DelayDiffEq.ji for module DelayDiffEq.
function residual!(res, xs, basis, τs, τc, σ)
f = Fun(basis, xs)
df = f'
resid = t -> df(t) - rhs(f(t), f, t, τs, τc, σ)
res .= coefficients(Fun(resid, basis, length(xs)))
return nothing
end
residual! (generic function with 1 method)
τs = 0.2
τc = 1.0 # compare 1.0
κ = 5.0
σ = x -> tanh(κ * x)
(::#37) (generic function with 1 method)
# example
res = optimize(0.6, 0.9, GoldenSection()) do period
return -prominence(period, τs, τc, σ, 20)
end
period = Optim.minimizer(res)
display(res)
Results of Optimization Algorithm * Algorithm: Golden Section Search * Search Interval: [0.600000, 0.900000] * Minimizer: 6.880456e-01 * Minimum: -1.616247e+06 * Iterations: 34 * Convergence: max(|x - x_upper|, |x - x_lower|) <= 2*(1.5e-08*|x|+2.2e-16): true * Objective Function Calls: 35
#period = 0.6880456129134734
n = 20
s = Fourier(0..period)
initial_f = Fun(x -> 0.25sin(2pi * x / period), s, n)
res = nlsolve(
(res, xs) -> residual!(res, xs, s, τs, τc, σ),
coefficients(initial_f))
u2 = Fun(s, res.zero)
display(res)
Results of Nonlinear Solver Algorithm * Algorithm: Trust-region with dogleg and autoscaling * Starting Point: [4.17764e-17, 0.25, 5.48067e-18, -1.58383e-17, 1.08582e-17, -5.33376e-17, -5.07818e-18, 1.2463e-18, -3.09692e-18, -2.22045e-17, -3.06162e-18, 9.78859e-18, 7.75504e-18, 3.11332e-17, -1.04506e-18, -1.38217e-17, -3.10254e-17, 4.44089e-17, -1.16039e-17, 1.09592e-17] * Zero: [-4.23347e-17, 0.0566217, 0.000303099, 1.73737e-20, 4.23792e-18, -0.000352873, 7.54974e-5, -3.37597e-19, -8.14445e-19, 4.10706e-6, -1.25653e-6, 3.64785e-19, 4.13686e-19, -5.69598e-8, 2.22776e-8, -7.81509e-19, -5.3754e-19, 8.35589e-10, -4.17137e-10, -1.31471e-19] * Inf-norm of residuals: 0.000000 * Iterations: 23 * Convergence: true * |x - x'| < 0.0e+00: false * |f(x)| < 1.0e-08: true * Function Calls (f): 24 * Jacobian Calls (df/dx): 17
prob = DDEProblem(
(u, h, p, t) -> rhs(u, t -> h(p, t), t, τs, τc, σ),
u2(0.), (p, t) -> u2(t), (0., 100.),
constant_lags = [τs, τs + τc])
sol = solve(prob, MethodOfSteps(Tsit5()));
figure(figsize = (12, 3))
t2 = linspace(0., 100., 10_000)
plot(t2, sol.(t2))
t1 = linspace(-10., 0., 1_000)
plot(t1, u2.(t1))
xlabel(L"time $t$")
ylabel(L"process $A(t)$")
xlim(-10, 100)
axvline(0, color="black", ls="--")
PyObject <matplotlib.lines.Line2D object at 0x7febfc92aba8>
$\Rightarrow$ the fast limit cycle is unstable for $\tau_c=1.0$ and stable for $\tau_c=1.2$
# example
res = optimize(1.0, 1.7, GoldenSection()) do period
return -prominence(period, τs, τc, σ, 20)
end
period = Optim.minimizer(res)
display(res)
Results of Optimization Algorithm * Algorithm: Golden Section Search * Search Interval: [1.000000, 1.700000] * Minimizer: 1.317599e+00 * Minimum: -5.570785e+06 * Iterations: 34 * Convergence: max(|x - x_upper|, |x - x_lower|) <= 2*(1.5e-08*|x|+2.2e-16): true * Objective Function Calls: 35
n = 20
#period = 1.317599
s = Fourier(0..period)
initial_f = Fun(x -> 0.25sin(2pi * x / period), s, n)
res = nlsolve(
(res, xs) -> residual!(res, xs, s, τs, τc, σ),
coefficients(initial_f),
iterations = 10,
ftol = 1e-5)
u1 = Fun(s, res.zero)
display(res)
Results of Nonlinear Solver Algorithm * Algorithm: Trust-region with dogleg and autoscaling * Starting Point: [3.61248e-18, 0.25, -2.75624e-17, -6.09514e-18, -7.69823e-18, -2.16611e-17, -1.37719e-17, 5.20587e-18, -8.06451e-18, 5.55112e-18, 8.04061e-18, 1.94301e-18, 1.86787e-17, -5.43382e-19, 2.69229e-17, -1.13746e-17, 1.21061e-17, -0.0, -8.9373e-18, 3.32642e-18] * Zero: [1.63883e-16, 0.201453, -0.00263917, -5.67875e-17, 6.30805e-17, -0.00703199, 0.00628418, 1.19685e-17, -5.65843e-17, 0.000344566, -0.00094484, -3.76048e-18, 2.70913e-18, 1.27718e-5, 0.000126291, -2.11326e-18, 6.37059e-18, -6.62078e-6, -1.67743e-5, 2.53448e-18] * Inf-norm of residuals: 0.000000 * Iterations: 3 * Convergence: true * |x - x'| < 0.0e+00: false * |f(x)| < 1.0e-05: true * Function Calls (f): 4 * Jacobian Calls (df/dx): 4
prob = DDEProblem(
(u, h, p, t) -> rhs(u, t -> h(p, t), t, τs, τc, σ),
u1(0.), (p, t) -> u1(t), (0., 100.),
constant_lags = [τs, τs + τc])
sol = solve(prob, MethodOfSteps(Tsit5()));
figure(figsize = (12, 3))
t2 = linspace(0., 100., 10_000)
plot(t2, sol.(t2))
t1 = linspace(-10., 0., 1_000)
plot(t1, u1.(t1))
xlabel(L"time $t$")
ylabel(L"process $A(t)$")
xlim(-10, 100)
axvline(0, color="black", ls="--")
PyObject <matplotlib.lines.Line2D object at 0x7febfc8ec518>
$\Rightarrow$ the slow limit cycle is stable from $\tau_c=1.0$ to $\tau_c=1.3$