From 14635192a26f1af1ab94c3d4c12799e139b7a6d0 Mon Sep 17 00:00:00 2001 From: Chris Rackauckas Date: Thu, 23 Jun 2022 15:35:36 -0400 Subject: [PATCH] format SciML Style --- .JuliaFormatter.toml | 1 + .github/workflows/FormatCheck.yml | 42 ++++ examples/User_Type.jl | 18 +- src/ODE.jl | 244 ++++++++++---------- src/common.jl | 85 +++---- src/runge_kutta.jl | 370 +++++++++++++++--------------- test/common.jl | 16 +- test/interface-tests.jl | 61 ++--- test/runtests.jl | 93 ++++---- 9 files changed, 497 insertions(+), 433 deletions(-) create mode 100644 .JuliaFormatter.toml create mode 100644 .github/workflows/FormatCheck.yml diff --git a/.JuliaFormatter.toml b/.JuliaFormatter.toml new file mode 100644 index 000000000..453925c3f --- /dev/null +++ b/.JuliaFormatter.toml @@ -0,0 +1 @@ +style = "sciml" \ No newline at end of file diff --git a/.github/workflows/FormatCheck.yml b/.github/workflows/FormatCheck.yml new file mode 100644 index 000000000..2a3517a0f --- /dev/null +++ b/.github/workflows/FormatCheck.yml @@ -0,0 +1,42 @@ +name: format-check + +on: + push: + branches: + - 'master' + - 'release-' + tags: '*' + pull_request: + +jobs: + build: + runs-on: ${{ matrix.os }} + strategy: + matrix: + julia-version: [1] + julia-arch: [x86] + os: [ubuntu-latest] + steps: + - uses: julia-actions/setup-julia@latest + with: + version: ${{ matrix.julia-version }} + + - uses: actions/checkout@v1 + - name: Install JuliaFormatter and format + # This will use the latest version by default but you can set the version like so: + # + # julia -e 'using Pkg; Pkg.add(PackageSpec(name="JuliaFormatter", version="0.13.0"))' + run: | + julia -e 'using Pkg; Pkg.add(PackageSpec(name="JuliaFormatter"))' + julia -e 'using JuliaFormatter; format(".", verbose=true)' + - name: Format check + run: | + julia -e ' + out = Cmd(`git diff --name-only`) |> read |> String + if out == "" + exit(0) + else + @error "Some files have not been formatted !!!" + write(stdout, out) + exit(1) + end' diff --git a/examples/User_Type.jl b/examples/User_Type.jl index 427eb4c09..6bac43a21 100644 --- a/examples/User_Type.jl +++ b/examples/User_Type.jl @@ -17,16 +17,16 @@ using Compat const DataFloat = Float64 mutable struct Vec2 - a::DataFloat - b::DataFloat + a::DataFloat + b::DataFloat end # This is the minimal set of function required on the type inorder to work with # the ODE module -Base.:/(x::Vec2, y::Real) = Vec2(x.a/y, x.b/y) +Base.:/(x::Vec2, y::Real) = Vec2(x.a / y, x.b / y) Base.:*(y::Real, x::Vec2) = Vec2(y * x.a, y * x.b) -Base.:*(x::Vec2, y::Real) = y*x -Base.:.*(y::Real, x::Vec2) = y*x +Base.:*(x::Vec2, y::Real) = y * x +Base.:.*(y::Real, x::Vec2) = y * x Base.:+(x::Vec2, y::Vec2) = Vec2(x.a + y.a, x.b + y.b) Base.:-(x::Vec2, y::Vec2) = Vec2(x.a - y.a, x.b - y.b) Base.norm(x::Vec2) = sqrt(x.a^2 + x.b^2) @@ -34,13 +34,13 @@ Base.zero(x::Type{Vec2}) = Vec2(zero(DataFloat), zero(DataFloat)) ODE.isoutofdomain(x::Vec2) = isnan(x.a) || isnan(x.b) # RHS function -f(t,y) = Vec2(-y.b, y.a) +f(t, y) = Vec2(-y.b, y.a) # Initial condtions start = Vec2(0.0, 0.1) # Time vector going from 0 to 2*PI in 0.01 increments -time = 0:0.1:4*pi +time = 0:0.1:(4 * pi) # Solve the ODE t, y = ode45(f, start, time) @@ -48,6 +48,6 @@ t, y = ode45(f, start, time) # Plot the solution a = map(y -> y.a, y) b = map(y -> y.b, y) -plot(t, a, label="a(t)") -plot(t, b, label="b(t)") +plot(t, a, label = "a(t)") +plot(t, b, label = "b(t)") legend() diff --git a/src/ODE.jl b/src/ODE.jl index a590085e2..8ea1f8728 100644 --- a/src/ODE.jl +++ b/src/ODE.jl @@ -11,10 +11,9 @@ using Reexport import DiffEqBase: solve -const warnkeywords = - (:save_idxs, :d_discontinuities, :isoutofdomain, :unstable_check, - :calck, :progress, - :internalnorm, :gamma, :beta1, :beta2, :qmax, :qmin, :qoldinit) +const warnkeywords = (:save_idxs, :d_discontinuities, :isoutofdomain, :unstable_check, + :calck, :progress, + :internalnorm, :gamma, :beta1, :beta2, :qmax, :qmin, :qoldinit) function __init__() global warnlist = Set(warnkeywords) @@ -44,7 +43,7 @@ export ODEjlAlgorithm # Butcher Tableaus, or more generally coefficient tables # see Hairer & Wanner 1992, p. 134, 166 -abstract type Tableau{Name, S, T<:Real} end +abstract type Tableau{Name, S, T <: Real} end # Name is the name of the tableau/method (a symbol) # S is the number of stages (an int) # T is the type of the coefficients @@ -69,12 +68,13 @@ abstract type Tableau{Name, S, T<:Real} end # | b_1 ... b_s this is the one used for stepping # | b'_1 ... b'_s this is the one used for error-checking -Base.eltype(b::Tableau{N,S,T}) where {N,S,T} = T +Base.eltype(b::Tableau{N, S, T}) where {N, S, T} = T order(b::Tableau) = b.order # Subtypes need to define a convert method to convert to a different # eltype with signature: -Base.convert(::Type{Tnew}, tab::Tableau) where {Tnew<:Real} = +function Base.convert(::Type{Tnew}, tab::Tableau) where {Tnew <: Real} error("Define convert method for concrete Tableau types") +end ############################################################################### ## HELPER FUNCTIONS @@ -82,33 +82,33 @@ Base.convert(::Type{Tnew}, tab::Tableau) where {Tnew<:Real} = # estimator for initial step based on book # "Solving Ordinary Differential Equations I" by Hairer et al., p.169 -function hinit(F, x0, t0::T, tend, p, reltol, abstol) where T +function hinit(F, x0, t0::T, tend, p, reltol, abstol) where {T} # Returns first step, direction of integration and F evaluated at t0 - tdir = sign(tend-t0) - tdir==0 && error("Zero time span") - tau = max(reltol*norm(x0, Inf), abstol) - d0 = norm(x0, Inf)/tau + tdir = sign(tend - t0) + tdir == 0 && error("Zero time span") + tau = max(reltol * norm(x0, Inf), abstol) + d0 = norm(x0, Inf) / tau f0 = F(t0, x0) - d1 = norm(f0, Inf)/tau + d1 = norm(f0, Inf) / tau if d0 < 1e-5 || d1 < 1e-5 h0 = 1e-6 else - h0 = 0.01*(d0/d1) + h0 = 0.01 * (d0 / d1) end - h0 = convert(T,h0) + h0 = convert(T, h0) # perform Euler step - x1 = x0 + tdir*h0*f0 - f1 = F(t0 + tdir*h0, x1) + x1 = x0 + tdir * h0 * f0 + f1 = F(t0 + tdir * h0, x1) # estimate second derivative - d2 = norm(f1 - f0, Inf)/(tau*h0) + d2 = norm(f1 - f0, Inf) / (tau * h0) if max(d1, d2) <= 1e-15 - h1 = max(T(10)^(-6), T(10)^(-3)*h0) + h1 = max(T(10)^(-6), T(10)^(-3) * h0) else - pow = -(2 + log10(max(d1, d2)))/(p + 1) + pow = -(2 + log10(max(d1, d2))) / (p + 1) h1 = 10^pow end - h1 = convert(T,h1) - return tdir*min(100*h0, h1, tdir*(tend-t0)), tdir, f0 + h1 = convert(T, h1) + return tdir * min(100 * h0, h1, tdir * (tend - t0)), tdir, f0 end # isoutofdomain takes the state and returns true if state is outside @@ -143,22 +143,24 @@ function make_consistent_types(fn, y0, tspan, btab::Tableau) # On time container: eltype Ty = typeof(y0) - Eyf = typeof(y0[1]/(tspan[end]-tspan[1])) + Eyf = typeof(y0[1] / (tspan[end] - tspan[1])) Et = eltype(tspan) - @assert Et<:Real - if !(Et<:AbstractFloat) + @assert Et <: Real + if !(Et <: AbstractFloat) Et = promote_type(Et, Float64) end # if all are Floats, make them the same - if Et<:AbstractFloat && Eyf<:AbstractFloat + if Et <: AbstractFloat && Eyf <: AbstractFloat Et = promote_type(Et, Eyf) Eyf = Et end - !isconcretetype(Et) && @warn("The eltype(tspan) is not a concrete type! Change type of tspan for better performance.") - !isconcretetype(Eyf) && @warn("The eltype(y0/tspan[1]) is not a concrete type! Change type of y0 and/or tspan for better performance.") + !isconcretetype(Et) && + @warn("The eltype(tspan) is not a concrete type! Change type of tspan for better performance.") + !isconcretetype(Eyf) && + @warn("The eltype(y0/tspan[1]) is not a concrete type! Change type of y0 and/or tspan for better performance.") btab_ = convert(Et, btab) return Et, Eyf, Ty, btab_ @@ -174,7 +176,7 @@ include("runge_kutta.jl") # with Adams-Bashforth-Moulton coefficients function ode_ms(F, x0, tspan, order::Integer; kwargs...) h = diff(tspan) - x = Vector{typeof(x0)}(undef,length(tspan)) + x = Vector{typeof(x0)}(undef, length(tspan)) x[1] = x0 if 1 <= order <= 4 @@ -182,27 +184,28 @@ function ode_ms(F, x0, tspan, order::Integer; kwargs...) else b = zeros(order, order) b[1:4, 1:4] = ms_coefficients4 - for s = 5:order - for j = 0:(s - 1) + for s in 5:order + for j in 0:(s - 1) # Assign in correct order for multiplication below # (a factor depending on j and s) .* (an integral of a polynomial with -(0:s), except -j, as roots) - p_int = polyint(poly(diagm(0 => -[0:j - 1; j + 1:s - 1]))) + p_int = polyint(poly(diagm(0 => -[0:(j - 1); (j + 1):(s - 1)]))) b[s, s - j] = ((-1)^j / factorial(j) - / factorial(s - 1 - j) * polyval(p_int, 1)) + / + factorial(s - 1 - j) * polyval(p_int, 1)) end end end # TODO: use a better data structure here (should be an order-element circ buffer) xdot = similar(x) - for i = 1:length(tspan)-1 + for i in 1:(length(tspan) - 1) # Need to run the first several steps at reduced order steporder = min(i, order) xdot[i] = F(tspan[i], x[i]) - x[i+1] = x[i] - for j = 1:steporder - x[i+1] += h[i]*b[steporder, j]*xdot[i-(steporder-1) + (j-1)] + x[i + 1] = x[i] + for j in 1:steporder + x[i + 1] += h[i] * b[steporder, j] * xdot[i - (steporder - 1) + (j - 1)] end end return vcat(tspan), x @@ -223,21 +226,21 @@ function fdjacobian(F, x::Number, t) ftx = F(t, x) # The 100 below is heuristic - dx = (x .+ (x==0))./100 - dFdx = (F(t,x+dx)-ftx)./dx + dx = (x .+ (x == 0)) ./ 100 + dFdx = (F(t, x + dx) - ftx) ./ dx return dFdx end function fdjacobian(F, x, t) ftx = F(t, x) - lx = max(length(x),1) + lx = max(length(x), 1) dFdx = zeros(eltype(x), lx, lx) - for j = 1:lx + for j in 1:lx # The 100 below is heuristic dx = zeros(eltype(x), lx) - dx[j] = (x[j] .+ (x[j]==0))./100 - dFdx[:,j] = (F(t,x+dx)-ftx)./dx[j] + dx[j] = (x[j] .+ (x[j] == 0)) ./ 100 + dFdx[:, j] = (F(t, x + dx) - ftx) ./ dx[j] end return dFdx end @@ -251,38 +254,37 @@ end # jacobian = G(t,y)::Function | nothing (FD) function ode23s(F, y0, tspan; reltol = 1.0e-5, abstol = 1.0e-8, - jacobian=nothing, - points=:all, - norm=LinearAlgebra.norm, - minstep=abs(tspan[end] - tspan[1])/1e18, - maxstep=abs(tspan[end] - tspan[1])/2.5, - initstep=0.) + jacobian = nothing, + points = :all, + norm = LinearAlgebra.norm, + minstep = abs(tspan[end] - tspan[1]) / 1e18, + maxstep = abs(tspan[end] - tspan[1]) / 2.5, + initstep = 0.0) # select method for computing the Jacobian if typeof(jacobian) == Function jac = jacobian else # fallback finite-difference - jac = (t, y)->fdjacobian(F, y, t) + jac = (t, y) -> fdjacobian(F, y, t) end # constants - d = 1/(2 + sqrt(2)) + d = 1 / (2 + sqrt(2)) e32 = 6 + sqrt(2) - # initialization t = tspan[1] tfinal = tspan[end] h = initstep - if h == 0. + if h == 0.0 # initial guess at a step size h, tdir, F0 = hinit(F, y0, t, tfinal, 3, reltol, abstol) else tdir = sign(tfinal - t) - F0 = F(t,y0) + F0 = F(t, y0) end h = tdir * min(abs(h), maxstep) @@ -290,55 +292,56 @@ function ode23s(F, y0, tspan; tout = [t] # first output time yout = [deepcopy(y)] # first output solution - - J = jac(t,y) # get Jacobian of F wrt y + J = jac(t, y) # get Jacobian of F wrt y while abs(t - tfinal) > 0 && minstep < abs(h) - if abs(t-tfinal) < abs(h) + if abs(t - tfinal) < abs(h) h = tfinal - t end - if size(J,1) == 1 - W = I - h*d*J + if size(J, 1) == 1 + W = I - h * d * J else # note: if there is a mass matrix M on the lhs of the ODE, i.e., # M * dy/dt = F(t,y) # we can simply replace eye(J) by M in the following expression # (see Sec. 5 in [SR97]) - W = lu( I - h*d*J ) + W = lu(I - h * d * J) end # approximate time-derivative of F - T = h*d*(F(t + h/100, y) - F0)/(h/100) + T = h * d * (F(t + h / 100, y) - F0) / (h / 100) # modified Rosenbrock formula - k1 = W\(F0 + T) - F1 = F(t + 0.5*h, y + 0.5*h*k1) - k2 = W\(F1 - k1) + k1 - ynew = y + h*k2 + k1 = W \ (F0 + T) + F1 = F(t + 0.5 * h, y + 0.5 * h * k1) + k2 = W \ (F1 - k1) + k1 + ynew = y + h * k2 F2 = F(t + h, ynew) - k3 = W\(F2 - e32*(k2 - F1) - 2*(k1 - F0) + T ) + k3 = W \ (F2 - e32 * (k2 - F1) - 2 * (k1 - F0) + T) - err = (abs(h)/6)*norm(k1 - 2*k2 + k3) # error estimate - delta = max(reltol*max(norm(y),norm(ynew)), abstol) # allowable error + err = (abs(h) / 6) * norm(k1 - 2 * k2 + k3) # error estimate + delta = max(reltol * max(norm(y), norm(ynew)), abstol) # allowable error # check if new solution is acceptable - if err <= delta - - if points==:specified || points==:all + if err <= delta + if points == :specified || points == :all # only points in tspan are requested # -> find relevant points in (t,t+h] - for toi in tspan[(tspan.>t) .& (tspan.<=t+h)] + for toi in tspan[(tspan .> t) .& (tspan .<= t + h)] # rescale to (0,1] - s = (toi-t)/h + s = (toi - t) / h # use interpolation formula to get solutions at t=toi push!(tout, toi) - push!(yout, y + h*( k1*s*(1-s)/(1-2*d) + k2*s*(s-2*d)/(1-2*d))) + push!(yout, + y + + h * (k1 * s * (1 - s) / (1 - 2 * d) + + k2 * s * (s - 2 * d) / (1 - 2 * d))) end end - if (points==:all) && (tout[end]!=t+h) + if (points == :all) && (tout[end] != t + h) # add the intermediate points push!(tout, t + h) push!(yout, ynew) @@ -349,30 +352,28 @@ function ode23s(F, y0, tspan; y = ynew F0 = F2 # use FSAL property - J = jac(t,y) # get Jacobian of F wrt y - # for new solution + J = jac(t, y) # get Jacobian of F wrt y + # for new solution end # update of the step size - h = tdir*min( maxstep, abs(h)*0.8*(delta/err)^(1/3) ) + h = tdir * min(maxstep, abs(h) * 0.8 * (delta / err)^(1 / 3)) end return tout, yout end - #ODEROSENBROCK Solve stiff differential equations, Rosenbrock method # with provided coefficients. -function oderosenbrock(F, x0, tspan, gamma, a, b, c; jacobian=nothing, kwargs...) - +function oderosenbrock(F, x0, tspan, gamma, a, b, c; jacobian = nothing, kwargs...) if typeof(jacobian) == Function G = jacobian else - G = (t, x)->fdjacobian(F, x, t) + G = (t, x) -> fdjacobian(F, x, t) end h = diff(tspan) - x = Vector{typeof(x0)}(undef,length(tspan)) + x = Vector{typeof(x0)}(undef, length(tspan)) x[1] = x0 solstep = 1 @@ -382,65 +383,68 @@ function oderosenbrock(F, x0, tspan, gamma, a, b, c; jacobian=nothing, kwargs... xs = x[solstep] dFdx = G(ts, xs) - jac = I/(gamma*hs) - dFdx + jac = I / (gamma * hs) - dFdx - g = Vector{typeof(x0)}(undef,size(a,1)) - g[1] = (jac \ F(ts + b[1]*hs, xs)) - x[solstep+1] = x[solstep] + b[1]*g[1] + g = Vector{typeof(x0)}(undef, size(a, 1)) + g[1] = (jac \ F(ts + b[1] * hs, xs)) + x[solstep + 1] = x[solstep] + b[1] * g[1] - for i = 2:size(a,1) + for i in 2:size(a, 1) dx = zero(x0) - dF = zero(x0/hs) - for j = 1:i-1 - dx += a[i,j]*g[j] - dF += c[i,j]*g[j] + dF = zero(x0 / hs) + for j in 1:(i - 1) + dx += a[i, j] * g[j] + dF += c[i, j] * g[j] end - g[i] = (jac \ (F(ts + b[i]*hs, xs + dx) + dF/hs)) - x[solstep+1] += b[i]*g[i] + g[i] = (jac \ (F(ts + b[i] * hs, xs + dx) + dF / hs)) + x[solstep + 1] += b[i] * g[i] end solstep += 1 end return vcat(tspan), x end - # Kaps-Rentrop coefficients const kr4_coefficients = (0.231, - [0 0 0 0 - 2 0 0 0 + [0 0 0 0 + 2 0 0 0 4.452470820736 4.16352878860 0 0 4.452470820736 4.16352878860 0 0], - [3.95750374663 4.62489238836 0.617477263873 1.28261294568], - [ 0 0 0 0 - -5.07167533877 0 0 0 - 6.02015272865 0.1597500684673 0 0 - -1.856343618677 -8.50538085819 -2.08407513602 0],) - -ode4s_kr(F, x0, tspan; jacobian=nothing, kwargs...) = oderosenbrock(F, x0, tspan, kr4_coefficients...; jacobian=jacobian, kwargs...) + [3.95750374663 4.62489238836 0.617477263873 1.28261294568], + [0 0 0 0 + -5.07167533877 0 0 0 + 6.02015272865 0.1597500684673 0 0 + -1.856343618677 -8.50538085819 -2.08407513602 0]) + +function ode4s_kr(F, x0, tspan; jacobian = nothing, kwargs...) + oderosenbrock(F, x0, tspan, kr4_coefficients...; jacobian = jacobian, kwargs...) +end # Shampine coefficients const s4_coefficients = (0.5, - [ 0 0 0 0 - 2 0 0 0 + [0 0 0 0 + 2 0 0 0 48/25 6/25 0 0 48/25 6/25 0 0], - [19/9 1/2 25/108 125/108], - [ 0 0 0 0 - -8 0 0 0 - 372/25 12/5 0 0 - -112/125 -54/125 -2/5 0],) - -ode4s_s(F, x0, tspan; jacobian=nothing, kwargs...) = - oderosenbrock(F, x0, tspan, s4_coefficients...; jacobian=jacobian, kwargs...) + [19 / 9 1 / 2 25 / 108 125 / 108], + [0 0 0 0 + -8 0 0 0 + 372/25 12/5 0 0 + -112/125 -54/125 -2/5 0]) + +function ode4s_s(F, x0, tspan; jacobian = nothing, kwargs...) + oderosenbrock(F, x0, tspan, s4_coefficients...; jacobian = jacobian, kwargs...) +end # Use Shampine coefficients by default (matching Numerical Recipes) -ode4s(F, x0, tspan; jacobian=nothing, kwargs...) = - ode4s_s(F, x0, tspan; jacobian=nothing, kwargs...) +function ode4s(F, x0, tspan; jacobian = nothing, kwargs...) + ode4s_s(F, x0, tspan; jacobian = nothing, kwargs...) +end -const ms_coefficients4 = [ 1 0 0 0 - -1/2 3/2 0 0 - 5/12 -4/3 23/12 0 - -9/24 37/24 -59/24 55/24] +const ms_coefficients4 = [1 0 0 0 + -1/2 3/2 0 0 + 5/12 -4/3 23/12 0 + -9/24 37/24 -59/24 55/24] ####### Common Interface Bindings diff --git a/src/common.jl b/src/common.jl index 2eb8d2e8b..c945917b7 100644 --- a/src/common.jl +++ b/src/common.jl @@ -1,27 +1,25 @@ -function solve( - prob::DiffEqBase.AbstractODEProblem{uType,tupType,isinplace}, - alg::AlgType, - timeseries=[], ts=[], ks=[]; - verbose=true, - save_timeseries=nothing, - saveat=eltype(tupType)[], reltol = 1e-5, abstol = 1e-8, - save_everystep=isempty(saveat), - dense = save_everystep && isempty(saveat), - save_start = save_everystep || isempty(saveat) || typeof(saveat) <: Number ? - true : prob.tspan[1] in saveat, - callback=nothing, - dtmin = abs(prob.tspan[2]-prob.tspan[1])/1e-9, - dtmax = abs(prob.tspan[2]-prob.tspan[1])/2.5, - timeseries_errors=true, dense_errors=false, - dt = 0.0, norm = LinearAlgebra.norm, - kwargs...) where {uType,tupType,isinplace,AlgType<:ODEjlAlgorithm} - +function solve(prob::DiffEqBase.AbstractODEProblem{uType, tupType, isinplace}, + alg::AlgType, + timeseries = [], ts = [], ks = []; + verbose = true, + save_timeseries = nothing, + saveat = eltype(tupType)[], reltol = 1e-5, abstol = 1e-8, + save_everystep = isempty(saveat), + dense = save_everystep && isempty(saveat), + save_start = save_everystep || isempty(saveat) || typeof(saveat) <: Number ? + true : prob.tspan[1] in saveat, + callback = nothing, + dtmin = abs(prob.tspan[2] - prob.tspan[1]) / 1e-9, + dtmax = abs(prob.tspan[2] - prob.tspan[1]) / 2.5, + timeseries_errors = true, dense_errors = false, + dt = 0.0, norm = LinearAlgebra.norm, + kwargs...) where {uType, tupType, isinplace, AlgType <: ODEjlAlgorithm} tType = eltype(tupType) if verbose warned = !isempty(kwargs) && check_keywords(alg, kwargs, warnlist) if !(typeof(prob.f) <: DiffEqBase.AbstractParameterizedFunction) && - typeof(alg) <: ode23s + typeof(alg) <: ode23s if DiffEqBase.has_tgrad(prob.f) @warn("Explicit t-gradient given to this stiff solver is ignored.") warned = true @@ -50,11 +48,14 @@ function solve( tspan = prob.tspan if typeof(saveat) <: Number - if (tspan[1]:saveat:tspan[end])[end] == tspan[end] - saveat_vec = convert(Vector{tType},collect(tType,tspan[1]+saveat:saveat:tspan[end])) - else - saveat_vec = convert(Vector{tType},collect(tType,tspan[1]+saveat:saveat:(tspan[end]-saveat))) - end + if (tspan[1]:saveat:tspan[end])[end] == tspan[end] + saveat_vec = convert(Vector{tType}, + collect(tType, (tspan[1] + saveat):saveat:tspan[end])) + else + saveat_vec = convert(Vector{tType}, + collect(tType, + (tspan[1] + saveat):saveat:(tspan[end] - saveat))) + end else saveat_vec = convert(Vector{tType}, collect(saveat)) end @@ -64,9 +65,9 @@ function solve( end if !isempty(saveat_vec) && saveat_vec[1] == tspan[1] - Ts = unique([saveat_vec;tspan[2]]) + Ts = unique([saveat_vec; tspan[2]]) else - Ts = unique([tspan[1];saveat_vec;tspan[2]]) + Ts = unique([tspan[1]; saveat_vec; tspan[2]]) end if save_everystep @@ -79,25 +80,25 @@ function solve( p = prob.p if isinplace - f = (t,u) -> (du = zero(u); prob.f(du,u,p,t); vec(du)) + f = (t, u) -> (du = zero(u); prob.f(du, u, p, t); vec(du)) elseif uType <: AbstractArray - f = (t,u) -> vec(prob.f(reshape(u,sizeu),p,t)) + f = (t, u) -> vec(prob.f(reshape(u, sizeu), p, t)) else - f = (t,u) -> prob.f(u,p,t) + f = (t, u) -> prob.f(u, p, t) end u0 = uType <: AbstractArray ? vec(prob.u0) : prob.u0 # Calling the solver, i.e. if the algorithm is ode45, # then AlgType(...) is ode45(...) - ts, timeseries_tmp = AlgType(f,u0,Ts; + ts, timeseries_tmp = AlgType(f, u0, Ts; norm = norm, - abstol=abstol, - reltol=reltol, - maxstep=dtmax, - minstep=dtmin, - initstep=dt, - points=points) + abstol = abstol, + reltol = reltol, + maxstep = dtmax, + minstep = dtmin, + initstep = dt, + points = points) if save_start start_idx = 1 # The index to start making the timeseries from @@ -109,15 +110,15 @@ function solve( # Reshape the result if needed if uType <: AbstractArray timeseries = Vector{uType}() - for i=start_idx:length(timeseries_tmp) - push!(timeseries,reshape(timeseries_tmp[i],sizeu)) + for i in start_idx:length(timeseries_tmp) + push!(timeseries, reshape(timeseries_tmp[i], sizeu)) end else timeseries = timeseries_tmp end - DiffEqBase.build_solution(prob,alg,ts,timeseries, - timeseries_errors = timeseries_errors, - dense_errors = dense_errors, - retcode = :Succss) + DiffEqBase.build_solution(prob, alg, ts, timeseries, + timeseries_errors = timeseries_errors, + dense_errors = dense_errors, + retcode = :Succss) end # function solve diff --git a/src/runge_kutta.jl b/src/runge_kutta.jl index 3e872d208..4b738c5dc 100644 --- a/src/runge_kutta.jl +++ b/src/runge_kutta.jl @@ -13,149 +13,156 @@ struct TableauRKExplicit{Name, S, T} <: Tableau{Name, S, T} # second for error calc. b::Matrix{T} c::Vector{T} - function TableauRKExplicit{Name, S, T}(order,a,b,c) where {Name, S, T} - @assert isa(S,Integer) - @assert isa(Name,Symbol) - @assert c[1]==0 + function TableauRKExplicit{Name, S, T}(order, a, b, c) where {Name, S, T} + @assert isa(S, Integer) + @assert isa(Name, Symbol) + @assert c[1] == 0 @assert istril(a) - @assert S==length(c)==size(a,1)==size(a,2)==size(b,2) - @assert size(b,1)==length(order) + @assert S == length(c) == size(a, 1) == size(a, 2) == size(b, 2) + @assert size(b, 1) == length(order) # consistency: - if T<:AbstractFloat - @assert norm(sum(a,dims=2)-c'',Inf) < 100*eps(T) + if T <: AbstractFloat + @assert norm(sum(a, dims = 2) - c'', Inf) < 100 * eps(T) else - @assert norm(sum(a,dims=2)-c'',Inf) < 100*eps(Float64) + @assert norm(sum(a, dims = 2) - c'', Inf) < 100 * eps(Float64) end - new{Name, S, T}(order,a,b,c) + new{Name, S, T}(order, a, b, c) end end function TableauRKExplicit(name::Symbol, order::(Tuple{Vararg{Int}}), - a::Matrix{T}, b::Matrix{T}, c::Vector{T}) where T - TableauRKExplicit{name,length(c),T}(order, a, b, c) + a::Matrix{T}, b::Matrix{T}, c::Vector{T}) where {T} + TableauRKExplicit{name, length(c), T}(order, a, b, c) end function TableauRKExplicit(name::Symbol, order::(Tuple{Vararg{Int}}), T::Type, - a::Matrix, b::Matrix, c::Vector) - TableauRKExplicit{name,length(c),T}(order, convert(Matrix{T},a), - convert(Matrix{T},b), convert(Vector{T},c) ) + a::Matrix, b::Matrix, c::Vector) + TableauRKExplicit{name, length(c), T}(order, convert(Matrix{T}, a), + convert(Matrix{T}, b), convert(Vector{T}, c)) end -conv_field(D,a::Array{T,N}) where {T,N} = convert(Array{D,N}, a) -function Base.convert(::Type{Tnew}, tab::TableauRKExplicit{Name,S,T}) where {Tnew<:Real,Name,S,T} +conv_field(D, a::Array{T, N}) where {T, N} = convert(Array{D, N}, a) +function Base.convert(::Type{Tnew}, + tab::TableauRKExplicit{Name, S, T}) where {Tnew <: Real, Name, S, T} # Converts the tableau coefficients to the new type Tnew newflds = () for n in fieldnames(typeof(tab)) - fld = getfield(tab,n) - if eltype(fld)==T + fld = getfield(tab, n) + if eltype(fld) == T newflds = tuple(newflds..., conv_field(Tnew, fld)) else newflds = tuple(newflds..., fld) end end - TableauRKExplicit{Name,S,Tnew}(newflds...) # TODO: could this be done more generically in a type-stable way? + TableauRKExplicit{Name, S, Tnew}(newflds...) # TODO: could this be done more generically in a type-stable way? end - isexplicit(b::TableauRKExplicit) = istril(b.a) # Test whether it's an explicit method -isadaptive(b::TableauRKExplicit) = size(b.b, 1)==2 - +isadaptive(b::TableauRKExplicit) = size(b.b, 1) == 2 # First same as last. Means ks[:,end]=ks_nextstep[:,1], c.f. H&W p.167 -isFSAL(btab::TableauRKExplicit) = btab.a[end,:]==btab.b[1,:] && btab.c[end]==1 # the latter is not needed really +isFSAL(btab::TableauRKExplicit) = btab.a[end, :] == btab.b[1, :] && btab.c[end] == 1 # the latter is not needed really ## Tableaus for explicit RK methods # Fixed step: -const bt_feuler = TableauRKExplicit(:feuler,(1,), Rational{Int64}, - zeros(Int,1,1), - [1][:,:], - [0] - ) -const bt_midpoint = TableauRKExplicit(:midpoint,(2,), Rational{Int64}, - [0 0 - 1//2 0], - [0 1], - [0, 1//2] - ) -const bt_heun = TableauRKExplicit(:heun,(2,), Rational{Int64}, - [0 0 - 1 0], - [1//2 1//2], +const bt_feuler = TableauRKExplicit(:feuler, (1,), Rational{Int64}, + zeros(Int, 1, 1), + [1][:, :], + [0]) +const bt_midpoint = TableauRKExplicit(:midpoint, (2,), Rational{Int64}, + [0 0 + 1//2 0], + [0 1], + [0, 1 // 2]) +const bt_heun = TableauRKExplicit(:heun, (2,), Rational{Int64}, + [0 0 + 1 0], + [1 // 2 1 // 2], [0, 1]) -const bt_rk4 = TableauRKExplicit(:rk4,(4,),Rational{Int64}, - [0 0 0 0 - 1//2 0 0 0 - 0 1//2 0 0 - 0 0 1 0], - [1//6 1//3 1//3 1//6], - [0, 1//2, 1//2, 1]) +const bt_rk4 = TableauRKExplicit(:rk4, (4,), Rational{Int64}, + [0 0 0 0 + 1//2 0 0 0 + 0 1//2 0 0 + 0 0 1 0], + [1 // 6 1 // 3 1 // 3 1 // 6], + [0, 1 // 2, 1 // 2, 1]) # Adaptive step: # Heun Euler https://en.wikipedia.org/wiki/Runge–Kutta_methods -const bt_rk21 = TableauRKExplicit(:heun_euler,(2,1), Rational{Int64}, - [0 0 - 1 0], - [1//2 1//2 - 1 0], - [0, 1]) +const bt_rk21 = TableauRKExplicit(:heun_euler, (2, 1), Rational{Int64}, + [0 0 + 1 0], + [1//2 1//2 + 1 0], + [0, 1]) # Bogacki–Shampine coefficients -const bt_rk23 = TableauRKExplicit(:bogacki_shampine,(2,3), Rational{Int64}, - [0 0 0 0 - 1//2 0 0 0 - 0 3//4 0 0 - 2//9 1//3 4//9 0], +const bt_rk23 = TableauRKExplicit(:bogacki_shampine, (2, 3), Rational{Int64}, + [0 0 0 0 + 1//2 0 0 0 + 0 3//4 0 0 + 2//9 1//3 4//9 0], [7//24 1//4 1//3 1//8 2//9 1//3 4//9 0], - [0, 1//2, 3//4, 1] - ) + [0, 1 // 2, 3 // 4, 1]) # Fehlberg https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta%E2%80%93Fehlberg_method -const bt_rk45 = TableauRKExplicit(:fehlberg,(4,5),Rational{Int64}, - [ 0 0 0 0 0 0 - 1//4 0 0 0 0 0 - 3//32 9//32 0 0 0 0 - 1932//2197 -7200//2197 7296//2197 0 0 0 - 439//216 -8 3680//513 -845//4104 0 0 - -8//27 2 -3544//2565 1859//4104 -11//40 0 ], - [25//216 0 1408//2565 2197//4104 -1//5 0 - 16//135 0 6656//12825 28561//56430 -9//50 2//55], - [0, 1//4, 3//8, 12//13, 1, 1//2]) +const bt_rk45 = TableauRKExplicit(:fehlberg, (4, 5), Rational{Int64}, + [0 0 0 0 0 0 + 1//4 0 0 0 0 0 + 3//32 9//32 0 0 0 0 + 1932//2197 -7200//2197 7296//2197 0 0 0 + 439//216 -8 3680//513 -845//4104 0 0 + -8//27 2 -3544//2565 1859//4104 -11//40 0], + [25//216 0 1408//2565 2197//4104 -1//5 0 + 16//135 0 6656//12825 28561//56430 -9//50 2//55], + [0, 1 // 4, 3 // 8, 12 // 13, 1, 1 // 2]) # Dormand-Prince https://en.wikipedia.org/wiki/Dormand%E2%80%93Prince_method -const bt_dopri5 = TableauRKExplicit(:dopri, (5,4), Rational{Int64}, - [0 0 0 0 0 0 0 - 1//5 0 0 0 0 0 0 - 3//40 9//40 0 0 0 0 0 - 44//45 -56//15 32//9 0 0 0 0 - 19372//6561 -25360//2187 64448//6561 -212//729 0 0 0 - 9017//3168 -355//33 46732//5247 49//176 -5103//18656 0 0 - 35//384 0 500//1113 125//192 -2187//6784 11//84 0], - [35//384 0 500//1113 125//192 -2187//6784 11//84 0 - 5179//57600 0 7571//16695 393//640 -92097//339200 187//2100 1//40], - [0, 1//5, 3//10, 4//5, 8//9, 1, 1] - ) +const bt_dopri5 = TableauRKExplicit(:dopri, (5, 4), Rational{Int64}, + [0 0 0 0 0 0 0 + 1//5 0 0 0 0 0 0 + 3//40 9//40 0 0 0 0 0 + 44//45 -56//15 32//9 0 0 0 0 + 19372//6561 -25360//2187 64448//6561 -212//729 0 0 0 + 9017//3168 -355//33 46732//5247 49//176 -5103//18656 0 0 + 35//384 0 500//1113 125//192 -2187//6784 11//84 0], + [35//384 0 500//1113 125//192 -2187//6784 11//84 0 + 5179//57600 0 7571//16695 393//640 -92097//339200 187//2100 1//40], + [0, 1 // 5, 3 // 10, 4 // 5, 8 // 9, 1, 1]) # Fehlberg 7(8) coefficients # Values from pag. 65, Fehlberg, Erwin. "Classical fifth-, sixth-, seventh-, and eighth-order Runge-Kutta formulas with stepsize control". # National Aeronautics and Space Administration. -const bt_feh78 = TableauRKExplicit(:feh78, (7,8), Rational{Int64}, - [ 0 0 0 0 0 0 0 0 0 0 0 0 0 - 2//27 0 0 0 0 0 0 0 0 0 0 0 0 - 1//36 1//12 0 0 0 0 0 0 0 0 0 0 0 - 1//24 0 1//8 0 0 0 0 0 0 0 0 0 0 - 5//12 0 -25//16 25//16 0 0 0 0 0 0 0 0 0 - 1//20 0 0 1//4 1//5 0 0 0 0 0 0 0 0 - -25//108 0 0 125//108 -65//27 125//54 0 0 0 0 0 0 0 - 31//300 0 0 0 61//225 -2//9 13//900 0 0 0 0 0 0 - 2 0 0 -53//6 704//45 -107//9 67//90 3 0 0 0 0 0 - -91//108 0 0 23//108 -976//135 311//54 -19//60 17//6 -1//12 0 0 0 0 - 2383//4100 0 0 -341//164 4496//1025 -301//82 2133//4100 45//82 45//164 18//41 0 0 0 - 3//205 0 0 0 0 -6//41 -3//205 -3//41 3//41 6//41 0 0 0 - -1777//4100 0 0 -341//164 4496//1025 -289//82 2193//4100 51//82 33//164 12//41 0 1 0], - [41//840 0 0 0 0 34//105 9//35 9//35 9//280 9//280 41//840 0 0 - 0 0 0 0 0 34//105 9//35 9//35 9//280 9//280 0 41//840 41//840], - [0, 2//27, 1//9, 1//6 , 5//12, 1//2 , 5//6 , 1//6 , 2//3 , 1//3 , 1 , 0, 1] - ) - +const bt_feh78 = TableauRKExplicit(:feh78, (7, 8), Rational{Int64}, + [0 0 0 0 0 0 0 0 0 0 0 0 0 + 2//27 0 0 0 0 0 0 0 0 0 0 0 0 + 1//36 1//12 0 0 0 0 0 0 0 0 0 0 0 + 1//24 0 1//8 0 0 0 0 0 0 0 0 0 0 + 5//12 0 -25//16 25//16 0 0 0 0 0 0 0 0 0 + 1//20 0 0 1//4 1//5 0 0 0 0 0 0 0 0 + -25//108 0 0 125//108 -65//27 125//54 0 0 0 0 0 0 0 + 31//300 0 0 0 61//225 -2//9 13//900 0 0 0 0 0 0 + 2 0 0 -53//6 704//45 -107//9 67//90 3 0 0 0 0 0 + -91//108 0 0 23//108 -976//135 311//54 -19//60 17//6 -1//12 0 0 0 0 + 2383//4100 0 0 -341//164 4496//1025 -301//82 2133//4100 45//82 45//164 18//41 0 0 0 + 3//205 0 0 0 0 -6//41 -3//205 -3//41 3//41 6//41 0 0 0 + -1777//4100 0 0 -341//164 4496//1025 -289//82 2193//4100 51//82 33//164 12//41 0 1 0], + [41//840 0 0 0 0 34//105 9//35 9//35 9//280 9//280 41//840 0 0 + 0 0 0 0 0 34//105 9//35 9//35 9//280 9//280 0 41//840 41//840], + [ + 0, + 2 // 27, + 1 // 9, + 1 // 6, + 5 // 12, + 1 // 2, + 5 // 6, + 1 // 6, + 2 // 3, + 1 // 3, + 1, + 0, + 1, + ]) ################################ # Fixed step Runge-Kutta methods @@ -170,11 +177,11 @@ ode4(fn, y0, tspan; kwargs...) = oderk_fixed(fn, y0, tspan, bt_rk4; kwargs...) function oderk_fixed(fn, y0, tspan, btab::TableauRKExplicit; kwargs...) # Non-arrays y0 treat as scalar fn_(t, y) = [fn(t, y[1])] - t,y = oderk_fixed(fn_, [y0], tspan, btab) + t, y = oderk_fixed(fn_, [y0], tspan, btab) return t, vcat_nosplat(y) end function oderk_fixed(fn, y0::AbstractVector, tspan, - btab_::TableauRKExplicit{N,S}, alias_u0=false) where {N,S} + btab_::TableauRKExplicit{N, S}, alias_u0 = false) where {N, S} # TODO: instead of AbstractVector use a Holy-trait # Needed interface: @@ -185,7 +192,7 @@ function oderk_fixed(fn, y0::AbstractVector, tspan, Et, Eyf, Ty, btab = make_consistent_types(fn, y0, tspan, btab_) dof = length(y0) - ys = Vector{Ty}(undef,length(tspan)) + ys = Vector{Ty}(undef, length(tspan)) allocate!(ys, y0, dof) if alias_u0 ys[1] = y0 @@ -195,16 +202,16 @@ function oderk_fixed(fn, y0::AbstractVector, tspan, tspan = convert(Vector{Et}, tspan) # work arrays: - ks = Vector{Ty}(undef,S) + ks = Vector{Ty}(undef, S) # allocate!(ks, y0, dof) # no need to allocate as fn is not in-place ytmp = similar(y0, Eyf, dof) - for i=1:length(tspan)-1 - dt = tspan[i+1]-tspan[i] - ys[i+1][:] = ys[i] - for s=1:S + for i in 1:(length(tspan) - 1) + dt = tspan[i + 1] - tspan[i] + ys[i + 1][:] = ys[i] + for s in 1:S calc_next_k!(ks, ytmp, ys[i], s, fn, tspan[i], dt, dof, btab) - for d=1:dof - ys[i+1][d] += dt * btab.b[s]*ks[s][d] + for d in 1:dof + ys[i + 1][d] += dt * btab.b[s] * ks[s][d] end end end @@ -226,18 +233,17 @@ ode78(fn, y0, tspan; kwargs...) = oderk_adapt(fn, y0, tspan, bt_feh78; kwargs... function oderk_adapt(fn, y0, tspan, btab::TableauRKExplicit; kwords...) # For y0 which don't support indexing. fn_ = (t, y) -> [fn(t, y[1])] - t,y = oderk_adapt(fn_, [y0], tspan, btab; kwords...) + t, y = oderk_adapt(fn_, [y0], tspan, btab; kwords...) return t, vcat_nosplat(y) end -function oderk_adapt(fn, y0::AbstractVector, tspan, btab_::TableauRKExplicit{N,S}; +function oderk_adapt(fn, y0::AbstractVector, tspan, btab_::TableauRKExplicit{N, S}; reltol = 1.0e-5, abstol = 1.0e-8, - norm=LinearAlgebra.norm, - minstep=abs(tspan[end] - tspan[1])/1e18, - maxstep=abs(tspan[end] - tspan[1])/2.5, - initstep=0, - points=:all, - alias_u0=false - ) where {N,S} + norm = LinearAlgebra.norm, + minstep = abs(tspan[end] - tspan[1]) / 1e18, + maxstep = abs(tspan[end] - tspan[1]) / 2.5, + initstep = 0, + points = :all, + alias_u0 = false) where {N, S} # Needed interface: # On components: # - note that the type of the components might change! @@ -246,12 +252,13 @@ function oderk_adapt(fn, y0::AbstractVector, tspan, btab_::TableauRKExplicit{N,S # For y0 which support indexing. Currently y0<:AbstractVector but # that could be relaxed with a Holy-trait. - !isadaptive(btab_) && error("Can only use this solver with an adaptive RK Butcher table") + !isadaptive(btab_) && + error("Can only use this solver with an adaptive RK Butcher table") Et, Eyf, Ty, btab = make_consistent_types(fn, y0, tspan, btab_) # parameters order = minimum(btab.order) timeout_const = 5 # after step reduction do not increase step for - # timeout_const steps + # timeout_const steps ## Initialization dof = length(y0) @@ -264,18 +271,18 @@ function oderk_adapt(fn, y0::AbstractVector, tspan, btab_::TableauRKExplicit{N,S if alias_u0 y = y0 # y at time t else - y = similar(y0, Eyf, dof) # y at time t + y = similar(y0, Eyf, dof) # y at time t end - y[:] = y0 + y[:] = y0 ytrial = similar(y0, Eyf, dof) # trial solution at time t+dt - yerr = similar(y0, Eyf, dof) # error of trial solution - ks = Vector{Ty}(undef,S) + yerr = similar(y0, Eyf, dof) # error of trial solution + ks = Vector{Ty}(undef, S) # allocate!(ks, y0, dof) # no need to allocate as fn is not in-place - ytmp = similar(y0, Eyf, dof) + ytmp = similar(y0, Eyf, dof) # output ys nsteps_fixed = length(tspan) # these are always output - ys = Vector{Ty}(undef,nsteps_fixed) + ys = Vector{Ty}(undef, nsteps_fixed) allocate!(ys, y0, dof) ys[1] = y0 @@ -290,16 +297,16 @@ function oderk_adapt(fn, y0::AbstractVector, tspan, btab_::TableauRKExplicit{N,S end # Time dt, tdir, ks[1] = hinit(fn, y, tstart, tend, order, reltol, abstol) # sets ks[1]=f0 - if initstep!=0 - dt = sign(initstep)==tdir ? initstep : error("initstep has wrong sign.") + if initstep != 0 + dt = sign(initstep) == tdir ? initstep : error("initstep has wrong sign.") end # Diagnostics dts = Et[] errs = Float64[] - steps = [0,0] # [accepted, rejected] + steps = [0, 0] # [accepted, rejected] ## Integration loop - islaststep = abs(t+dt-tend)<=eps(tend) ? true : false + islaststep = abs(t + dt - tend) <= eps(tend) ? true : false timeout = 0 # for step-control iter = 2 # the index into tspan and ys while true @@ -307,26 +314,28 @@ function oderk_adapt(fn, y0::AbstractVector, tspan, btab_::TableauRKExplicit{N,S rk_embedded_step!(ytrial, yerr, ks, ytmp, y, fn, t, dt, dof, btab) # Check error and find a new step size: err, newdt, timeout = stepsize_hw92!(dt, tdir, y, ytrial, yerr, order, timeout, - dof, abstol, reltol, maxstep, norm) + dof, abstol, reltol, maxstep, norm) - if err<=1 # accept step + if err <= 1 # accept step # diagnostics - steps[1] +=1 + steps[1] += 1 push!(dts, dt) push!(errs, err) # Output: f0 = ks[1] - f1 = isFSAL(btab) ? ks[S] : fn(t+dt, ytrial) - if points==:specified + f1 = isFSAL(btab) ? ks[S] : fn(t + dt, ytrial) + if points == :specified # interpolate onto given output points - while iter-1= tdir*tend - dt = tend-t + if tdir * (t + dt + dt / 100) >= tdir * tend + dt = tend - t islaststep = true # next step is the last, if it succeeds end - elseif abs(newdt)0 + newdt = min(maxstep, tdir * dt * max(facmin, fac * (1 / err)^(T(1 / (order + 1))))) # Eq 4.13 modified + if timeout > 0 newdt = min(newdt, dt) timeout -= 1 end - return err, tdir*newdt, timeout + return err, tdir * newdt, timeout end -function calc_next_k!(ks::Vector, ytmp::Ty, y, s, fn, t, dt, dof, btab) where Ty +function calc_next_k!(ks::Vector, ytmp::Ty, y, s, fn, t, dt, dof, btab) where {Ty} # Calculates the next ks and puts it into ks[s] # - ks and ytmp are modified inside this function. @@ -450,49 +460,51 @@ function calc_next_k!(ks::Vector, ytmp::Ty, y, s, fn, t, dt, dof, btab) where Ty # On y0 container: setindex!, getindex, fn ytmp[:] = y - for ss=1:s-1, d=1:dof - ytmp[d] += dt * ks[ss][d] * btab.a[s,ss] + for ss in 1:(s - 1), d in 1:dof + ytmp[d] += dt * ks[ss][d] * btab.a[s, ss] end - ks[s] = fn(t + btab.c[s]*dt, ytmp)::Ty + ks[s] = fn(t + btab.c[s] * dt, ytmp)::Ty nothing end # Helper functions: -function allocate!(vec::Vector{T}, y0, dof) where T +function allocate!(vec::Vector{T}, y0, dof) where {T} # Allocates all vectors inside a Vector{Vector} using the same # kind of container as y0 has and element type eltype(eltype(vec)). - for s=1:length(vec) + for s in 1:length(vec) vec[s] = similar(y0, eltype(T), dof) end end function index_or_push!(vec, i, val) # Fills in the vector until there is no space, then uses push! # instead. - if length(vec)>=i + if length(vec) >= i vec[i] = val else push!(vec, val) end end -vcat_nosplat(y) = eltype(y[1])[el[1] for el in y] # Does vcat(y...) without the splatting +vcat_nosplat(y) = eltype(y[1])[el[1] for el in y] # Does vcat(y...) without the splatting -function hermite_interp!(y, tquery,t,dt,y0,y1,f0,f1) +function hermite_interp!(y, tquery, t, dt, y0, y1, f0, f1) # For dense output see Hairer & Wanner p.190 using Hermite # interpolation. Updates y in-place. # # f_0 = f(x_0 , y_0) , f_1 = f(x_0 + h, y_1 ) # this is O(3). TODO for higher order. - theta = (tquery-t)/dt - for i=1:length(y0) - y[i] = ((1-theta)*y0[i] + theta*y1[i] + theta*(theta-1) * - ((1-2*theta)*(y1[i]-y0[i]) + (theta-1)*dt*f0[i] + theta*dt*f1[i]) ) + theta = (tquery - t) / dt + for i in 1:length(y0) + y[i] = ((1 - theta) * y0[i] + theta * y1[i] + + theta * (theta - 1) * + ((1 - 2 * theta) * (y1[i] - y0[i]) + (theta - 1) * dt * f0[i] + + theta * dt * f1[i])) end nothing end -function hermite_interp(tquery,t,dt,y0,y1,f0,f1) +function hermite_interp(tquery, t, dt, y0, y1, f0, f1) # Returns the y instead of in-place y = similar(y0) - hermite_interp!(y,tquery,t,dt,y0,y1,f0,f1) + hermite_interp!(y, tquery, t, dt, y0, y1, f0, f1) return y end diff --git a/test/common.jl b/test/common.jl index d58b10209..03cc91105 100644 --- a/test/common.jl +++ b/test/common.jl @@ -1,30 +1,32 @@ using ODE, DiffEqBase, DiffEqProblemLibrary, Test -using DiffEqProblemLibrary.ODEProblemLibrary: importodeproblems; importodeproblems() +using DiffEqProblemLibrary.ODEProblemLibrary: importodeproblems; +importodeproblems(); import DiffEqProblemLibrary.ODEProblemLibrary: prob_ode_linear, prob_ode_2Dlinear using LinearAlgebra -dt=1/2^(4) +dt = 1 / 2^(4) -algs = [ode23(),ode45(),ode78(),ode4(),ode4ms(),ode4s()] # no ode23s +algs = [ode23(), ode45(), ode78(), ode4(), ode4ms(), ode4s()] # no ode23s # Check for errors prob = prob_ode_linear for alg in algs - sol =solve(prob,alg;dt=dt,abstol=1e-6,reltol=1e-3) + sol = solve(prob, alg; dt = dt, abstol = 1e-6, reltol = 1e-3) @test typeof(sol[2]) <: Number end prob = prob_ode_2Dlinear for alg in algs - sol =solve(prob,alg;dt=dt,dtmin=eps(),abstol=1e-6,reltol=1e-3) - @test size(sol[2]) == (4,2) + sol = solve(prob, alg; dt = dt, dtmin = eps(), abstol = 1e-6, reltol = 1e-3) + @test size(sol[2]) == (4, 2) end # Check that warnings are displayed for alg in algs - sol = solve(prob,alg;dt=dt,abstol=1e-6,reltol=1e-3, dense=true, verbose=true, + sol = solve(prob, alg; dt = dt, abstol = 1e-6, reltol = 1e-3, dense = true, + verbose = true, save_idxs = true, progress = true, beta1 = 1.23, gamma = 0.5) end diff --git a/test/interface-tests.jl b/test/interface-tests.jl index a6ef8d510..e45a27152 100644 --- a/test/interface-tests.jl +++ b/test/interface-tests.jl @@ -7,36 +7,38 @@ import Base: +, -, *, / using LinearAlgebra -const delta0 = 0. -const V0 = 1. -const g0 = 0. +const delta0 = 0.0 +const V0 = 1.0 +const g0 = 0.0 # define custom type ... struct CompSol <: AbstractMatrix{ComplexF64} - rho::Matrix{ComplexF64} - x::Float64 - p::Float64 + rho::Matrix{ComplexF64} + x::Float64 + p::Float64 - CompSol(r,x,p) = new(copy(r),x,p) + CompSol(r, x, p) = new(copy(r), x, p) end -Base.getindex(A::CompSol,i...) = A.rho[i...] -Base.setindex!(A,x,i...) = setindex!(A.rho,x,i...) +Base.getindex(A::CompSol, i...) = A.rho[i...] +Base.setindex!(A, x, i...) = setindex!(A.rho, x, i...) # ... which has to support the following operations # to work with odeX -LinearAlgebra.norm(y::CompSol, p::Float64) = maximum([LinearAlgebra.norm(y.rho, p) abs(y.x) abs(y.p)]) +function LinearAlgebra.norm(y::CompSol, p::Float64) + maximum([LinearAlgebra.norm(y.rho, p) abs(y.x) abs(y.p)]) +end LinearAlgebra.norm(y::CompSol) = LinearAlgebra.norm(y::CompSol, 2.0) -+(y1::CompSol, y2::CompSol) = CompSol(y1.rho+y2.rho, y1.x+y2.x, y1.p+y2.p) --(y1::CompSol, y2::CompSol) = CompSol(y1.rho-y2.rho, y1.x-y2.x, y1.p-y2.p) -*(y1::CompSol, s::Real) = CompSol(y1.rho*s, y1.x*s, y1.p*s) -*(s::Real, y1::CompSol) = y1*s -/(y1::CompSol, s::Real) = CompSol(y1.rho/s, y1.x/s, y1.p/s) ++(y1::CompSol, y2::CompSol) = CompSol(y1.rho + y2.rho, y1.x + y2.x, y1.p + y2.p) +-(y1::CompSol, y2::CompSol) = CompSol(y1.rho - y2.rho, y1.x - y2.x, y1.p - y2.p) +*(y1::CompSol, s::Real) = CompSol(y1.rho * s, y1.x * s, y1.p * s) +*(s::Real, y1::CompSol) = y1 * s +/(y1::CompSol, s::Real) = CompSol(y1.rho / s, y1.x / s, y1.p / s) ### new for PR #68 -Base.abs(y::CompSol) = LinearAlgebra.norm(y, 2.) # TODO not needed anymore once https://github.com/JuliaLang/julia/pull/11043 is in current stable julia -Base.zero(::Type{CompSol}) = CompSol(complex(zeros(2,2)), 0., 0.) +Base.abs(y::CompSol) = LinearAlgebra.norm(y, 2.0) # TODO not needed anymore once https://github.com/JuliaLang/julia/pull/11043 is in current stable julia +Base.zero(::Type{CompSol}) = CompSol(complex(zeros(2, 2)), 0.0, 0.0) ODE.isoutofdomain(y::CompSol) = any(isnan, vcat(y.rho[:], y.x, y.p)) ################################################################################ @@ -44,32 +46,33 @@ ODE.isoutofdomain(y::CompSol) = any(isnan, vcat(y.rho[:], y.x, y.p)) # define RHSs of differential equations # delta, V and g are parameters function rhs(t, y, delta, V, g) - H = [[-delta/2 V]; [V delta/2]] + H = [[-delta / 2 V]; [V delta / 2]] - rho_dot = -im*H*y.rho + im*y.rho*H - x_dot = y.p - p_dot = -y.x + rho_dot = -im * H * y.rho + im * y.rho * H + x_dot = y.p + p_dot = -y.x - return CompSol( rho_dot, x_dot, p_dot) + return CompSol(rho_dot, x_dot, p_dot) end # inital conditons -rho0 = zeros(2,2); -rho0[1,1]=1.; -y0 = CompSol(complex(rho0), 2., 1.) +rho0 = zeros(2, 2); +rho0[1, 1] = 1.0; +y0 = CompSol(complex(rho0), 2.0, 1.0) # solve ODEs -endt = 2*pi; +endt = 2 * pi; -t,y1 = ODE.ode45((t,y)->rhs(t, y, delta0, V0, g0), y0, [0., endt]) # used as reference +t, y1 = ODE.ode45((t, y) -> rhs(t, y, delta0, V0, g0), y0, [0.0, endt]) # used as reference print("Testing interface for scalar-like state... ") for solver in solvers # these only work with some Array-like interface defined: if solver in [ODE.ode23s, ODE.ode4s_s, ODE.ode4s_kr] continue end - t,y2 = solver((t,y)->rhs(t, y, delta0, V0, g0), y0, range(0., stop=endt, length=500)) - @test LinearAlgebra.norm(y1[end]-y2[end])<0.15 + t, y2 = solver((t, y) -> rhs(t, y, delta0, V0, g0), y0, + range(0.0, stop = endt, length = 500)) + @test LinearAlgebra.norm(y1[end] - y2[end]) < 0.15 end println("ok.") diff --git a/test/runtests.jl b/test/runtests.jl index 3e5837f0d..d8a2139fc 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -3,33 +3,33 @@ using ODE, Test, LinearAlgebra tols = [5e-2, 1e-2, 1e-2] solvers = [ - ## Non-stiff - # fixed step - ODE.ode1, - ODE.ode2_midpoint, - ODE.ode2_heun, - ODE.ode4, - ODE.ode4ms, - ODE.ode5ms, - # adaptive -# ODE.ode21, # this fails on Travis with 0.4?! TODO revert once fixed. - ODE.ode23, - ODE.ode45_dp, - ODE.ode45_fe, - ODE.ode78, + ## Non-stiff + # fixed step + ODE.ode1, + ODE.ode2_midpoint, + ODE.ode2_heun, + ODE.ode4, + ODE.ode4ms, + ODE.ode5ms, + # adaptive + # ODE.ode21, # this fails on Travis with 0.4?! TODO revert once fixed. + ODE.ode23, + ODE.ode45_dp, + ODE.ode45_fe, + ODE.ode78, - ## Stiff - # fixed-step - ODE.ode4s_s, - ODE.ode4s_kr, - # adaptive - ODE.ode23s] + ## Stiff + # fixed-step + ODE.ode4s_s, + ODE.ode4s_kr, + # adaptive + ODE.ode23s] dtypes = [Float32, Float64, BigFloat] for solver in solvers println("using $solver") - for (tol,T) in zip(tols, dtypes) - if solver==ODE.ode45_fe && T==Float32 + for (tol, T) in zip(tols, dtypes) + if solver == ODE.ode45_fe && T == Float32 # For some reason ode45_fe hangs with Float32! continue end @@ -37,62 +37,61 @@ for solver in solvers # dy # -- = 6 ==> y = 6t # dt - t,y=solver((t,y)->T(6), T(0), T[0:.1:1;]) - @test maximum(abs.(y-6t)) < tol - @test eltype(t)==T + t, y = solver((t, y) -> T(6), T(0), T[0:0.1:1;]) + @test maximum(abs.(y - 6t)) < tol + @test eltype(t) == T # dy # -- = 2t ==> y = t.^2 # dt - t,y=solver((t,y)->T(2)*t, T(0), T[0:.001:1;]) - @test maximum(abs.(y-t.^2)) < tol - @test eltype(t)==T + t, y = solver((t, y) -> T(2) * t, T(0), T[0:0.001:1;]) + @test maximum(abs.(y - t .^ 2)) < tol + @test eltype(t) == T # dy # -- = y ==> y = y0*e.^t # dt - t,y=solver((t,y)->y, T(1), T[0:.001:1;]) - @test maximum(abs.(y-ℯ.^t)) < tol - @test eltype(t)==T + t, y = solver((t, y) -> y, T(1), T[0:0.001:1;]) + @test maximum(abs.(y - ℯ .^ t)) < tol + @test eltype(t) == T - t,y=solver((t,y)->y, T(1), T[1:-.001:0;]) - @test maximum(abs.(y .- ℯ .^ (t.-1))) < tol - @test eltype(t)==T + t, y = solver((t, y) -> y, T(1), T[1:-0.001:0;]) + @test maximum(abs.(y .- ℯ .^ (t .- 1))) < tol + @test eltype(t) == T # dv dw # -- = -w, -- = v ==> v = v0*cos(t) - w0*sin(t), w = w0*cos(t) + v0*sin(t) # dt dt # # y = [v, w] - t,y=solver((t,y)->[-y[2]; y[1]], T[1, 2], T[0:.001:2*pi;]) + t, y = solver((t, y) -> [-y[2]; y[1]], T[1, 2], T[0:0.001:(2 * pi)]) ys = hcat(y...)' # convert Vector{Vector{T}} to Matrix{T} - @test maximum(abs.(ys-[cos.(t)-2*sin.(t) 2*cos.(t)+sin.(t)])) < tol - @test eltype(t)==T + @test maximum(abs.(ys - [cos.(t) - 2 * sin.(t) 2 * cos.(t) + sin.(t)])) < tol + @test eltype(t) == T end end # Test negative starting times ODE.ode23s -@test length(ODE.ode23s((t,y)->[-y[2]; y[1]], [1., 2.], [-5., 0])[1]) > 1 - +@test length(ODE.ode23s((t, y) -> [-y[2]; y[1]], [1.0, 2.0], [-5.0, 0])[1]) > 1 # rober testcase from http://www.unige.ch/~hairer/testset/testset.html let println("ROBER test case") function f(t, y) ydot = similar(y) - ydot[1] = -0.04*y[1] + 1.0e4*y[2]*y[3] - ydot[3] = 3.0e7*y[2]*y[2] + ydot[1] = -0.04 * y[1] + 1.0e4 * y[2] * y[3] + ydot[3] = 3.0e7 * y[2] * y[2] ydot[2] = -ydot[1] - ydot[3] ydot end - t = [0., 1e11] - t,y = ode23s(f, [1.0, 0.0, 0.0], t; abstol=1e-8, reltol=1e-8, - maxstep=1e11/10, minstep=1e11/1e18) + t = [0.0, 1e11] + t, y = ode23s(f, [1.0, 0.0, 0.0], t; abstol = 1e-8, reltol = 1e-8, + maxstep = 1e11 / 10, minstep = 1e11 / 1e18) refsol = [0.2083340149701255e-07, - 0.8333360770334713e-13, - 0.9999999791665050] # reference solution at tspan[2] - @test LinearAlgebra.norm(refsol-y[end], Inf) < 2e-10 + 0.8333360770334713e-13, + 0.9999999791665050] # reference solution at tspan[2] + @test LinearAlgebra.norm(refsol - y[end], Inf) < 2e-10 end include("interface-tests.jl") include("common.jl")