I have the following system of equations obtained by implementing Sympletic Euler time scheme to wave equation. I want to model this in Fenics. Here 'u' is the displacement and 'p' is corresponding velocity.

$$ \dot{u} = p $$ $$ u(x, 0) = u_i$$ $$ p(x, 0) = 0$$ $$ u(x, t) = 0 \quad \text{on} \quad \partial \Omega$$

enter image description here

I read in a documentation that the velocity 'p' should be obtained first and thereafter 'u' should be computed. Also, that it must be ensured that we never have to explicitly compute and store the inverse of the mass matrix M.

Here, $M = uvdx$ and $K = \nabla u \nabla v dx$ .

Is this the correct way to implement the above system of equations? Can anyone kindly guide on how this can be corrected / refined.

u = interpolate(ui,V)           
V = FunctionSpace(mesh,"CG",2)

u_p = function(V)
u_new = function(V)
p = interpolate(Constant(0.0),V)

while t =< T
    bc.apply(M,u)        # apply boundary conditions    
    solve(M, u_p.vector(), u.vector)    
    p_new = p.vector() - dt*K*u_p.vector() 
    t += dt              # solve for next timestep
    u_new.vector()[:] = u.vector() + dt*p_new.vector()
    u.assign(u_new)      # update old vector
    p.vector()[:] = p_new
  • $\begingroup$ Just ran the code: Getting an error with the statement bc.apply(M,u) $\endgroup$
    – CRG
    Nov 7 '15 at 4:55