DP8 dense output about differentialequations.jl

[ChrisRackauckas]

Details:

The Bogacki-Shampine pair is documented by RKSuite where the coefficients taken from, following the paper. The code shows how to do the interpolation, but this ported "overshoots" the edge leading to little gaps.

Edit

Fixed. Dr. Bogacki sent me different coefficients which worked perfectly.

[ChrisRackauckas]

The DP8 code is derived from Hairer's dop853. The dense output also overshoots. The relevant code from there is:

               CONT(J)=Y(I)
               YDIFF=K5(I)-Y(I)
               CONT(J+NRD)=YDIFF
               BSPL=H*K1(I)-YDIFF
               CONT(J+NRD*2)=BSPL
               CONT(J+NRD*3)=YDIFF-H*K4(I)-BSPL
               CONT(J+NRD*4)=D41*K1(I)+D46*K6(I)+D47*K7(I)+D48*K8(I)
     &                  +D49*K9(I)+D410*K10(I)+D411*K2(I)+D412*K3(I)
               CONT(J+NRD*5)=D51*K1(I)+D56*K6(I)+D57*K7(I)+D58*K8(I)
     &                  +D59*K9(I)+D510*K10(I)+D511*K2(I)+D512*K3(I)
               CONT(J+NRD*6)=D61*K1(I)+D66*K6(I)+D67*K7(I)+D68*K8(I)
     &                  +D69*K9(I)+D610*K10(I)+D611*K2(I)+D612*K3(I)
               CONT(J+NRD*7)=D71*K1(I)+D76*K6(I)+D77*K7(I)+D78*K8(I)
     &                  +D79*K9(I)+D710*K10(I)+D711*K2(I)+D712*K3(I)
   62       CONTINUE 
C ---     THE NEXT THREE FUNCTION EVALUATIONS
            DO 51 I=1,N 
  51           Y1(I)=Y(I)+H*(A141*K1(I)+A147*K7(I)+A148*K8(I)
     &            +A149*K9(I)+A1410*K10(I)+A1411*K2(I)+A1412*K3(I)
     &            +A1413*K4(I))
            CALL FCN(N,X+C14*H,Y1,K10,RPAR,IPAR)
            DO 52 I=1,N 
  52           Y1(I)=Y(I)+H*(A151*K1(I)+A156*K6(I)+A157*K7(I)
     &            +A158*K8(I)+A1511*K2(I)+A1512*K3(I)+A1513*K4(I)
     &            +A1514*K10(I))
            CALL FCN(N,X+C15*H,Y1,K2,RPAR,IPAR)
            DO 53 I=1,N 
  53           Y1(I)=Y(I)+H*(A161*K1(I)+A166*K6(I)+A167*K7(I)
     &            +A168*K8(I)+A169*K9(I)+A1613*K4(I)+A1614*K10(I)
     &            +A1615*K2(I))
            CALL FCN(N,X+C16*H,Y1,K3,RPAR,IPAR)
            NFCN=NFCN+3 
C ---     FINAL PREPARATION
            DO 63 J=1,NRD
               I=ICOMP(J)
               CONT(J+NRD*4)=H*(CONT(J+NRD*4)+D413*K4(I)+D414*K10(I)
     &            +D415*K2(I)+D416*K3(I))
               CONT(J+NRD*5)=H*(CONT(J+NRD*5)+D513*K4(I)+D514*K10(I)
     &            +D515*K2(I)+D516*K3(I))
               CONT(J+NRD*6)=H*(CONT(J+NRD*6)+D613*K4(I)+D614*K10(I)
     &            +D615*K2(I)+D616*K3(I))
               CONT(J+NRD*7)=H*(CONT(J+NRD*7)+D713*K4(I)+D714*K10(I)
     &            +D715*K2(I)+D716*K3(I))

However, this makes no sense because earlier the saving is

  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)
      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)
  35  K5(I)=Y(I)+H*K4(I)

which means K5(I)=Y(I)+H*K4(I), so YDIFF=K5(I)-Y(I) is simply H*K4 which is dt times k13, or the rate of the 13th stage which is the combination previous stages used for the update. Then BSPL=H*K1(I)-YDIFF or k1 - k13 if written as rates. But then this means that the next thing saved is BSPL = YDIFF-H*K4(I)-BSPL which when written as rates and substituted is just -BSPL since YDIFF-H*K4(I) = 0? There is something weird here if it's written like this.

[ChrisRackauckas]

Update on BS5: I discussed this directly with Peter Bogacki and Lawrence Shampine and got updated coefficients which work (also got lots of ideas for new algorithms/APIs/interfaces). Still working on DP8

[ChrisRackauckas]

K4 in the dop853 is updated as f(t,utmp) and so the above simplification was wrong. This was fixed and the tests now pass for all dense!