diff --git a/src/adams_methods.jl b/src/adams_methods.jl index a5593fc7..1a727041 100644 --- a/src/adams_methods.jl +++ b/src/adams_methods.jl @@ -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] @@ -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 @@ -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 @@ -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) @@ -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= 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]=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]) @@ -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 @@ -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 @@ -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 @@ -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 @@ -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 @@ -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