Second order time discretization problem

mathematics
2

You can convert to a first order system but then formally eliminate the extra velocity unknown to avoid solving in a mixed function space. The implicit midpoint rule gives you

\frac{u^{n+1} - u^n}{\Delta t} = \frac{1}{2}(\dot{u}^{n+1} + \dot{u}^{n})~~~\Rightarrow~~~\dot{u}^{n+1} = 2(\frac{u^{n+1} - u^n}{\Delta t}) - \dot{u}^n

which can be substituted into the left-hand side of

\frac{\dot{u}^{n+1} - \dot{u}^n}{\Delta t} = F((u^{n+1} + u^{n})/2)

to get a problem for u^{n+1} alone, given both u^n and \dot{u}^n from the previous time step. The following example demonstrates this for the scalar wave equation, and can be observed to converge with second-order accuracy in L^\infty(L^2) as h\sim\Delta t\to 0:

from dolfin import *
N = 64
T = 1.0
c = Constant(0.1)
Dt = Constant(T/N)
mesh = UnitIntervalMesh(N)
V = FunctionSpace(mesh,"CG",1)

# Set up time integration scheme:
u = Function(V)
u_old = Function(V)
u_mid = 0.5*(u + u_old)
udot_old = Function(V)
v = TestFunction(V)
udot = Constant(2/Dt)*u - Constant(2/Dt)*u_old - udot_old
uddot = (udot - udot_old)/Dt

# Weak problem residual for wave equation:
res = uddot*v*dx + (c**2)*inner(grad(u_mid),grad(v))*dx

# Define exact solution and set initial conditions:
t = Constant(0)
tv = variable(t)
x = SpatialCoordinate(mesh)[0]
x_shift = x-c*tv
u_ex = conditional(lt(x_shift,0.5),
                   conditional(gt(x_shift,0),sin(2*pi*(x_shift))**4,
                               Constant(0)),Constant(0))
udot_ex = diff(u_ex,tv)
u_old.assign(project(u_ex,V))
udot_old.assign(project(udot_ex,V))

# Time stepping loop:
lhsForm = derivative(res,u)
rhsForm = -res
du = Function(V)
for i in range(0,N):
    solve(lhsForm==rhsForm,du)
    u.assign(u+du)
    udot_old.assign(udot)
    u_old.assign(u)

# Check error:
t.assign(T)
e = u - u_ex
print(sqrt(assemble(e*e*dx)))

I’ll also remark that the implicit midpoint rule can be obtained as a special case of generalized-\alpha where the spectral radius for infinite time step (which all other parameters can be derived from) is set to 1.

2 Likes