Fixed bug in runge-kutta

When intital time step overstepped tspan[end] then it didn't work correctly.
mauro3 · Jun 10, 2015 · 4c5683a · 4c5683a
1 parent 037d19c
commit 4c5683a
Show file tree

Hide file tree

Showing 2 changed files with 7 additions and 6 deletions.
diff --git a/src/ODE.jl b/src/ODE.jl
@@ -62,6 +62,7 @@ Base.convert{Tnew<:Real}(::Type{Tnew}, tab::Tableau) = error("Define convert met
 # estimator for initial step based on book
 # "Solving Ordinary Differential Equations I" by Hairer et al., p.169
 function hinit(F, x0, t0, tend, p, reltol, abstol)
+    # 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)
@@ -84,7 +85,7 @@ function hinit(F, x0, t0, tend, p, reltol, abstol)
         pow = -(2. + log10(max(d1, d2)))/(p + 1.)
         h1 = 10.^pow
     end
-    return tdir*min(100*h0, h1), tdir, f0
+    return tdir*min(100*h0, h1, tdir*(tend-t0)), tdir, f0
 end
 
 # isoutofdomain takes the state and returns true if state is outside

diff --git a/src/runge_kutta.jl b/src/runge_kutta.jl
@@ -286,7 +286,7 @@ function oderk_adapt{N,S}(fn, y0::AbstractVector, tspan, btab_::TableauRKExplici
     steps = [0,0]  # [accepted, rejected]
 
     ## Integration loop
-    laststep = 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,7 +307,7 @@ function oderk_adapt{N,S}(fn, y0::AbstractVector, tspan, btab_::TableauRKExplici
             f1 = isFSAL(btab) ? ks[S] : fn(t+dt, ytrial)
             if points==:specified
                 # interpolate onto given output points
-                while iter-1<nsteps_fixed && (tdir*tspan[iter]<tdir*(t+dt) || laststep) # output at all new times which are < t+dt
+                while iter-1<nsteps_fixed && (tdir*tspan[iter]<tdir*(t+dt) || islaststep) # output at all new times which are < t+dt
                     hermite_interp!(ys[iter], tspan[iter], t, dt, y, ytrial, f0, f1) # TODO: 3rd order only!
                     iter += 1
                 end
@@ -328,7 +328,7 @@ function oderk_adapt{N,S}(fn, y0::AbstractVector, tspan, btab_::TableauRKExplici
             ks[1] = f1 # load ks[1]==f0 for next step
 
             # Break if this was the last step:
-            laststep && break
+            islaststep && break
 
             # Swap bindings of y and ytrial, avoids one copy
             y, ytrial = ytrial, y
@@ -340,13 +340,13 @@ function oderk_adapt{N,S}(fn, y0::AbstractVector, tspan, btab_::TableauRKExplici
             # Hit end point exactly if next step within 1% of end:
             if tdir*(t+dt*1.01) >= tdir*tend
                 dt = tend-t
-                laststep = true # next step is the last, if it succeeds
+                islaststep = true # next step is the last, if it succeeds
             end
         elseif abs(newdt)<minstep  # minimum step size reached, break
             println("Warning: dt < minstep.  Stopping.")
             break
         else # redo step with smaller dt
-            laststep = false
+            islaststep = false
             steps[2] +=1
             dt = newdt
             timeout = timeout_const