Implement stiff solver

Author: BenLauwens

src/coroutines.jl

@@ -1,4 +1,4 @@
-type Coroutine <: DiscreteProcess
+mutable struct Coroutine <: DiscreteProcess
   bev :: BaseEvent
   fsm :: FiniteStateMachine
   target :: AbstractEvent

src/odes/QSS.jl

@@ -11,10 +11,10 @@ struct QSS{T} <: Integrator
   function QSS{T}(model::Model, t::Float64, x₀::Vector{Float64}, p::Vector{Float64};
                   order::Number=4, Δrel::Float64=1e-6, Δabs::Float64=1e-6) where T
     qss = new(UInt8(order), model, p, Vector{Taylor1{Float64}}(), Vector{Float64}(), Δrel, Δabs)
-    t₀ = t + Taylor1(Float64, order)
+    t₀ = t + Taylor1(Float64, order + 1)
     for q₀ in x₀
       push!(qss.t, t)
-      push!(qss.q, q₀ + Taylor1(zeros(Float64, order+1)))
+      push!(qss.q, q₀ + Taylor1(zeros(Float64, order + 1)))
     end
     for i in 1:order-1
       q = deepcopy(qss.q)
@@ -33,7 +33,7 @@ function Continuous(model::Model, env::Environment, x₀::Vector{Float64}, p::Ve
 end
 
 function initial_values(qss::QSS, t::Float64)
-  t₀ = t + Taylor1(Float64, qss.order+1)
+  t₀ = t + Taylor1(Float64, qss.order + 1)
   x₀ = Vector{Taylor1{Float64}}()
   for (i, f) in enumerate(qss.model.f)
     push!(x₀, integrate(f(t₀, qss.q, qss.p), qss.q[i][1]))
@@ -45,10 +45,10 @@ function step(var::Variable, cont::Continuous, qss::QSS)
   t = now(environment(var))
   n = length(qss.model.f)
   i = var.id
-  t₀ = t + Taylor1(Float64, qss.order+1)
+  t₀ = t + Taylor1(Float64, qss.order + 1)
   x₀ = advance_time(var, t)
   update_quantized_state(qss, var, t)
-  Δt = compute_next_time(var.x, max(qss.Δrel*x₀, qss.Δabs))
+  Δt = compute_next_time(qss, var)
   reset(var)
   schedule(var, Δt)
   for j in filter(j->qss.model.deps[j,i], 1:n)
@@ -60,7 +60,7 @@ function step(var::Variable, cont::Continuous, qss::QSS)
       advance_time(qss, k, t)
     end
     dep.x = integrate(qss.model.f[j](t₀, qss.q, qss.p), x₀)
-    Δt = recompute_next_time(qss, dep.x, qss.q[j], max(qss.Δrel*x₀, qss.Δabs))
+    Δt = recompute_next_time(qss, dep)
     reset(dep)
     schedule(dep, Δt)
   end
@@ -79,19 +79,82 @@ function update_quantized_state(qss::QSS{non_stiff}, var::Variable, t::Float64)
   qss.t[i] = t
 end
 
-function update_quantized_state(qss::QSS{stiff}, vars::Vector{Variable}, i::UInt, t::Float64)
+function update_quantized_state(qss::QSS{stiff}, var::Variable, t::Float64)
+  i = var.id
+  t₀ = t + Taylor1(Float64, qss.order + 1)
+  x₀ = evaluate(var.x)
+  Δq = max(qss.Δrel*x₀, qss.Δabs)
   for (j, istrue) in enumerate(qss.model.deps[i, :])
     istrue && advance_time(qss, j, t)
   end
-  q₋ = deepcopy(qss.q)
+  q̲  = deepcopy(qss.q)
+  q̲[i] = Taylor1(zeros(order+1))+x₀-Δq
+  x̲ = integrate(qss.model.f[i](t₀, q̲, qss.p), x₀)
+  for k in 1:order-1
+    q̲[i] = x̲-Δq
+    x̲ = integrate(qss.model.f[i](t₀, q̲, qss.p), x₀)
+  end
+  q̲[i] = x̲-Δq
+  q̲[i][end] = 0.0
+  q̅ = deepcopy(qss.q)
+  q̅[i] = Taylor1(zeros(order+1))+x₀+Δq
+  x̅ = integrate(qss.model.f[i](t₀, q̅, qss.p), x₀)
+  for k in 1:order-1
+    q̅[i] = x̅+Δq
+    x̅ = integrate(qss.model.f[i](t₀, q̅, qss.p), x₀)
+  end
+  q̅[i] = x̅+Δq
+  q̅[i][end] = 0.0
+  if x̲[end] * x̅[end] > 0.0
+    if x̅[end] > 0.0
+      var.x = deepcopy(x̅)
+      q = deepcopy(q̅)
+    else
+      var.x = deepcopy(x̲)
+      q = deepcopy(q̲)
+    end
+  else
+    q̃ = brent(nth_derivative, x₀-Δq, x₀+Δq, qss, i; xtol=min(Δq/100, 1e-7))
+    qss.q[i] = q̃ + Taylor1(zeros(Float64, order + 1))
+    var.x = integrate(qss.model.f[i](t₀, qss.q, qss.p), x₀)
+    for k in 1:order-1
+      qss.q[i] = deepcopy(var.x)
+      qss.q[i][1] = q̃
+      var.x = integrate(qss.model.f[i](t₀, qss.q, qss.p), x₀)
+    end
+    qss.q[i][end] = 0.0
+  end
 end
 
-function compute_next_time(x::Taylor1, Δq::Float64)
-  (abs(Δq/x[end]))^(1.0/x.order)
+function nth_derivative(q₀::Float64, qss::QSS{stiff}, i::UInt)
+  t₀ = t + Taylor1(Float64, qss.order + 1)
+  q = deepcopy(qss.q)
+  q[i][2:end] = 0.0
+  q[i][1] = q₀
+  q[i] = integrate(qss.model.f[i](t₀, q, qss.p), q₀)
+  for k in 1:order-1
+    q[i] = integrate(qss.model.f[i](t₀, q, qss.p), q₀)
+  end
+  q[i][end]
+end
+
+function compute_next_time(qss::QSS{non_stiff}, var::Variable)
+  x₀ = evaluate(var.x)
+  Δq = max(qss.Δrel*x₀, qss.Δabs)
+  (abs(Δq/var.x[end]))^(1.0/qss.order)
 end
 
-function recompute_next_time(::QSS{non_stiff}, x::Taylor1{Float64}, q::Taylor1{Float64}, Δq::Float64)
-  p = (x-q).coeffs
+function compute_next_time(qss::QSS{stiff}, var::Variable)
+  x₀ = evaluate(var.x)
+  Δq = max(qss.Δrel*x₀, qss.Δabs)
+  (abs(Δq/var.x[end]))^(1.0/qss.order)
+end
+
+function recompute_next_time(qss::QSS{non_stiff}, var::Variable)
+  i = var.id
+  x₀ = evaluate(var.x)
+  Δq = max(qss.Δrel*x₀, qss.Δabs)
+  p = (var.x-qss.q[i]).coeffs
   p[1] -= Δq
   neg = roots(p)
   p[1] += 2Δq

src/odes/roots.jl

@@ -137,9 +137,78 @@ function roots(coeff::Vector{Float64}) :: Vector{Complex{Float64}}
         mat[2:n-1, 1:n-2] = eye(Float64, n-2)
         mat[:, n-1] = - coeff[1:n-1] / coeff[n]
         res[1:n-1] = eigvals(mat)
-        #res[1:n-1] = roots(Poly(coeff))
       end
     end
   end
   res
 end
+
+function brent(f::Function, x0::Number, x1::Number, args...;
+               xtol::AbstractFloat=1e-7, ytol=2eps(Float64),
+               maxiter::Integer=50)
+  EPS = eps(Float64)
+  y0 = f(x0,args...)
+  y1 = f(x1,args...)
+  if abs(y0) < abs(y1)
+      # Swap lower and upper bounds.
+      x0, x1 = x1, x0
+      y0, y1 = y1, y0
+  end
+  x2 = x0
+  y2 = y0
+  x3 = x2
+  bisection = true
+  for _ in 1:maxiter
+    # x-tolerance.
+    if abs(x1-x0) < xtol
+      return x1
+    end
+
+    # Use inverse quadratic interpolation if f(x0)!=f(x1)!=f(x2)
+    # and linear interpolation (secant method) otherwise.
+    if abs(y0-y2) > ytol && abs(y1-y2) > ytol
+      x = x0*y1*y2/((y0-y1)*(y0-y2)) +
+          x1*y0*y2/((y1-y0)*(y1-y2)) +
+          x2*y0*y1/((y2-y0)*(y2-y1))
+    else
+      x = x1 - y1 * (x1-x0)/(y1-y0)
+    end
+
+    # Use bisection method if satisfies the conditions.
+    delta = abs(2EPS*abs(x1))
+    min1 = abs(x-x1)
+    min2 = abs(x1-x2)
+    min3 = abs(x2-x3)
+    if (x < (3x0+x1)/4 && x > x1) ||
+       (bisection && min1 >= min2/2) ||
+       (!bisection && min1 >= min3/2) ||
+       (bisection && min2 < delta) ||
+       (!bisection && min3 < delta)
+      x = (x0+x1)/2
+      bisection = true
+    else
+      bisection = false
+    end
+
+    y = f(x,args...)
+    # y-tolerance.
+    if abs(y) < ytol
+      return x
+    end
+    x3 = x2
+    x2 = x1
+    if sign(y0) != sign(y)
+      x1 = x
+      y1 = y
+    else
+      x0 = x
+      y0 = y
+    end
+    if abs(y0) < abs(y1)
+        # Swap lower and upper bounds.
+      x0, x1 = x1, x0
+      y0, y1 = y1, y0
+    end
+  end
+  error("Max iteration exceeded")
+end

src/processes.jl

@@ -1,4 +1,4 @@
-type Process <: DiscreteProcess
+mutable struct Process <: DiscreteProcess
   bev :: BaseEvent
   task :: Task
   target :: AbstractEvent

test/continuous.jl

@@ -27,7 +27,7 @@ end
 sim = Simulation()
 cont = @continuous diffeq(sim, [0.0, 20.0], [2020.0, 0.0]; stiff=false, order=4, Δrel=1e-16, Δabs=1e-6)
 @process report(sim, cont)
-@time run(sim, 500)
+@time run(sim, 2000)
 
 @model function bouncing_ball(t, x, p, dx)
   g = 9.81

test/events.jl

@@ -1,6 +1,6 @@
 using SimJulia
 
-type TestException <: Exception end
+struct TestException <: Exception end
 
 function test_callback_event(ev::Event)
   println("Hi $ev has value $(value(ev))")

test/simulations.jl

@@ -1,6 +1,6 @@
 using SimJulia
 
-type TestException <: Exception end
+struct TestException <: Exception end
 
 function test_callback(ev::AbstractEvent)
   println("Hi I timed out at $(now(environment(ev)))")