Line data    Source code

       1             : struct RK43StepperState{N, T}
       2             :     # NOTE: nothing to store (solver stores Hermite piece)
       3         941 :     RK43StepperState{N, T}() where {N, T} = new()
       4             : end
       5             : 
       6           9 : Base.copy(::RK43StepperState{N, T}) where {N, T} = RK43StepperState{N, T}()
       7             : 
       8           2 : function Base.show(io::IO, ::RK43StepperState)
       9           2 :     print(io, "RK43")
      10             : end
      11             : 
      12             : # TODO: rename RK43Stepper, (discriminating) type: not NTuple (see TODO below)
      13             : 
      14     1283875 : function _step!(solver::Solver{N, T, RK43StepperState{N, T}}, dy,
      15             :                 t::T, h::T) where {N, T}
      16             :     # _step!(...) ::Maybe{SVector{N, T}} where ...
      17             : 
      18     1283875 :     fdy = guard(dy, solver)
      19             : 
      20             :     # NOTE: Do we have to store the state???
      21             : 
      22     1283875 :     y  = solver.y
      23     1283875 :     k0 = solver.dy * h
      24     2567750 :     k1 = fdy(t + h / 2, y + k0 / 2) * h
      25     2337168 :     k2 = fdy(t + h / 2, y + k1 / 2) * h
      26     2336269 :     k3 = fdy(t + h, y + k2) * h
      27             : 
      28     1688465 :     issuccess(fdy) || return Maybe{SVector{N, T}}(state(fdy))
      29             : 
      30      879285 :     solver.ynew = y + (k0 / 6 + k1 / 3 + k2 / 3 + k3 / 6)
      31     1758570 :     solver.dynew = fdy(t + h, solver.ynew)
      32             : 
      33      879285 :     issuccess(fdy) || return Maybe{SVector{N, T}}(state(fdy))
      34             : 
      35      879285 :     yerr = (k3 - solver.dynew * h) / 6
      36             : 
      37             :     # NOTE: Nothing to "commit" for RK43, we have all in y, dy, ynew, dynew
      38             :     # GENERALLY: solver.stepper = ...stepper_state...
      39             : 
      40      879285 :     Maybe{SVector{N, T}}(yerr, Ok)
      41             : end
      42             : 
      43           0 : order0(::RK43StepperState) = 4
      44           0 : order1(::RK43StepperState) = 4 # TODO: [alg] Is this correct for this tableau??
      45             : 
      46      213796 : function _interpolate(solver::Solver{N, T, RK43StepperState{N, T}},
      47             :                      t::T) where {N, T}
      48      213796 :     t0, t1 = solver.t, solver.tnew
      49             : 
      50      213796 :     (t == t0) && return (solver.y, solver.dy)
      51      211324 :     (t == t1) && return (solver.ynew, solver.dynew)
      52             : 
      53      207480 :     h = t1 - t0
      54             : 
      55      207480 :     t = (t - t0) / h
      56           0 :     @massert isfinite(t)
      57             : 
      58      207480 :     cy0, cdy0, cy1, cdy1 = HermiteSpline.hermitebasis(t)
      59      207480 :     y = solver.y*cy0 + solver.dy*(cdy0*h) +
      60             :         solver.ynew*cy1 + solver.dynew*(cdy1*h)
      61             : 
      62      207480 :     ey0, edy0, ey1, edy1 = HermiteSpline.hermitederivativebasis(t)
      63      207480 :     dy = solver.y*(ey0/h) + solver.dy*edy0 +
      64             :          solver.ynew*(ey1/h) + solver.dynew*edy1
      65             : 
      66      207480 :     y, dy
      67             : end