Butcher's Tableau to Diagonally Implicit Runge Kutta method

Question

I'm trying to implement a function that would input a butcher's tableau for a Runge-Kutta method of order $s$, an ode and its jacobian, as well as the initial values and return the next step. If we have a butcher tableau, $$ \begin{array}{c|cc} \alpha & \beta \\ \hline & \gamma \\ \end{array} $$ we'd like to convert that into a diagonally implicit butcher tableau $$ y_{n+1} = y_n + \sum_{i=1}^s \gamma_i \kappa_i $$ with \begin{align} \kappa_1 &= f\left(t_0 + \alpha_1 h, y_0 + \beta_{1j}\kappa_j\right) \\ \kappa_2 &= f\left(t_0 + \alpha_2 h, y_0 + h\sum_{j=1}^2 \beta_{2j}\kappa_j\right)\\ &\vdots \\ \kappa_s &= f\left(t_0 + \alpha_s h, y_0 + h\sum_{j=1}^s \beta_{sj}\kappa_j\right). \end{align} Keep in mind for the i-th step we need to solve for $\kappa_i$, where we use Newton Method to solve for a vector $\delta_i$ s.t. $F(\delta_i) = 0$, $$ F(\delta_i) = \delta_i - f\left(t_0 + \alpha_i h, y_0 + \beta_{ii}\delta_i + h\sum_{j=1}^{i-1} \beta_{ij}\kappa_j\right) $$ where DF is our Jacobian matrix.

def dirk(f, Df, t0, y0, h, alpha, beta, gamma):
    #---# Parameters #-------------------------------------------#
    #f(t,y) : lambda function of the ODE y'(t) = f(t, y)
    #Df(t,y): lambda function for the Jacobian of f
    #t0     : initial time
    #y0     : initial coordinate
    #alpha  : Butcher tableau vector of size s (weights of y step)
    #beta   : Butcher tableau matrix of size s*s (weights of steps)
    #gamma  : Butcher tableau vector of size s (weights of k)
    #------------------------------------------------------------#

    n = len(beta)
    m = len(f(t0, y0))
    K = np.zeros(len(alpha))
    if beta[0][0] != 0:
        K0,_ = newton(lambda delta: delta - f(t0 + alpha[0]*h, y0+h*beta[0][0]*delta), lambda delta: np.eye(len(y0))-Df(t0+alpha[0]*h,y0+h*beta[0][0]*delta), f(t0, y0), 1e-8,1000) 
    else:
        K0 = y0 + h*f(t0, y0)
    K = np.insert(K, 0, K0)
    for i in range(1, n):
        sum = 0
        for j in range(0, i):
            sum += beta[i][j]*K[i-1]
        if beta[i][i] != 0:
            F =lambda delta: delta - f(t0+alpha[i]*h, y0+h*(delta*beta[i][i]+sum))
            DF = lambda delta: np.eye(m)-Df(t0+alpha[i]*h,y0+h*(beta[i][i]*delta+sum))
            sol,_ = newton(F, DF, f(t0, y0), 1e-8,1000)
            print(sol)
        else:
            sol = f(t0 + alpha[i]*h, y0+h*sum)
        
        K = np.insert(K, i, sol)
    phi = 0
    for i in range(n):
        phi += gamma[i]*K[i]
    return y0 + h*phi

def newton(F,DF,x0,eps,K):
    #---# Parameters #-------------------------------------------#
    #F   : lambda function F(delta)
    #DF  : our lambda function of the Jacobian of F
    #x0  : our initial point f(t0, y0)
    #eps : our tolerance for 0
    #K   : maximum number of steps
    #------------------------------------------------------------# 
  
    k  = 0
    x  = x0.copy().astype(float) 
    delta = np.zeros([len(x0)])
    Fx = F(x)
    while Fx.dot(Fx) > eps*eps and k<K:
        delta[:] = np.linalg.solve(DF(x), Fx)
        x[:] -= delta[:]  
        Fx = F(x)
        k += 1
    return [x,k]

We can also use this function evolve to find the value for a given $T>t_0$ over $N$ steps so $h=(T-t_0)/N$.

def evolve(phi, f,Df, t0,y0, T,N):
    #---# Parameters #-------------------------------------------#
    #phi: our method
    #f(t, y): lambda function
    #Df(t, y): Jacobian of our function
    #t0: starting time
    #y0: initial coordinates
    #T: final time
    #N: number of steps
    #------------------------------------------------------------# 

    h = T/N
    y = np.zeros( [N+1,len(y0)] )
    y[0] = y0
    t = 0
    for i in range(N):
        y[i+1] = phi(f,Df, t,y[i], h)
        t = t+h
    return y

I'm testing on the code

stepper = lambda f,Df,t0,y0,h: dirk(f,Df,t0,y0,h,aCN,bCN,gCN)
y = evolve(stepper,
           lambda t,y: array([ -y[0],
                                y[0]-exp(-y[2]),
                                1 ]),
           lambda t,y: array([[-1,0,0],[1,0,exp(-y[2])],[0,0,0]]),
           0,array([1.,1.,0.]), 1,10)
print("%d " % len(y), end="")
print("%1.3e  %1.3e  %1.3e  " % tuple(y[0]),   end="")
print("%1.3e  %1.3e  %1.3e  " % tuple(y[-1]),   end="")

I should be getting the vector (0.3667, 1.001, 1.000) but instead I get (0.9955, 0.9955, -0.004536). I assume my code is breaking down inside the dirk function, particularly the input for Newton's Method.

Answer 1

Your big $F$ function and its derivative are slightly wrong. $F$ is only wrong in the formula, it should be $$ F(\delta_i) = \delta_i - f\left(t_0 + \alpha_i h, y_0 + h\beta_{ii}\delta_i + h\sum_{j=1}^{i-1} \beta_{ij}\kappa_j\right) $$ This is correct in the code. However, the derivative should be $$ F'(\delta_i) = I - h\beta_{ii}\partial_yf\left(t_0 + \alpha_i h, y_0 + h\beta_{ii}\delta_i + h\sum_{j=1}^{i-1} \beta_{ij}\kappa_j\right) $$ The factor $h\beta_{ii}$ is missing in the code. In both places.

Also, there is something strange happening with the K array. The insert operation is not doing what you think it is doing. Why not just declare K = np.zeros([n,m], dtype=float) and then assign K[0]=K0, K[i]=sol.

With these changes I get for y[-1] the values

3.67339302e-01  1.00054511e+00  1.00000000e+00  

with the 3rd order DIRK table from https://www.math.auckland.ac.nz/~butcher/ODE-book-2008/Tutorials/IRK.pdf

---

((Note that the format specifier %1.3e has 1 as the total length, thus is ignored. Typical formats leading to mostly straight columns are 8.5 or 12.8.))