Second order time discretization problem

mathematics
1

I found the elastodynamic example in fenics tutorial:
\nabla \cdot \sigma + \rho b = \rho \ddot{u}
with
\sigma =\lambda \text{tr}(\varepsilon)\mathbb{1} + 2\mu\varepsilon and \varepsilon = (\nabla u + (\nabla u)^T)/2
I am going to use centered finite difference schemes instead of generalized -\alpha method as follows:
\ddot{u} = \frac{u^{i+2} - 2u^{i+1} + u^{i} }{\Delta^2{t}}
I assumed u^{i} is my initial condition which is displacement at time {i} but I am confused about how can I define u^{i+1} as my initial condition as well.

\quad

P.S. one of the other option comes to my mind is imposing change of variables, like:
w = \dot{u}
\dot{w} = \frac{w^{i+1}-w^{i}}{\Delta{t}}

And solve this with a mixed formulation.

\quad
Thank you in advance.

1 Like
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
3

Good morning,
Please I have the same problem I don’t know if you can help me I also have to discretize my equation I have chosen the implicit and explicit form of Crank-Nicolson however I have made change of variable to have the form of first order and apply the implicit and explicit form of Crank-Nicolson but I’m a little lost and I could find a document that has determined equation has the form a and l but I don’t know how we can declare the initial condition
please can you help me