# Singly diagonally implicit Runge-Kutta methods for time stepping

Hello,

I am a bit confused about the implementation of SDIRK method in ngsolve when i studied the SDIRK method in i-tutorial unit 3.1

# introduction of SDIRK method

Let us consider more sophisticated time integration methods, SDIRK methods. To simplify presentation we set f=0.

SDIRK methods for the semi-discrete problem

\frac{d}{dt} u = M^{-1} F(u) = -M^{-1} \cdot A u

are of the form:

u^{n+1} = u^n + \Delta t M^{-1} \sum_{i=0}^{s-1} b_i k_i

with k_i = - A u_i , where u_i is the solution to

(M + a_{ii} \Delta t A) u_i = M u^n - \Delta t \sum_{i=0}^{i-1} a_{ij} k_j, \quad i=0,..,s-1.

The coefficients a,b and c are stored in the butcher tableau:

\begin{array}{c|cccc} c_0 & a_{00} & 0 & \ddots& \\ c_1 & a_{10} & a_{11} & 0 & \ddots \\ \vdots & \vdots & \vdots & \ddots & 0\\ c_{s-1} & a_{s-1,0} & a_{s-1,1} & \dots& a_{s-1,s-1} \\ \hline & b_{0} & b_1 & \dots& b_{s-1} \\ \end{array}

For an SDIRK method we have a^\ast = a_{ii},~i=0,..,s-1.

# implementation

• butcher tableau
class sdirk5: #order 4
stages = 5
a = [[1/4, 0, 0, 0, 0],
[1/2, 1/4, 0, 0, 0],
[17/50,-1/25, 1/4, 0, 0],
[371/1360, -137/2720, 15/544, 1/4,0],
[25/24, -49/48, 125/16, -85/12, 1/4]]
b = [25/24, -49/48, 125/16, -85/12, 1/4]
c = [1/4, 3/4, 11/20, 1/2, 1]
astar = 1/4

• time stepping function
def TimeStepping_app4(m, a, butchertab, dt,
initial_cond = None, t0 = 0, tend = 2,
nsamples = 10):
mstar.AsVector().data = m.mat.AsVector() + butchertab.astar * dt * a.mat.AsVector()
invmstar = mstar.Inverse(freedofs=fes.FreeDofs())
invmass = m.mat.Inverse(freedofs=fes.FreeDofs())

rhsi, Mu0, ui = gfu.vec.CreateVector(),gfu.vec.CreateVector(),gfu.vec.CreateVector()
k = [gfu.vec.CreateVector() for i in range(butchertab.stages)]

if initial_cond:
gfu.Set(initial_cond)
cnt = 0; time = t0
sample_int = int(floor(tend / dt / nsamples)+1)
gfuD = GridFunction(gfu.space)
gfut = GridFunction(gfu.space,multidim=0)
while time < tend - 0.5 * dt:
Mu0.data = m.mat * gfu.vec
for i in range(butchertab.stages):
rhsi.data = Mu0
for j in range(0,i):
rhsi.data += dt * butchertab.a[i][j] * k[j]
t.Set(time+butchertab.c[i]*dt)
gfu.Set(uD,BND)
ui.data = gfu.vec; rhsi.data -= mstar * ui
ui.data += invmstar * rhsi
k[i].data = - a.mat * ui
t.Set(time+dt)
gfu.Set(uD,BND)
Mu0.data -= m.mat * gfu.vec
for i in range(0,butchertab.stages):
Mu0.data += dt * butchertab.b[i] * k[i]
gfu.vec.data += invmass * Mu0
print("\r",time,end="")
if cnt % sample_int == 0:
cnt +=


# questions

• inconsistency between formulas and codes
when we use the results of previous sub-steps u_{j} to compute u_{i}, i think the formula should be of the following form (formula in the tutotial made a small mistake in \Sigma)
(M + a_{ii} \Delta t A) u_i = M u^n - \Delta t \sum_{j=0}^{i-1} a_{ij} k_j, \quad i=0,..,s-1.

The implementation of the formula is as follows (i add some comments with #====)

Mu0.data = m.mat * gfu.vec
for i in range(butchertab.stages):
#====now we have rhs = Mu^{n}====
rhsi.data = Mu0
for j in range(0,i):
#====why not rhsi.data -= dt * butchertab.a[i][j] * k[j]====
rhsi.data += dt * butchertab.a[i][j] * k[j]
t.Set(time+butchertab.c[i]*dt)
#====time-dependent boundary conditions====
gfu.Set(uD,BND)
ui.data = gfu.vec
rhsi.data -= mstar * ui
#====homogeneous + BC====
ui.data += invmstar * rhsi
k[i].data = - a.mat * ui
t.Set(time+dt)


Is there a sign problem? (rhsi.data += dt * butchertab.a[i][j] * k[j] for M u^n - \Delta t \sum_\limits{j=0}^{i-1} a_{ij} k_j, \quad i=0,..,s-1.)

Hi,

indeed the sign in the writen formula is wrong. The code should however be ok.
Thanks for the careful check an report!

Best,
Christoph

Thanks for your quick reply. In this tutorial, i think there is indeed only one small error with the subscript of sum in the formula((M + a_{ii} \Delta t A) u_i = M u^n - \Delta t \sum\limits_{i=0}^{i-1} a_{ij} k_j, \quad i=0,..,s-1. should be (M + a_{ii} \Delta t A) u_i = M u^n - \Delta t \sum\limits_{j=0}^{i-1} a_{ij} k_j, \quad i=0,..,s-1.) and a sign error in the corresponding code.