From 0f1aaf6dda9cf58c119b59b66841f9c35f4087ee Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Sun, 29 Jul 2018 23:07:17 -0700 Subject: [PATCH] 07: --- src/ODE.jl | 13 +++++-------- src/common.jl | 2 +- src/runge_kutta.jl | 13 ++++++------- 3 files changed, 12 insertions(+), 16 deletions(-) diff --git a/src/ODE.jl b/src/ODE.jl index 5f59c7d90..570ae88de 100644 --- a/src/ODE.jl +++ b/src/ODE.jl @@ -147,7 +147,6 @@ function make_consistent_types(fn, y0, tspan, btab::Tableau) Eyf = typeof(y0[1]/(tspan[end]-tspan[1])) Et = eltype(tspan) - @show Et @assert Et<:Real if !(Et<:AbstractFloat) Et = promote_type(Et, Float64) @@ -176,7 +175,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)}(length(tspan)) + x = Vector{typeof(x0)}(undef,length(tspan)) x[1] = x0 if 1 <= order <= 4 @@ -289,10 +288,8 @@ function ode23s(F, y0, tspan; h = tdir * min(abs(h), maxstep) y = y0 - tout = Vector{typeof(t)}(undef, 1) - tout[1] = t # first output time - yout = Vector{typeof(y0)}(undef, 1) - yout[1] = deepcopy(y) # first output solution + tout = [t] # first output time + yout = [deepcopy(y)] # first output solution J = jac(t,y) # get Jacobian of F wrt y @@ -376,7 +373,7 @@ function oderosenbrock(F, x0, tspan, gamma, a, b, c; jacobian=nothing, kwargs... end h = diff(tspan) - x = Vector{typeof(x0)}(length(tspan)) + x = Vector{typeof(x0)}(undef,length(tspan)) x[1] = x0 solstep = 1 @@ -388,7 +385,7 @@ function oderosenbrock(F, x0, tspan, gamma, a, b, c; jacobian=nothing, kwargs... jac = I/(gamma*hs) - dFdx - g = Vector{typeof(x0)}(size(a,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] diff --git a/src/common.jl b/src/common.jl index f000d431a..08cdd070b 100644 --- a/src/common.jl +++ b/src/common.jl @@ -13,7 +13,7 @@ function solve( 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 = Base.vecnorm, + dt = 0.0, norm = LinearAlgebra.norm, kwargs...) where {uType,tupType,isinplace,AlgType<:ODEjlAlgorithm} tType = eltype(tupType) diff --git a/src/runge_kutta.jl b/src/runge_kutta.jl index 41032a75e..d8e49b3c0 100644 --- a/src/runge_kutta.jl +++ b/src/runge_kutta.jl @@ -22,9 +22,9 @@ struct TableauRKExplicit{Name, S, T} <: Tableau{Name, S, T} @assert size(b,1)==length(order) # consistency: if T<:AbstractFloat - @assert norm(sum(a,2)-c'',Inf) < 100*eps(T) + @assert norm(sum(a,dims=2)-c'',Inf) < 100*eps(T) else - @assert norm(sum(a,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) end @@ -185,13 +185,13 @@ function oderk_fixed(fn, y0::AbstractVector, tspan, Et, Eyf, Ty, btab = make_consistent_types(fn, y0, tspan, btab_) dof = length(y0) - ys = Vector{Ty}(length(tspan)) + ys = Vector{Ty}(undef,length(tspan)) allocate!(ys, y0, dof) ys[1] = deepcopy(y0) tspan = convert(Vector{Et}, tspan) # work arrays: - ks = Vector{Ty}(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 @@ -242,7 +242,6 @@ 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") - @show tspan Et, Eyf, Ty, btab = make_consistent_types(fn, y0, tspan, btab_) # parameters order = minimum(btab.order) @@ -261,13 +260,13 @@ function oderk_adapt(fn, y0::AbstractVector, tspan, btab_::TableauRKExplicit{N,S 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}(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) # output ys nsteps_fixed = length(tspan) # these are always output - ys = Vector{Ty}(nsteps_fixed) + ys = Vector{Ty}(undef,nsteps_fixed) allocate!(ys, y0, dof) ys[1] = y0