Line data    Source code

       1             : """
       2             :     h = guessinitailstepsize(solver, dy)
       3             : 
       4             : Determine initial step size (evaluates `dy` once).
       5             : 
       6             : !!! warning "precondition"
       7             :     Attributes `y, dy, t, tfinal, tdir` are initialized.
       8             : 
       9             : Does one additional Euler step (one evaluation of `dy`), see
      10             : *Hairer. Solving Ordinary Differential Equations I. p. 169*.
      11             : 
      12             : Code adapted from
      13             : [JuliaDiffEq:ODE.jl](https://github.com/JuliaDiffEq/ODE.jl/blob/master/src/ODE.jl)
      14             : """
      15       12261 : function guessinitialstepsize(solver::Solver{N, T, S}, dy) where {N, T, S}
      16       12261 :     @massert issuccess(solver.state)
      17             : 
      18       12261 :     t, tfinal, tdir = solver.t, solver.tfinal, solver.tdir
      19             : 
      20       12400 :     normy = norm(solver.y)
      21       12261 :     τ = max.(solver.opts.rtol * normy, solver.opts.atol)
      22             :     # Note: This is not exactly scaled_norm() (Hairer: same norm).
      23             :     #       Should be good enough.
      24             : 
      25       12261 :     d0 = normy / τ
      26       12761 :     d1 = norm(solver.dy) / τ
      27             : 
      28             : 
      29       24383 :     h0 =
      30             :         if ((d0 < T(1e-5)) || (d1 < T(1e-5)))
      31         514 :             T(1e-6)
      32             :         else
      33       24008 :             h0 = (d0 / d1) * T(0.01)
      34             :         end
      35             : 
      36             :     # perform Euler step
      37             : 
      38       12261 :     fdy = guard(dy, solver)
      39             : 
      40       12261 :     y1 = solver.y + tdir * h0 * solver.dy
      41       24522 :     dy1 = fdy(t + tdir * h0, y1)
      42             : 
      43       12261 :     issuccess(fdy) || return tdir * h0 / 2 # leave to step size control
      44             : 
      45             :     # estimate second derivative
      46             : 
      47       12733 :     d2 = norm(dy1 - solver.dy) / (τ * h0)
      48             : 
      49       12261 :     h1 =
      50             :         if (max(d1, d2) <= T(1e-15))
      51         466 :             max(T(1e-6), T(1e-3) * h0)
      52             :         else
      53       24056 :             (T(0.01) / max(d1, d2)) ^ (1 / order1(solver.stepper))
      54             :         end
      55             : 
      56       12261 :     return tdir * min(min(100 * h0, h1), tdir * (tfinal - t))
      57             : end
      58             : 
      59       12261 : function initialize!(solver::Solver{N, T, S},
      60             :                      dy, t0::T,  t1::T, y0::SVector{N, T}) where {N, T, S}
      61       12261 :     solver.state = Ok
      62             : 
      63       12261 :     (t0 != t1) || @warn "initialized solver with t0=t1=$(t0)"
      64             : 
      65       12261 :     solver.t = t0
      66       12261 :     solver.tfinal = t1
      67       12261 :     solver.tdir = (t1 - t0) >= 0 ? +1 : -1
      68       12261 :     solver.y = y0
      69       12261 :     solver.nsteps = 0
      70       12261 :     solver.rejected = false
      71             : 
      72       12261 :     fdy = guard(dy, solver)
      73       24522 :     solver.dy = fdy(t0, y0)
      74             : 
      75       12261 :     if !issuccess(fdy)
      76           0 :         @warn "initialization failed: could not evaluate dy($(t0),$(y0))"
      77           0 :         return (solver.state = Failed)
      78             :     end
      79             : 
      80       12261 :     solver.ynew = solver.y
      81       12261 :     solver.dynew = solver.dy
      82       12261 :     solver.tnew = solver.t
      83             : 
      84       12261 :     hinit = solver.opts.h0
      85       24522 :     (hinit > 0) || (hinit = guessinitialstepsize(solver, dy))
      86             : 
      87       12261 :     solver.absh = abs(hinit)
      88           0 :     Ok
      89             : end