# Butcher's Tableau to Diagonally Implicit Runge Kutta method

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.

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.))

• Thanks, stupid on my part for the derivative. As for the array I realised the values of K I was calling were wrong because K was 1D whereas I thought it was 2D. Commented Mar 2, 2021 at 15:18