我通过将辛欧拉时间方案实施到波动方程获得了以下方程组。我想在 Fenics 中对此进行建模。这里'u'是位移，'p'是相应的速度。

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

我在文档中读到应该首先获得速度“p”，然后应该计算“u”。此外，必须确保我们永远不必显式计算和存储质量矩阵 M 的逆矩阵。

这里，和。$M = uvdx$$K = \nabla u \nabla v dx$

这是实现上述方程组的正确方法吗？任何人都可以指导如何纠正/完善这一点。

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