Verner 9 Runge Kutta Interpolator Issue

New to Julia
1

Hello, I am trying to implement the vern9 runge kutta method for a side project, and so far, everything has been good, except for the interpolation of the states for a dense output. For some reason, when compared with a true data, the interpolated values are off by quite a lot. I am wondering if my understanding of how to implement the interpolator is just wrong then. It seems straightforward from the source code online, and also from the vern 1993 paper as well. Any help would be greatly appreciated!

def vern9pre0(theta, r):

    theta2 = theta * theta
    
    b1  = theta  * evalpoly(theta, r["r011"], r["r012"], r["r013"], r["r014"], r["r015"], r["r016"], r["r017"], r["r018"], r["r019"])
    b8  = theta2 * evalpoly(theta, r["r082"], r["r083"], r["r084"], r["r085"], r["r086"], r["r087"], r["r088"], r["r089"])
    b9  = theta2 * evalpoly(theta, r["r092"], r["r093"], r["r094"], r["r095"], r["r096"], r["r097"], r["r098"], r["r099"])
    b10 = theta2 * evalpoly(theta, r["r102"], r["r103"], r["r104"], r["r105"], r["r106"], r["r107"], r["r108"], r["r109"])
    b11 = theta2 * evalpoly(theta, r["r112"], r["r113"], r["r114"], r["r115"], r["r116"], r["r117"], r["r118"], r["r119"])
    b12 = theta2 * evalpoly(theta, r["r122"], r["r123"], r["r124"], r["r125"], r["r126"], r["r127"], r["r128"], r["r129"])
    b13 = theta2 * evalpoly(theta, r["r132"], r["r133"], r["r134"], r["r135"], r["r136"], r["r137"], r["r138"], r["r139"])
    b14 = theta2 * evalpoly(theta, r["r142"], r["r143"], r["r144"], r["r145"], r["r146"], r["r147"], r["r148"], r["r149"])
    b15 = theta2 * evalpoly(theta, r["r152"], r["r153"], r["r154"], r["r155"], r["r156"], r["r157"], r["r158"], r["r159"])
    b17 = theta2 * evalpoly(theta, r["r172"], r["r173"], r["r174"], r["r175"], r["r176"], r["r177"], r["r178"], r["r179"])
    b18 = theta2 * evalpoly(theta, r["r182"], r["r183"], r["r184"], r["r185"], r["r186"], r["r187"], r["r188"], r["r189"])
    b19 = theta2 * evalpoly(theta, r["r192"], r["r193"], r["r194"], r["r195"], r["r196"], r["r197"], r["r198"], r["r199"])
    b20 = theta2 * evalpoly(theta, r["r202"], r["r203"], r["r204"], r["r205"], r["r206"], r["r207"], r["r208"], r["r209"])
    b21 = theta2 * evalpoly(theta, r["r212"], r["r213"], r["r214"], r["r215"], r["r216"], r["r217"], r["r218"], r["r219"])
    b22 = theta2 * evalpoly(theta, r["r222"], r["r223"], r["r224"], r["r225"], r["r226"], r["r227"], r["r228"], r["r229"])
    b23 = theta2 * evalpoly(theta, r["r232"], r["r233"], r["r234"], r["r235"], r["r236"], r["r237"], r["r238"], r["r239"])
    b24 = theta2 * evalpoly(theta, r["r242"], r["r243"], r["r244"], r["r245"], r["r246"], r["r247"], r["r248"], r["r249"])
    b25 = theta2 * evalpoly(theta, r["r252"], r["r253"], r["r254"], r["r255"], r["r256"], r["r257"], r["r258"], r["r259"])
    b26 = theta2 * evalpoly(theta, r["r262"], r["r263"], r["r264"], r["r265"], r["r266"], r["r267"], r["r268"], r["r269"])
    
    return (b1, b8, b9, b10, b11, b12, b13, b14, b15,
            b17, b18, b19, b20, b21, b22, b23, b24, b25, b26)

def ode_interpolant(theta, dt, y0, y1, k, r):
    (b1, b8, b9, b10, b11, b12, b13, b14, b15,
     b17, b18, b19, b20, b21, b22, b23, b24, b25, b26) = vern9pre0(theta, r)

    return y0 + dt * (
        k[0]  * b1  + k[1]  * b8  + k[2]  * b9  + k[3]  * b10 +
        k[4]  * b11 + k[5]  * b12 + k[6]  * b13 + k[7]  * b14 +
        k[8]  * b15 + k[10] * b17 + k[11] * b18 + k[12] * b19 +
        k[13] * b20 + k[14] * b21 + k[15] * b22 + k[16] * b23 +
        k[17] * b24 + k[18] * b25 + k[19] * b26
    )
2

Are you sure those are the right k’s? It won’t be the same k’s as the integration.

3

Yeah, the k values are mapped according to what’s in the source code. I’m pretty sure the mistake is on my end, but I wanted to be sure of my understanding of the theory.

4

So yeah if you mapped the ks correctly to their right indexes then this is at least right in theory, sans any typos I didn’t catch.

1 Like