Wave Equation

67. Wave Equation#

We consider the wave equation as typical example of second order hyperbolic equation. In contrast to the heat equation, we has a second order time derivative:

\[ \frac{\partial^2 u}{\partial t^2}(x,t) - \Delta u(x,t) = f(x,t) \qquad \forall \, x \in \Omega, \forall t \in (0,T) \]

Since it is second order in time, we need two initial conditions, for the function itself and for the time derivative:

\[\begin{eqnarray*} u(x,0) & = & u_0(x) \\ \frac{\partial u}{\partial t} (x,0) & = & v_0(x) \end{eqnarray*}\]

We proceed similar as for the heat equation

  • variational formulation in space

  • Galerkin approximation by the spatial basis

and obtain the second order ODE:

\[ M \ddot u(t) + A u(t) = f(t) \qquad \forall \, t \in (0,T) \]

with initial conditions

\[ u(0) = u_0 \qquad \text{and} \qquad \dot u(0) = v_0 \]

67.1. Newmark time-stepping method#

We reduce the second order ODE to a first order system. For this, we introduce the velocity function \(v = \dot u\):

\[\begin{eqnarray*} \dot u & = & v \\ M \dot v & = & f - A u \end{eqnarray*}\]

Now we apply the trapezoidal rule for time-discretization:

\[\begin{eqnarray*} \frac{u^j - u^{j-1}}{\tau} & = & \frac{v^j + v^{j-1}}{2} \\ M \frac{v^j - v^{j-1}}{\tau} & = & \frac{f^j+f^{j-1} - A (u^j+u^{j-1})}{2} \end{eqnarray*}\]

Expressing \(u^j\) from the first equation, and inserting it into the second we obtain

\[ M \frac{v^j - v^{j-1}}{\tau} = \frac{f^j+f^{j-1} - A (u^{j-1}+\tfrac{\tau}{2}(v^j+v^{j-1})+u^{j-1})}{2}, \]

and shifting the unknowns to the left we obtain a linear system for the update of velocity:

\[ (M + \tfrac{\tau^2}{4} A) (v^{j}-v^{j-1}) = \tau \left( \tfrac{1}{2} (f^j + f^{j-1}) - A u^{j-1} - \tfrac{\tau}{2} A v^{j-1} \right) \]

Although by introducing the system of ODEs, the size of the linear system to solve did not increase.

from ngsolve import *
from ngsolve.webgui import Draw
from netgen.occ import unit_square
from time import sleep

mesh = Mesh(unit_square.GenerateMesh(maxh=0.05))

tau = 0.01
tend = 1
u0 = exp(-100*( (x-0.5)**2 + (y-0.5)**2))
v0 = 0

fes = H1(mesh, order=3)
u,v = fes.TnT()
mform = u*v*dx
aform = grad(u)*grad(v)*dx

m = BilinearForm(mform).Assemble()
a = BilinearForm(aform).Assemble()
mstar = BilinearForm(mform+tau**2/4*aform).Assemble()
mstarinv = mstar.mat.Inverse()
f = LinearForm(fes).Assemble()

gfu = GridFunction(fes)
gfv = GridFunction(fes)
gfu.Set(u0)
gfv.Set(u0)

scene = Draw (gfu, deformation=True)
sleep(2)

for j in range(int(tend/tau)):
    gfu.vec.data += tau/2 * gfv.vec
    gfv.vec.data += tau * mstarinv * (f.vec - a.mat * gfu.vec)
    gfu.vec.data += tau/2 * gfv.vec

    scene.Redraw()
    sleep(0.1)
Contents