Skip to content

Commit

Permalink
another batch of changes to ode113
Browse files Browse the repository at this point in the history
  • Loading branch information
@obiajulu
obiajulu committed Aug 6, 2016
1 parent 12c592b commit 73b1ea9
Showing 1 changed file with 95 additions and 92 deletions.
187 changes: 95 additions & 92 deletions src/adams_methods.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -211,21 +211,19 @@ function ode113(F, y0::AbstractVector, tspan, order::Integer=4;
t = [tstart]

# initialzation of y[n], which will have variable length as well
y = Array(Ty, 1)
# TODO: fix after https://github.com/JuliaLang/julia/issues/16667
# y[1] = deepcopy(y0)
y[1] = copy(y0)
#y_n = deepcopy(y0)
y_n = copy(y0)

# initialzation of dy[n], which will have variable length as well
_dy = F(t[1],y[1])
Tdy = typeof(_dy)
dy = Tdy[_dy]
dy_n = F(t[1],y0)
Tdy = typeof(dy_n)

# Declare work variables for t_np1, y_np1 and dy_np1 which are for
# for t_{n+1}, y(t_{n+1}), and dy(t_{n+1})
t_np1 = tspan[1]
y_np1 = similar(y0, Eyf, dof)
dy_np1 = similar(dy[1])
dy_np1 = similar(dy_n)

# variables for dense output
tout = [tstart]
Expand All @@ -248,14 +246,14 @@ function ode113(F, y0::AbstractVector, tspan, order::Integer=4;
end

n = 1 #`n` counts the time steps. We start at n=1, the first step
steporder=1 # `steporder` is the order of the AB method at the current step
stageorder=1 # `stageorder`+1 is the order of the AB method at the current step
# Note: that due to Julia starting indexing at 1, whereas Hairer
# indices his variable for order, k, from 0, we have
# k+1 = steporder in what follows
# k+1 = stageorder in what follows

# initialization of arrays used in estimation y_np1 given dy[n],y[n], dt and
# list order timestep are initialized to a maxorder x maxorder dimension
c,g,b,ϕstar, ϕ = intialize_coef_arrays(max_order+3,Tdy,Et,y0)
c,g,b, ϕstar_nm1, ϕstar_n, ϕ_n, ϕ_np1 = intialize_coef_arrays(max_order+3,Tdy,Et)

if max_order <= 13
γ_star = γ_star_array
Expand Down Expand Up @@ -298,20 +296,19 @@ function ode113(F, y0::AbstractVector, tspan, order::Integer=4;
#---################################################################

# Calculation of g_{j}(n)coefficients for j=0,...,k+1;
g_coefs!(g,c,dt,t,t_np1,n,min(n,steporder+2))
g_coefs!(g,c,dt,t,t_np1,n,min(n,stageorder+2))

# Calculation of ϕstar_{j}(n), and ϕ_{j}(n) coef. for j=0,...,k;
ϕ_and_ϕstar_coefs!(b, ϕ, ϕstar,dy[n],t,t_np1,n,1:min(max(n-1,1),steporder+1))
#Note that 1<= min(max(n-1,1),steporder+1) <= n-1 and steporder+1
ϕ_and_ϕstar_coefs!(b, ϕstar_nm1, ϕstar_n, ϕ_n,dy_n,t,t_np1,n,1:min(max(n-1,1),stageorder+1))
#Note that 1<= min(max(n-1,1),stageorder+1) <= n-1 and stageorder+1

# Calculate next value, y_np1, with order = k = steporder + 1. See
# Calculate predictor with explicit method of order `stageorder`. See
# equation (5.5) in Hairer et al, Vol 1, where we note that here
# y_np1 = y_{k}(t_{n+1})

#y_np1 = copy(y[n])
copy!(y_np1,y[n])
for j= 1:steporder
y_np1 += dt*g[j]ϕstar[2,j]
copy!(y_np1,y_n)
for j= 1:stageorder
y_np1 += dt*g[j]*ϕstar_n[j]
end

# Compute current estimate for derivative
Expand All @@ -326,20 +323,20 @@ function ode113(F, y0::AbstractVector, tspan, order::Integer=4;
# desired. As least one loop is necessary, for the sake of the
# adaptive step selection
CE_loops = 1
if n>=steporder+1
if n>=stageorder # try with n >= stageorder
for loop_counter=1:CE_loops
# (C)orrection of y_np1
# Add new coefficients based on this prediction
# computation of ϕ, ϕstar using recurrence relations
ϕ_coefs!(ϕ, ϕstar,dy_np1,1:steporder+1)
ϕ_coefs!(ϕstar_n,ϕ_np1,dy_np1,1:stageorder+1)

# Calculation of g_{steporder + 1}(n) coefficients using
# Calculation of g_{stageorder + 1}(n) coefficients using
# recurrence relations has already been done above, and does
# not depend on dy_np1 or y_np1, so no update necessary

# We correct the computation of y_np1 using y_np1
# as part of our prediction Hairer equation 5.7
y_np1 += dt*g[steporder+1]*ϕ[3,steporder+1]
y_np1 += dt*g[stageorder+1]*ϕ_np1[stageorder+1]

# (E)valuate estimate for next derivative again
dy_np1=F(t_np1,y_np1)
Expand All @@ -352,13 +349,13 @@ function ode113(F, y0::AbstractVector, tspan, order::Integer=4;
# ->if no, use this `next_dt`` for next attempt at this step
#---################################################################

#Don't allow increases of steps until n >= steporder + 2
#Don't allow increases of steps until n >= stageorder + 2
if n==1
successful_step = true
next_dt = dt/4
elseif n<steporder+2
yerr_k = dt*(g[steporder+1] - g[steporder])*ϕ[3,steporder+1]
err_k, next_dt,timeout_VS = stepsize_hw92!(dt, tdir, y[n], y_np1, yerr_k, steporder,
elseif n<stageorder+2
yerr_k = dt*(g[stageorder+1] - g[stageorder])*ϕ_np1[stageorder+1]
err_k, next_dt,timeout_VS = stepsize_hw92!(dt, tdir, y_n, y_np1, yerr_k, stageorder,
timeout_VS, dof, abstol, reltol, maxstep, norm)
if err_k <= 1
successful_step = true
Expand All @@ -372,13 +369,13 @@ function ode113(F, y0::AbstractVector, tspan, order::Integer=4;

# We need ϕ_{k+1}(n+1) and g_{k+1}(n);the latter has already
# been calculated above, the former still needs calculation
ϕ_coefs!(ϕ, ϕstar,dy_np1,steporder+2)
ϕ_coefs!(ϕstar_n,ϕ_np1,dy_np1,stageorder+2)

# Calculate ||LE_{k+1}|| and next stepsize

yerr_kp1= dt*(g[steporder+2] - g[steporder+1])*ϕ[3,steporder+2]
err_kp1, next_dt, timeout_VS = stepsize_hw92!(dt,tdir,y[n],y_np1,
yerr_kp1, steporder, timeout_VS,dof,
yerr_kp1= dt*(g[stageorder+2] - g[stageorder+1])*ϕ_np1[stageorder+2]
err_kp1, next_dt, timeout_VS = stepsize_hw92!(dt,tdir,y_n,y_np1,
yerr_kp1, stageorder, timeout_VS,dof,
abstol, reltol, maxstep, norm)

# Hairer equation 7.4: If ||LE_{k+1}|| =< 1
Expand Down Expand Up @@ -406,31 +403,37 @@ function ode113(F, y0::AbstractVector, tspan, order::Integer=4;
if tdir*(t_np1 + next_dt) >= tdir*(tfinal)
next_dt = tfinal-t_np1
end
push!(t,copy(t_np1))

#--#####################################################################
# Interpolation for dense output
#--#####################################################################

# interpolate specified points
# TODO: less expensive interpolation. Core problem is that current
# frame work does not allow for inexpensive interpolation for uneven
# step sizes.

while iter <= length(tspan) && tdir*tspan[iter]<tdir*t_np1
h_iter = tspan[iter]-t[n]
dt_iter = tspan[iter]-t[n]

# Calculation of g_{j}(n)coefficients for j=0,...,k+1;
g_coefs!(g,c,dt,t,tspan[iter],n,min(n,steporder+1))
g_coefs!(g,c,dt,t,tspan[iter],n,min(n,stageorder+1))

# Calculation of ϕstar_{j}(n), and ϕ_{j}(n) coef. for j=0,...,k;
ϕ_and_ϕstar_coefs!(b, ϕ, ϕstar,dy[n],t,tspan[iter],n,1:min(max(n-1,1),steporder+1))
ϕ_and_ϕstar_coefs!(b, ϕstar_nm1, ϕstar_n, ϕ_n,dy_n,t,tspan[iter],n,1:min(max(n-1,1),stageorder+1))

yout_iter = copy(y[n])
for j= 1:steporder
yout_iter += h_iter*g[j]*ϕstar[2,j]
yout_iter = copy(y_n)
for j= 1:stageorder
yout_iter += dt_iter*g[j]*ϕstar_n[j]
end

#CE iteration for more precise interpolation
if n>=steporder+1
# CE iteration for more precise interpolation
if n>=stageorder+1
dy_iter=F(tspan[iter],yout_iter)
ϕ_coefs!(ϕ, ϕstar,dy_iter,1:steporder+1)
ϕ_coefs!(ϕstar_n,ϕ_np1,dy_iter,1:stageorder+1)

yout_iter += h_iter*g[steporder+1]*ϕ[3,steporder+1]
yout_iter += dt_iter*g[stageorder+1]*ϕ_np1[stageorder+1]
end

push!(tout,tspan[iter])
Expand All @@ -449,57 +452,52 @@ function ode113(F, y0::AbstractVector, tspan, order::Integer=4;
end

# restore values of coefficients for choosen step size
g_coefs!(g,c,dt,t,t_np1,n,min(n,steporder+1))
ϕ_and_ϕstar_coefs!(b, ϕ, ϕstar,dy[n],t,t_np1,n,1:min(max(n-1,1),steporder+1))
ϕ_coefs!(ϕ, ϕstar,dy_np1,1:steporder+1)

# push work variables
push!(y,copy(y_np1))
push!(dy,copy(dy_np1))
push!(t,copy(t_np1))
g_coefs!(g,c,dt,t,t_np1,n,min(n,stageorder+1))
ϕ_and_ϕstar_coefs!(b, ϕstar_nm1, ϕstar_n, ϕ_n, dy_n,t,t_np1,n,1:min(max(n-1,1),stageorder+1))
ϕ_coefs!(ϕstar_n,ϕ_np1,dy_np1,1:stageorder+1)

#--#####################################################################
# Variable order selection for next step
# Note run until a success step was taken for n-th step
# Not run until a success step was taken for n-th step
# See Hairer pg 424 for indepeth description
#--#####################################################################

current_steporder = steporder
current_stageorder = stageorder
if n<=order+1
# Need to run the first few steps at reduced order for it requires
# the last 'order' number of points
steporder = n==1 ? 1 : n-1
stageorder = n==1 ? 1 : n-1

elseif adaptive_order && steporder>1
elseif adaptive_order && stageorder>1
# First, compute truncation errors for implicit orders: k-1,k
# against current k+1

yerr_km1 = dt*(g[steporder] - g[steporder-1])*ϕ[3,steporder]
err_km1 = stepsize_hw92!(dt, tdir, y[n], y_np1, yerr_km1, steporder,
yerr_km1 = dt*(g[stageorder] - g[stageorder-1])*ϕ_np1[stageorder]
err_km1 = stepsize_hw92!(dt, tdir, y_n, y_np1, yerr_km1, stageorder,
timeout_VS,dof, abstol, reltol, maxstep, norm)[1]

yerr_k = dt*(g[steporder+1] - g[steporder])*ϕ[3,steporder+1]
err_k = stepsize_hw92!(dt, tdir, y[n], y_np1, yerr_k, steporder,
yerr_k = dt*(g[stageorder+1] - g[stageorder])*ϕ_np1[stageorder+1]
err_k = stepsize_hw92!(dt, tdir, y_n, y_np1, yerr_k, stageorder,
timeout_VS, dof, abstol, reltol, maxstep, norm)[1]

# Decrement steporder if LE_{k} and LE_{k-1} < LE_{k+1}
# Decrement stageorder if LE_{k} and LE_{k-1} < LE_{k+1}
if max(err_km1,err_k) <= err_kp1
steporder = max(order,steporder - 1)
stageorder = max(order,stageorder - 1)

elseif consec_const_steps >= steporder && timeout_VO==0 && n>steporder + 3
elseif consec_const_steps >= stageorder && timeout_VO==0 && n>stageorder + 3
# Given enough consecutive steps taken at a constant step size
# we compute truncation error for order implicit k+2
ϕ_and_ϕstar_coefs!(b, ϕ, ϕstar,dy[n],t,t_np1,n,steporder+2)
ϕ_coefs!(ϕ, ϕstar,dy_np1,steporder+3)
ϕ_and_ϕstar_coefs!(b, ϕstar_nm1, ϕstar_n, ϕ_n,dy_n,t,t_np1,n,stageorder+2)
ϕ_coefs!(ϕstar_n,ϕ_np1,dy_np1,stageorder+3)

yerr_kp2 = dt*γ_star[steporder + 3]*ϕ[3,steporder+3]
err_kp2 = stepsize_hw92!(dt,tdir,y[n],y_np1,yerr_kp2,steporder,
yerr_kp2 = dt*γ_star[stageorder + 3]*ϕ_np1[stageorder+3]
err_kp2 = stepsize_hw92!(dt,tdir,y_n,y_np1,yerr_kp2,stageorder,
timeout_VS, dof, abstol, reltol, maxstep, norm)[1]

# Increment steporder if LE_{k+2} < LE_{k+1} and max(LE_{k},
# Increment stageorder if LE_{k+2} < LE_{k+1} and max(LE_{k},
# LE_{k-1}) > LE_{k+1}
if err_kp2 < err_kp1
steporder = min(max_order,steporder + 1)
stageorder = min(max_order,stageorder + 1)
timeout_VO = timeout_VO_const # reset timout_VO
end
end
Expand All @@ -513,15 +511,19 @@ function ode113(F, y0::AbstractVector, tspan, order::Integer=4;

# Shift data from this step, from the "current step" position to the
# "previous step" position for ϕstar and ϕ
# Note: ϕ[2,:] and ϕ_star[2,:] corresponds to ϕ_n and ϕ*_n, and
# Note: ϕ_n[:] and ϕ_star[2,:] corresponds to ϕ_n and ϕ*_n, and
# ϕ_star[1,:] to ϕ*_{n-1}, etc
# TODO: implement with DataStructures.CircularBuffer
ϕstar[1,1:min(max(n-1,1),current_steporder+1)] =
ϕstar[2,1:min(max(n-1,1),current_steporder+1)]
ϕstar_nm1[1:min(max(n-1,1),current_stageorder+1)] =
ϕstar_n[1:min(max(n-1,1),current_stageorder+1)]

# Update stepsize variables
prev_dt = dt
dt = next_dt

# push work variables
y_n = y_np1
dy_n = dy_np1

# Increment our step
n= n+1
end #while loop over steps
Expand Down Expand Up @@ -558,11 +560,11 @@ end
# Helper functions for calculating coefficients from recurence relationships
################################################################################
#Based on Hairer Lemma 5.1 on page 399
#Calculates c_j,q and g_n,j for n and j =1:steporder
#Calculates c_j,q and g_n,j for n and j =1:stageorder
"compute c_{j,q} and g_{n,j} coefficients for ode113"
function g_coefs!(g,c,dt,t,t_np1,n,steporder)
for j = 1:steporder
for q = 1:(steporder)-(j-1)
function g_coefs!(g,c,dt,t,t_np1,n,stageorder)
for j = 1:stageorder
for q = 1:((stageorder)-(j-1))
if j ==1
c[j,q]=1/q
elseif j==2
Expand All @@ -580,17 +582,16 @@ end
#Based on Hairer equation 5.9. Calculates b, ϕ, ϕstar for
#n and j in range_of_indices
"compute ϕ_{n,j} and ϕ*_{n,j} coefficients for ode113"
function ϕ_and_ϕstar_coefs!(b, ϕ, ϕstar,dy_n,t,t_np1,n,range_of_indices)
function ϕ_and_ϕstar_coefs!(b, ϕstar_nm1, ϕstar_n, ϕ_n,dy_n,t,t_np1,n,range_of_indices)
for j in range_of_indices #0:k-1
if j== 1#0
if j== 1
b[j] = 1
ϕ[2,j] = copy(dy_n)
ϕstar[2,j] = copy(dy_n)
else j>1#0
#Note: since n>=j, when this runs, n>=2 (so n-1>=1) and
ϕ_n[j] = copy(dy_n)
ϕstar_n[j] = copy(dy_n)
else
b[j] = b[j-1]*(t_np1-t[n-(j-1)+1])/(t[n]-t[n-(j-1)])
ϕ[2,j] = ϕ[2,j-1]- ϕstar[1,j-1]
ϕstar[2,j] = b[j]*ϕ[2,j]
ϕ_n[j] = ϕ_n[j-1]- ϕstar_nm1[j-1]
ϕstar_n[j] = b[j]*ϕ_n[j]
end
end
return nothing
Expand All @@ -599,13 +600,13 @@ end

#Based on Hairer equation 5.9. Just calculates the ϕ coefficients
"compute ϕ_{n+1,j} coefficients for ode113"
function ϕ_coefs!(ϕ, ϕstar,dy_n,range_of_indices)
function ϕ_coefs!(ϕstar_n,ϕ_np1,dy_n,range_of_indices)
for j in range_of_indices
if j== 1#0
ϕ[3,j] = copy(dy_n)
else j>1#0
ϕ_np1[j] = copy(dy_n)
else
#Note: since n>=j, when this runs, n>=2 (so n-1>=1) and
ϕ[3,j] = ϕ[3,j-1] - ϕstar[2,j-1]
ϕ_np1[j] = ϕ_np1[j-1] - ϕstar_n[j-1]
end
end
return nothing
Expand All @@ -615,11 +616,13 @@ end
# Helper functions initializing/resizing arrays
################################################################################
"used to initialize working arrays used in ode1113"
function intialize_coef_arrays(dim,T,Et,y0)
c = zeros(Et,dim,dim)
g = zeros(dim)
b = zeros(dim)
@compat ϕstar = Array{T,2}(2,dim)
@compat ϕ = Array{T,2}(3,dim)
return c,g,b,ϕstar, ϕ
function intialize_coef_arrays(max_order,T,Et)
c = zeros(Et,max_order,max_order)
g = zeros(max_order)
b = zeros(max_order)
@compat ϕstar_nm1 = Array{T,2}(1,max_order)
@compat ϕstar_n = Array{T,2}(1,max_order)
@compat ϕ_n = Array{T,2}(1,max_order)
@compat ϕ_np1 = Array{T,2}(1,max_order)
return c,g,b, ϕstar_nm1, ϕstar_n, ϕ_n, ϕ_np1
end

0 comments on commit 73b1ea9

Please sign in to comment.