Implementation of DP8.

src/dp5/helpers.jl

@@ -1,37 +1,5 @@
 function do_step!(solver, h)
 
-    # define constants
-    c2=0.2
-    c3=0.3
-    c4=0.8
-    c5=8.0/9.0
-    a21=0.2
-    a31=3.0/40.0
-    a32=9.0/40.0
-    a41=44.0/45.0
-    a42=-56.0/15.0
-    a43=32.0/9.0
-    a51=19372.0/6561.0
-    a52=-25360.0/2187.0
-    a53=64448.0/6561.0
-    a54=-212.0/729.0
-    a61=9017.0/3168.0
-    a62=-355.0/33.0
-    a63=46732.0/5247.0
-    a64=49.0/176.0
-    a65=-5103.0/18656.0
-    a71=35.0/384.0
-    a73=500.0/1113.0
-    a74=125.0/192.0
-    a75=-2187.0/6784.0
-    a76=11.0/84.0
-    e1=71.0/57600.0
-    e3=-71.0/16695.0
-    e4=71.0/1920.0
-    e5=-17253.0/339200.0
-    e6=22.0/525.0
-    e7=-1.0/40.0
-
     ####### First 6 stages
 
     solver.y1 .= solver.y .+ h .* a21 .* solver.k1

src/dp5/mod.jl

@@ -2,6 +2,7 @@ module  DP5
 
 using ..DormandPrince: DormandPrince, DP5Solver, Vars, Consts, Options, Report
 
+include("params.jl")
 include("helpers.jl")
 include("checks.jl")
 include("solver.jl")

src/dp5/params.jl

@@ -0,0 +1,31 @@
+# define constants
+const c2=0.2
+const c3=0.3
+const c4=0.8
+const c5=8.0/9.0
+const a21=0.2
+const a31=3.0/40.0
+const a32=9.0/40.0
+const a41=44.0/45.0
+const a42=-56.0/15.0
+const a43=32.0/9.0
+const a51=19372.0/6561.0
+const a52=-25360.0/2187.0
+const a53=64448.0/6561.0
+const a54=-212.0/729.0
+const a61=9017.0/3168.0
+const a62=-355.0/33.0
+const a63=46732.0/5247.0
+const a64=49.0/176.0
+const a65=-5103.0/18656.0
+const a71=35.0/384.0
+const a73=500.0/1113.0
+const a74=125.0/192.0
+const a75=-2187.0/6784.0
+const a76=11.0/84.0
+const e1=71.0/57600.0
+const e3=-71.0/16695.0
+const e4=71.0/1920.0
+const e5=-17253.0/339200.0
+const e6=22.0/525.0
+const e7=-1.0/40.0

src/dp8/checks.jl

@@ -0,0 +1,34 @@
+# include("types.jl")
+
+# make enums for every error type and return that instead of printing errors
+function check_max_allowed_steps(options::Options{T}) where T
+    if options.maximum_allowed_steps < 0
+        return false
+    else
+        return true
+    end
+end
+
+function check_uround(options::Options{T}) where T
+    if (options.uround <= 1e-35) || (options.uround >= 1.0) 
+        return false
+    else
+        return true
+    end
+end
+
+function check_beta(options::Options{T}) where T
+    if options.beta > 0.2
+        return false
+    else
+        return true
+    end
+end
+
+function check_safety_factor(options::Options{T}) where T
+    if (options.safety_factor >= 1.0) || (options.safety_factor <= 1e-4)
+        return false
+    else
+        return true
+    end
+end

src/dp8/helpers.jl

@@ -0,0 +1,211 @@
+function do_step!(solver, h)
+
+    ####### The 12 stages
+    # DO 22 I=1,N 
+    # 22  Y1(I)=Y(I)+H*A21*K1(I)  
+    #     CALL FCN(N,X+C2*H,Y1,K2,RPAR,IPAR)
+    solver.y1 .= solver.y .+ h .* a21 .* solver.k1
+    solver.f(solver.vars.x + c2 * h, solver.y1, solver.k2)
+    #     DO 23 I=1,N 
+    # 23  Y1(I)=Y(I)+H*(A31*K1(I)+A32*K2(I))  
+    #     CALL FCN(N,X+C3*H,Y1,K3,RPAR,IPAR)
+    solver.y1 .= solver.y .+ h .* (a31 .* solver.k1 .+ a32 .* solver.k2)
+    solver.f(solver.vars.x + c3 * h, solver.y1, solver.k3)
+    #     DO 24 I=1,N 
+    # 24  Y1(I)=Y(I)+H*(A41*K1(I)+A43*K3(I))  
+    #     CALL FCN(N,X+C4*H,Y1,K4,RPAR,IPAR)
+    solver.y1 .= solver.y .+ h .* (a41 .* solver.k1 .+ a43 .* solver.k3)
+    solver.f(solver.vars.x + c4 * h, solver.y1, solver.k4)
+    #     DO 25 I=1,N 
+    # 25  Y1(I)=Y(I)+H*(A51*K1(I)+A53*K3(I)+A54*K4(I))
+    #     CALL FCN(N,X+C5*H,Y1,K5,RPAR,IPAR)
+    solver.y1 .= solver.y .+ h .* (a51 .* solver.k1 .+ a53 .* solver.k3 .+ a54 .* solver.k4)
+    solver.f(solver.vars.x + c5*h, solver.y1, solver.k5)
+    #     DO 26 I=1,N 
+    # 26  Y1(I)=Y(I)+H*(A61*K1(I)+A64*K4(I)+A65*K5(I))
+    #     CALL FCN(N,X+C6*H,Y1,K6,RPAR,IPAR)
+    solver.y1 .= solver.y .+ h .* (a61 .* solver.k1 .+ a64 .* solver.k4 .+ a65 .* solver.k5)
+    solver.f(solver.vars.x + c6*h, solver.y1, solver.k6)
+    #     DO 27 I=1,N 
+    # 27  Y1(I)=Y(I)+H*(A71*K1(I)+A74*K4(I)+A75*K5(I)+A76*K6(I))
+    #     CALL FCN(N,X+C7*H,Y1,K7,RPAR,IPAR)
+    solver.y1 .= solver.y .+ h .* (
+        a71 .* solver.k1 .+ a74 .* solver.k4 .+ a75 .* solver.k5 .+ a76 .* solver.k6
+    )
+    solver.f(solver.vars.x + c7*h, solver.y1, solver.k7)
+    #     DO 28 I=1,N 
+    # 28  Y1(I)=Y(I)+H*(A81*K1(I)+A84*K4(I)+A85*K5(I)+A86*K6(I)+A87*K7(I))  
+    #     CALL FCN(N,X+C8*H,Y1,K8,RPAR,IPAR)
+    solver.y1 .= solver.y .+ h .* (
+        a81 .* solver.k1 .+ a84 .* solver.k4 .+ a85 .* solver.k5 .+ a86 .* solver.k6 .+ a87 .* solver.k7
+    )
+    solver.f(solver.vars.x + c8*h, solver.y1, solver.k8)
+    #     DO 29 I=1,N 
+    # 29  Y1(I)=Y(I)+H*(A91*K1(I)+A94*K4(I)+A95*K5(I)+A96*K6(I)+A97*K7(I)
+    #    &   +A98*K8(I))
+    #     CALL FCN(N,X+C9*H,Y1,K9,RPAR,IPAR)
+    solver.y1 .= solver.y .+ h .* (
+        a91 .* solver.k1 .+ a94 .* solver.k4 .+ a95 .* solver.k5 .+ 
+        a96 .* solver.k6 .+ a97 .* solver.k7 .+ a98 .* solver.k8
+    )
+    solver.f(solver.vars.x + c9*h, solver.y1, solver.k9)
+    #     DO 30 I=1,N 
+    # 30  Y1(I)=Y(I)+H*(A101*K1(I)+A104*K4(I)+A105*K5(I)+A106*K6(I)
+    #    &   +A107*K7(I)+A108*K8(I)+A109*K9(I))
+    #     CALL FCN(N,X+C10*H,Y1,K10,RPAR,IPAR)
+    solver.y1 .= solver.y .+ h .* (
+        a101 .* solver.k1 .+ a104 .* solver.k4 .+ a105 .* solver.k5 .+ 
+        a106 .* solver.k6 .+ a107 .* solver.k7 .+ a108 .* solver.k8 .+ a109 .* solver.k9
+    )
+    solver.f(solver.vars.x + c10*h, solver.y1, solver.k10)
+    #     DO 31 I=1,N 
+    # 31  Y1(I)=Y(I)+H*(A111*K1(I)+A114*K4(I)+A115*K5(I)+A116*K6(I)
+    #    &   +A117*K7(I)+A118*K8(I)+A119*K9(I)+A1110*K10(I))
+    #     CALL FCN(N,X+C11*H,Y1,K2,RPAR,IPAR)
+    solver.y1 .= solver.y .+ h .* (
+        a111 .* solver.k1 .+ a114 .* solver.k4 .+ a115 .* solver.k5 .+ 
+        a116 .* solver.k6 .+ a117 .* solver.k7 .+ a118 .* solver.k8 .+ 
+        a119 .* solver.k9 .+ a1110 .* solver.k10
+    )
+    solver.f(solver.vars.x + c11*h, solver.y1, solver.k2)
+    #     XPH=X+H
+    xph = solver.vars.x + h
+    #     DO 32 I=1,N 
+    # 32  Y1(I)=Y(I)+H*(A121*K1(I)+A124*K4(I)+A125*K5(I)+A126*K6(I)
+    #    &   +A127*K7(I)+A128*K8(I)+A129*K9(I)+A1210*K10(I)+A1211*K2(I))
+    #     CALL FCN(N,XPH,Y1,K3,RPAR,IPAR)
+    solver.y1 .= solver.y .+ h .* (
+        a121 .* solver.k1 .+ a124 .* solver.k4 .+ a125 .* solver.k5 .+ 
+        a126 .* solver.k6 .+ a127 .* solver.k7 .+ a128 .* solver.k8 .+ 
+        a129 .* solver.k9 .+ a1210 .* solver.k10 .+ a1211 .* solver.k2
+    )
+    solver.f(xph, solver.y1, solver.k3)
+    #     NFCN=NFCN+11
+    #     DO 35 I=1,N 
+    #     K4(I)=B1*K1(I)+B6*K6(I)+B7*K7(I)+B8*K8(I)+B9*K9(I)
+    #    &   +B10*K10(I)+B11*K2(I)+B12*K3(I)   
+    #     K5(I)=Y(I)+H*K4(I)
+    solver.k4 .= (
+        b1 .* solver.k1 .+ b6 .* solver.k6 .+ b7 .* solver.k7 .+ b8 .* solver.k8 .+ 
+        b9 .* solver.k9 .+ b10 .* solver.k10 .+ b11 .* solver.k2 .+ b12 .* solver.k3
+    )
+    solver.k5 .= solver.y .+ h .* solver.k4
+
+end
+
+function error_estimation(solver)
+
+    err, err2 = mapreduce(+, 
+        solver.consts.atol_iter, 
+        solver.consts.rtol_iter, 
+        solver.y,
+        solver.k1,
+        solver.k3,
+        solver.k4,
+        solver.k5,
+        solver.k6,
+        solver.k7,
+        solver.k8,
+        solver.k9,
+        solver.k10,
+    ) do atoli, rtoli, yi, k1i, k2i, k3i, k4i, k5i, k6i, k7i, k8i, k9i, k10i
+
+        #     DO 42 I=1,N 
+        #     SK=ATOL(I)+RTOL(I)*MAX(ABS(Y(I)),ABS(K5(I)))
+        #     ERRI=K4(I)-BHH1*K1(I)-BHH2*K9(I)-BHH3*K3(I)
+        #     ERR2=ERR2+(ERRI/SK)**2
+        #     ERRI=ER1*K1(I)+ER6*K6(I)+ER7*K7(I)+ER8*K8(I)+ER9*K9(I)
+        # &      +ER10*K10(I)+ER11*K2(I)+ER12*K3(I)
+        # 42    ERR=ERR+(ERRI/SK)**2
+
+        sk = atoli + rtoli*max(abs(yi), abs(k5i))
+        erri = k4i - bhh1*k1i - bhh2*k9i - bhh3*k3i
+        err2 = (abs(erri)/sk)^2
+
+        erri = er1*k1i + er6*k6i + er7*k7i + er8*k8i + er9*k9i + er10*k10i + er11*k2i + er12*k3i
+        err = (abs(erri)/sk)^2
+
+        (err, err2)
+
+    end
+
+    # DENO=ERR+0.01D0*ERR2
+    # IF (DENO.LE.0.D0) DENO=1.D0
+    # ERR=ABS(H)*ERR*SQRT(1.D0/(N*DENO))
+
+    den0 = err + 0.01*err2
+    den0 = if den0 <= 0.0
+        1.0
+    else
+        den0
+    end
+
+    err = abs(solver.vars.h)*err*sqrt(1.0/(length(solver.y) * den0))
+
+    return err 
+end
+
+function estimate_second_derivative(solver, h)
+
+    der2 = mapreduce(+, solver.consts.atol_iter, solver.consts.rtol_iter, solver.k2, solver.k1, solver.y) do atoli, rtoli, f1i, f0i, yi
+        sk = atoli + rtoli*abs(yi)
+        ((f1i-f0i)/sk)^2
+    end
+
+    der2 = sqrt(der2)/h
+
+    return der2
+
+end
+
+function stiffness_detection!(solver, naccpt, h)
+    if (mod(naccpt, solver.options.stiffness_test_activation_step) == 0) || (solver.vars.iasti > 0)
+        #stnum = 0.0
+        #stden = 0.0
+
+        stnum, stden = mapreduce(.+, solver.k2, solver.k6, solver.y1, solver.ysti) do k2i, k6i, y1i, ystii
+            #stnum = abs(k2i-k6i)^2 
+            #stden = abs(y1i-ystii)^2
+            abs(k2i-k6i)^2, abs(y1i-ystii)^2
+            # stnum, stden
+        end
+
+        if stden > 0.0
+            solver.vars.hlamb = h*sqrt(stnum/stden)
+        else
+            solver.vars.hlamb = Inf
+        end
+
+
+        if solver.vars.hlamb > 3.25
+            solver.vars.iasti += 1
+            if solver.vars.iasti == 15
+                @debug "The problem seems to become stiff at $x" 
+            end
+        else 
+            solver.vars.nonsti += 1
+            if solver.vars.nonsti == 6
+                solver.vars.iasti = 0
+            end
+        end
+    end
+end
+
+function euler_first_guess(solver, hmax, posneg)
+
+    dnf, dny = mapreduce(.+, solver.consts.atol_iter, solver.consts.rtol_iter, solver.k1, solver.y) do atoli, rtoli, f0i, yi
+        sk = atoli + rtoli*abs(yi)
+        abs(f0i/sk)^2, abs(yi/sk)^2 # dnf, dny
+    end
+
+
+    if (dnf <= 1.0e-10) || (dny <= 1.0e-10)
+        h = 1.0e-6
+    else
+        h = 0.01*sqrt(dny/dnf)
+    end
+    h = min(h, hmax)
+    h = h * Base.sign(posneg)
+
+    return h, dnf
+end

src/dp8/mod.jl

@@ -0,0 +1,10 @@
+module  DP8
+
+using ..DormandPrince: DormandPrince, DP8Solver, Vars, Consts, Options, Report
+
+include("params.jl")
+include("helpers.jl")
+include("checks.jl")
+include("solver.jl")
+
+end