Skip to content
This repository has been archived by the owner on Sep 9, 2024. It is now read-only.

Commit

Permalink
07:
Browse files Browse the repository at this point in the history
  • Loading branch information
ChrisRackauckas committed Jul 30, 2018
1 parent 910498c commit 0f1aaf6
Show file tree
Hide file tree
Showing 3 changed files with 12 additions and 16 deletions.
13 changes: 5 additions & 8 deletions src/ODE.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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]

Expand Down
2 changes: 1 addition & 1 deletion src/common.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
13 changes: 6 additions & 7 deletions src/runge_kutta.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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)
Expand All @@ -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

Expand Down

0 comments on commit 0f1aaf6

Please sign in to comment.