Waves in blood vessels

The flow of blood in our bodies is basically fluid flow in a network of pipes. Unlike rigid pipes, the walls in the blood vessels are elastic and will increase their diameter when the pressure rises. The elastic forces will then push the wall back and accelerate the fluid. This interaction between the flow of blood and the deformation of the vessel wall results in waves traveling along our blood vessels.

A model for one-dimensional waves along blood vessels can be derived from averaging the fluid flow over the cross section of the blood vessels. Let \(x\) be a coordinate along the blood vessel and assume that all cross sections are circular, though with different radius \(R(x,t)\). The main quantities to compute is the cross section area \(A(x,t)\), the averaged pressure \(P(x,t)\), and the total volume flux \(Q(x,t)\). The area of this cross section is

\[A(x,t) = 2\pi\int_{0}^{R(x,t)} rdr,\]

Let \(v_x(x,t)\) be the velocity of blood averaged over the cross section at point \(x\). The volume flux, being the total volume of blood passing a cross section per time unit, becomes

\[Q(x,t) = A(x,t)v_x(x,t) \thinspace\]

Mass balance and Newton’s second law lead to the PDEs

(1)\[ \frac{\partial A}{\partial t} + \frac{\partial Q}{\partial x} = 0,\]

(2)\[ \frac{\partial Q}{\partial t} + \frac{\gamma+2}{\gamma + 1} \frac{\partial}{\partial x}\left(\frac{Q^2}{A}\right) + \frac{A}{\varrho}\frac{\partial P}{\partial x} = -2\pi (\gamma + 2)\frac{\mu}{\varrho}\frac{Q}{A},\]

where \(\gamma\) is a parameter related to the velocity profile, \(\varrho\) is the density of blood, and \(\mu\) is the dynamic viscosity of blood.

We have three unknowns \(A\), \(Q\), and \(P\), and two equations (1) and (2). A third equation is needed to relate the flow to the deformations of the wall. A common form for this equation is

(3)\[ \frac{\partial P}{\partial t} + \frac{1}{C} \frac{\partial Q}{\partial x} =0,\]

where \(C\) is the compliance of the wall, given by the constitutive relation

\[C = \frac{\partial A}{\partial P} + \frac{\partial A}{\partial t},\]

which require a relationship between \(A\) and \(P\). One common model is to view the vessel wall, locally, as a thin elastic tube subject to an internal pressure. This gives the relation

\[P=P_0 + \frac{\pi h E}{(1-\nu^2)A_0}(\sqrt{A} - \sqrt{A_0}),\]

where \(P_0\) and \(A_0\) are corresponding reference values when the wall is not deformed, \(h\) is the thickness of the wall, and \(E\) and \(\nu\) are Young’s modulus and Poisson’s ratio of the elastic material in the wall. The derivative becomes

\[C = \frac{\partial A}{\partial P} = \frac{2(1-\nu^2)A_0}{\pi h E}\sqrt{A_0} + 2\left(\frac{(1-\nu^2)A_0}{\pi h E}\right)^2(P-P_0) {\thinspace .}\]

Another (nonlinear) deformation model of the wall, which has a better fit with experiments, is

\[P = P_0\exp{(\beta (A/A_0 - 1))},\]

where \(\beta\) is some parameter to be estimated. This law leads to

\[C = \frac{\partial A}{\partial P} = \frac{A_0}{\beta P} {\thinspace .}\]

Reduction to standard wave equation. It is not uncommon to neglect the viscous term on the right-hand side of (2) and also the quadratic term with \(Q^2\) on the left-hand side. The reduced equations (2) and (3) form a first-order linear wave equation system:

\[C\frac{\partial P}{\partial t} = - \frac{\partial Q}{\partial x},\]

\[\frac{\partial Q}{\partial t} = - \frac{A}{\varrho}\frac{\partial P}{\partial x} {\thinspace .}\]

These can be combined into standard 1D wave equation PDE by differentiating the first equation with respect \(t\) and the second with respect to \(x\),

\[\frac{\partial}{\partial t}\left( CC\frac{\partial P}{\partial t} \right) = \frac{\partial}{\partial x}\left( \frac{A}{\varrho}\frac{\partial P}{\partial x}\right),\]

which can be approximated by

\[\frac{\partial^2 Q}{\partial t^2} = c^2\frac{\partial^2 Q}{\partial x^2},\quad c = \sqrt{\frac{A}{\varrho C}},\]

where the \(A\) and \(C\) in the expression for \(c\) are taken as constant reference values.

Electromagnetic waves

Light and radio waves are governed by standard wave equations arising from Maxwell’s general equations. When there are no charges and no currents, as in a vacuum, Maxwell’s equations take the form

\[\begin{split}\nabla\cdot\pmb{E} &= 0,\\ \nabla\cdot\pmb{B} &= 0,\\ \nabla\times\pmb{E} &= -\frac{\partial\pmb{B}}{\partial t},\\ \nabla\times\pmb{B} &= \mu_0\epsilon_0\frac{\partial\pmb{E}}{\partial t},\end{split}\]

where \(\epsilon_0=8.854187817620\cdot 10^{-12}\) (F/m) is the permittivity of free space, also known as the electric constant, and \(\mu_0=1.2566370614\cdot 10^{-6}\) (H/m) is the permeability of free space, also known as the magnetic constant. Taking the curl of the two last equations and using the identity

\[\nabla\times (\nabla\times \pmb{E}) = \nabla(\nabla \cdot \pmb{E}) - \nabla^2\pmb{E} = - \nabla^2\pmb{E}\hbox{ when }\nabla\cdot\pmb{E}=0,\]

immediately gives the wave equation governing the electric and magnetic field:

\[\frac{\partial^2\pmb{E}}{\partial t^2} = c^2\frac{\partial^2\pmb{E}}{\partial x^2},\]

\[\frac{\partial^2\pmb{E}}{\partial t^2} = c^2\frac{\partial^2\pmb{E}}{\partial x^2},\]

with \(c=1/\sqrt{\mu_0\epsilon_0}\) as the velocity of light. Each component of \(\pmb{E}\) and \(\pmb{B}\) fulfills a wave equation and can hence be solved independently.

Exercises (4)

Exercise 14: Simulate waves on a non-homogeneous string

Simulate waves on a string that consists of two materials with different density. The tension in the string is constant, but the density has a jump at the middle of the string. Experiment with different sizes of the jump and produce animations that visualize the effect of the jump on the wave motion.

Hint. According to the section Waves on a string, the density enters the mathematical model as \(\varrho\) in \(\varrho u_{tt} = Tu_{xx}\), where \(T\) is the string tension. Modify, e.g., the wave1D_u0_sv.py code to incorporate the tension and two density values. Make a mesh function rho with density values at each spatial mesh point. A value for the tension may be 150 N. Corresponding density values can be computed from the wave velocity estimations in the guitar function in the wave1D_u0_sv.py file.

Filename: wave1D_u0_sv_discont.py.

Exercise 15: Simulate damped waves on a string

Formulate a mathematical model for damped waves on a string. Use data from the section Running a case, and tune the damping parameter so that the string is very close to the rest state after 15 s. Make a movie of the wave motion. Filename: wave1D_u0_sv_damping.py.

Exercise 16: Simulate elastic waves in a rod

A hammer hits the end of an elastic rod. The exercise is to simulate the resulting wave motion using the model (8.8) from the section Elastic waves in a rod. Let the rod have length \(L\) and let the boundary \(x=L\) be stress free so that \(\sigma_{xx}=0\), implying that \(\partial u/\partial x=0\). The left end \(x=0\) is subject to a strong stress pulse (the hammer), modeled as

\[\begin{split}\sigma_{xx}(t) = \left\lbrace\begin{array}{ll} S,& 0 < t \leq t_s,\\ 0, & t > t_s \end{array}\right.\end{split}\]

The corresponding condition on \(u\) becomes \(u_x= S/E\) for \(t\leq t_s\) and zero afterwards (recall that \(\sigma_{xx} = Eu_x\)). This is a non-homogeneous Neumann condition, and you will need to approximate this condition and combine it with the scheme (the ideas and manipulations follow closely the handling of a non-zero initial condition \(u_t=V\) in wave PDEs or the corresponding second-order ODEs for vibrations). Filename: wave_rod.py.

Exercise 17: Simulate spherical waves

Implement a model for spherically symmetric waves using the method described in the section Spherical waves. The boundary condition at \(r=0\) must be \(\partial u/\partial r=0\), while the condition at \(r=R\) can either be \(u=0\) or a radiation condition as described in Problem 20: Implement open boundary conditions. The \(u=0\) condition is sufficient if \(R\) is so large that the amplitude of the spherical wave has become insignificant. Make movie(s) of the case where the source term is located around \(r=0\) and sends out pulses

\[\begin{split}f(r,t) = \left\lbrace\begin{array}{ll} Q\exp{(-\frac{r^2}{2\Delta r^2})}\sin\omega t,& \sin\omega t\geq 0\\ 0, & \sin\omega t < 0 \end{array}\right.\end{split}\]

Here, \(Q\) and \(\omega\) are constants to be chosen.

Hint. Use the program wave1D_u0_sv.py as a starting point. Let solver compute the \(v\) function and then set \(u=v/r\). However, \(u=v/r\) for \(r=0\) requires special treatment. One possibility is to compute u[1:] = v[1:]/r[1:] and then set u[0]=u[1]. The latter makes it evident that \(\partial u/\partial r = 0\) in a plot.

Filename: wave1D_spherical.py.

Exercise 18: Explain why numerical noise occurs

The experiments performed in Exercise 8: Send pulse waves through a layered medium shows considerable numerical noise in the form of non-physical waves, especially for \(s_f=4\) and the plug pulse or the half a “cosinehat” pulse. The noise is much less visible for a Gaussian pulse. Run the case with the plug and half a “cosinehat” pulses for \(s_f=1\), \(C=0.9, 0.25\), and \(N_x=40,80,160\). Use the numerical dispersion relation to explain the observations. Filename: pulse1D_analysis.pdf.

Exercise 19: Investigate harmonic averaging in a 1D model

Harmonic means are often used if the wave velocity is non-smooth or discontinuous. Will harmonic averaging of the wave velocity give less numerical noise for the case \(s_f=4\) in Exercise 8: Send pulse waves through a layered medium? Filenames: pulse1D_harmonic.pdf, pulse1D_harmonic.py.

Problem 20: Implement open boundary conditions

To enable a wave to leave the computational domain and travel undisturbed through the boundary \(x=L\), one can in a one-dimensional problem impose the following condition, called a radiation condition or open boundary condition:

(4)\[ \frac{\partial u}{\partial t} + c\frac{\partial u}{\partial x} = 0{\thinspace .}\]

The parameter \(c\) is the wave velocity.

Show that (4) accepts a solution \(u = g_R(x-ct)\), but not \(u = g_L(x+ct)\). This means that (4) will allow any right-going wave \(g_R(x-ct)\) to pass through the boundary.

A corresponding open boundary condition for a left-going wave through \(x=0\) is

(5)\[ \frac{\partial u}{\partial t} + c\frac{\partial u}{\partial x} = 0{\thinspace .}\]

The condition (4) can be discretized by centered differences at the spatial end point \(i=N_x\), corresponding to \(x=x_R\):

(6)\[ [D_{2t}u + cD_{2x}u =0]^n_{N_x}\]\[ {\thinspace .}\]

Eliminate the fictitious value \(u_{N_x+1}^n\) by using the discrete equation at the same point. The equation for the first step, \(u_i^1\), is in principal affected, but we can then use the condition \(u_{N_x}=0\) since the wave has not yet reached the right boundary.

A corresponding open boundary condition for a left-going wave through \(x=0\) is

(7)\[ \frac{\partial u}{\partial t} - c\frac{\partial u}{\partial x} = 0{\thinspace .}\]

Implement a solver that incorporates the conditions (6) and (7). Start with some peak-shaped Gaussian function in the middle of the domain as \(I(x)\) and demonstrate that waves travel undisturbed out of the domain at \(x=L\) and \(x=0\). Make a nose test for checking that the surface is flat after a certain time.

Remark. The condition (4) works perfectly in 1D when \(c\) is known. In 2D and 3D, however, the condition reads \(u_t + c_x u_x + c_y u_y=0\), where \(c_x\) and \(c_y\) are the wave speeds in the \(x\) and \(y\) directions. Estimating these components (i.e., the direction of the wave) is often challenging. Other methods are normally used in 2D and 3D to let waves move out of a computational domain. Filename: wave1D_open_BC.py.

Problem 21: Earthquake-generated tsunami over a subsea hill

A subsea earthquake leads to an immediate lift of the water surface, see Figure Sketch of initial water surface due to a subsea earthquake. The lifted water surface splits into two tsunamis, one traveling to the right and one to the left, as depicted in Figure An initial surface elevation is split into two waves. Since tsunamis are normally very long waves, compared to the depth, with a small amplitude, compared to the wave length, the wave equation model described in the section The linear shallow water equations is relevant:

\[\eta_{tt} = (gH(x)\eta_x)_x,\]

where \(g\) is the acceleration of gravity, and \(H(x)\) is the still water depth.

Sketch of initial water surface due to a subsea earthquake

An initial surface elevation is split into two waves

To simulate the right-going tsunami, we can impose a symmetry boundary at \(x=0\): \(\partial\eta\ \partial x =0\). We then simulate the wave motion in \([0,L]\). Unless the ocean ends at \(x=L\), the waves should travel undisturbed through the boundary \(x=L\). A radiation condition as explained in Problem 20: Implement open boundary conditions can be used for this purpose. Alternatively, one can just stop the simulations before the wave hits the boundary at \(x=L\). In that case it does not matter what kind of boundary condition we use at \(x=L\). Imposing \(\eta =0\) and stopping the simulations when \(|\eta_i^n| > \epsilon\), \(i=N_x-1\), is a possibility (\(\epsilon\) is a small parameter).

The shape of the initial surface can be taken as a Gaussian function,

\[I(x;I_0,I_a,I_m,I_s) = I_0 + I_a\exp{\left(-\left(\frac{x-I_m}{I_s}\right)^2\right)},\]

with \(I_m=0\) reflecting the location of the peak of \(I(x)\) and \(I_s\) being a measure of the width of the function \(I(x)\) (\(I_s\) is \(\sqrt{2}\) times the standard deviation of the familiar normal distribution curve).

Now we extend the problem with a hill at the sea bottom, see Figure Sketch of an earthquake-generated tsunami passing over a subsea hill. The wave speed \(c=\sqrt{gH(x)} = \sqrt{g(H_0-B(x))}\) will then be reduced in the shallow water above the hill.

Sketch of an earthquake-generated tsunami passing over a subsea hill

One possible form of the hill is a Gaussian function,

(8)\[ B(x;B_0,B_a,B_m,B_s) = B_0 + B_a\exp{\left(-\left(\frac{x-B_m}{B_s}\right)^2\right)},\]

but many other shapes are also possible, e.g., a “cosine hat” where

(9)\[ B(x; B_0, B_a, B_m, B_s) = B_0 + B_a\cos{\left( \pi\frac{x-B_m}{2B_s}\right)},\]

when \(x\in [B_m - B_s, B_m + B_s]\) while \(B=B_0\) outside this interval.

Also an abrupt construction may be tried:

(10)\[ B(x; B_0, B_a, B_m, B_s) = B_0 + B_a,\]

for \(x\in [B_m - B_s, B_m + B_s]\) while \(B=B_0\) outside this interval.

The wave1D_dn_vc.py program can be used as starting point for the implementation. Visualize both the bottom topography and the water surface elevation in the same plot. Allow for a flexible choice of bottom shape: (8), (9), (10), or \(B(x)=B_0\) (flat).

The purpose of this problem is to explore the quality of the numerical solution \(\eta^n_i\) for different shapes of the bottom obstruction. The “cosine hat” and the box-shaped hills have abrupt changes in the derivative of \(H(x)\) and are more likely to generate numerical noise than the smooth Gaussian shape of the hill. Investigate if this is true. Filenames: tsunami1D_hill.py, tsunami1D_hill.pdf.

Problem 22: Earthquake-generated tsunami over a 3D hill

This problem extends Problem 21: Earthquake-generated tsunami over a subsea hill to a three-dimensional wave phenomenon, governed by the 2D PDE (9.5). We assume that the earthquake arise from a fault along the line \(x=0\) in the \(xy\)-plane so that the initial lift of the surface can be taken as \(I(x)\) in Problem 21: Earthquake-generated tsunami over a subsea hill. That is, a plane wave is propagating to the right, but will experience bending because of the bottom.

The bottom shape is now a function of \(x\) and \(y\). An “elliptic” Gaussian function in two dimensions, with its peak at \((B_{mx}, B_{my})\), generalizes (8):

(11)\[ B(x;B_0,B_a,B_{mx}, B_{my} ,B_s, b) = B_0 + B_a\exp{\left(-\left(\frac{x-B_{mx}}{B_s}\right)^2 -\left(\frac{y-B_{my}}{bB_s}\right)^2\right)},\]

where \(b\) is a scaling parameter: \(b=1\) gives a circular Gaussian function with circular contour lines, while \(b\neq 1\) gives an elliptic shape with elliptic contour lines.

The “cosine hat” (9) can also be generalized to

(12)\[ B(x; B_0, B_a, B_{mx}, B_{my}, B_s) = B_0 + B_a\cos{\left( \pi\frac{x-B_{mx}}{2B_s}\right)} \cos{\left( \pi\frac{y-B_{my}}{2B_s}\right)},\]

when \(0 \leq \sqrt{x^2+y^2} \leq B_s\) and \(B=B_0\) outside this circle.

A box-shaped obstacle means that

(13)\[ B(x; B_0, B_a, B_m, B_s, b) = B_0 + B_a\]

for \(x\) and \(y\) inside a rectangle

\[B_{mx}-B_s \leq x \leq B_{mx} + B_s,\quad B_{my}-bB_s \leq y \leq B_{my} + bB_s,\]

and \(B=B_0\) outside this rectangle. The \(b\) parameter controls the rectangular shape of the cross section of the box.

Note that the initial condition and the listed bottom shapes are symmetric around the line \(y=B_{my}\). We therefore expect the surface elevation also to be symmetric with respect to this line. This means that we can halve the computational domain by working with \([0,L_x]\times [0, B_{my}]\). Along the upper boundary, \(y=B_{my}\), we must impose the symmetry condition \(\partial \eta/\partial n=0\). Such a symmetry condition (\(-\eta_x=0\)) is also needed at the \(x=0\) boundary because the initial condition has a symmetry here. At the lower boundary \(y=0\) we also set a Neumann condition (which becomes \(-\eta_y=0\)). The wave motion is to be simulated until the wave hits the reflecting boundaries where \(\partial\eta/\partial n =\eta_x =0\) (one can also set \(\eta =0\) - the particular condition does not matter as long as the simulation is stopped before the wave is influenced by the boundary condition).

Visualize the surface elevation. Investigate how different hill shapes, different sizes of the water gap above the hill, and different resolutions \(\Delta x = \Delta y = h\) and \(\Delta t\) influence the numerical quality of the solution. Filenames: tsunami2D_hill.py, tsunami2D_hill.pdf.

Problem 23: Investigate Matplotlib for visualization

Play with native Matplotlib code for visualizing 2D solutions of the wave equation with variable wave velocity. See if there are effective ways to visualize both the solution and the wave velocity. Filename: tsunami2D_hill_mpl.py.

Problem 24: Investigate visualization packages

Create some fancy 3D visualization of the water waves and the subsea hill in Problem 22: Earthquake-generated tsunami over a 3D hill. Try to make the hill transparent. Possible visualization tools are

• Mayavi
• Paraview
• OpenDX

Filename: tsunami2D_hill_viz.py.

Problem 25: Implement loops in compiled languages

Extend the program from Problem 22: Earthquake-generated tsunami over a 3D hill such that the loops over mesh points, inside the time loop, are implemented in compiled languages. Consider implementations in Cython, Fortran via f2py, C via Cython, C via f2py, C/C++ via Instant, and C/C++ via scipy.weave. Perform efficiency experiments to investigate the relative performance of the various implementations. It is often advantageous to normalize CPU times by the fastest method on a given mesh. Filename: tsunami2D_hill_compiled.py.

Exercise 26: Simulate seismic waves in 2D

The goal of this exercise is to simulate seismic waves using the PDE model (8.16) in a 2D \(xz\) domain with geological layers. Introduce \(m\) horizontal layers of thickness \(h_i\), \(i=0,\ldots,m-1\). Inside layer number \(i\) we have a vertical wave velocity \(c_{z,i}\) and a horizontal wave velocity \(c_{h,i}\). Make a program for simulating such 2D waves. Test it on a case with 3 layers where

\[c_{z,0}=c_{z,1}=c_{z,2},\quad c_{h,0}=c_{h,2},\quad c_{h,1} \ll c_{h,0} {\thinspace .}\]

Let \(s\) be a localized point source at the middle of the Earth’s surface (the upper boundary) and investigate how the resulting wave travels through the medium. The source can be a localized Gaussian peak that oscillates in time for some time interval. Place the boundaries far enough from the expanding wave so that the boundary conditions do not disturb the wave. Then the type of boundary condition does not matter, except that we physically need to have \(p=p_0\), where \(p_0\) is the atmospheric pressure, at the upper boundary. Filename: seismic2D.py.

Project 27: Model 3D acoustic waves in a room

The equation for sound waves in air is derived in the section Sound waves in liquids and gases and reads

\[p_{tt} = c^2\nabla^2 p,\]

where \(p(x,y,z,t)\) is the pressure and \(c\) is the speed of sound, taken as 340 m/s.

However, sound is absorbed in the air due to relaxation of molecules in the gas. A model for simple relaxation, valid for gases consisting only of one type of molecules, is a term \(c^2\tau_s\nabla^2 p_t\) in the PDE, where \(\tau_s\) is the relaxation time. If we generate sound from, e.g., a loudspeaker in the room, this sound source must also be added to the governing equation.

The PDE with the mentioned type of damping and source then becomes

\[p_tt = c^2\nabla^p + c^2\tau_s\nabla^2 p_t + f,\]

where \(f(x,y,z,t)\) is the source term.

The walls can absorb some sound. A possible model is to have a “wall layer” (thicker than the physical wall) outside the room where \(c\) is changed such that some of the wave energy is reflected and some is absorbed in the wall. The absorption of energy can be taken care of by adding a damping term \(bp_t\) in the equation:

\[p_tt + bp_t = c^2\nabla^p + c^2\tau_s\nabla^2 p_t + f{\thinspace .}\]

Typically, \(b=0\) in the room and \(b>0\) in the wall. A discontinuity in \(b\) or \(c\) will give rise to reflections. It can be wise to use a constant \(c\) in the wall to control reflections because of the discontinuity between \(c\) in the air and in the wall, while \(b\) is gradually increased as we go into the wall to avoid reflections because of rapid changes in \(b\). At the outer boundary of the wall the condition \(p=0\) or \(\partial p/\partial n=0\) can be imposed. The waves should anyway be approximately dampened to \(p=0\) this far out in the wall layer.

There are two strategies for discretizing the \(\nabla^2 p_t\) term: using a center difference between times \(n+1\) and \(n-1\) (if the equation is sampled at level \(n\)), or use a one-sided difference based on levels \(n\) and \(n-1\). The latter has the advantage of not leading to any equation system, while the former is second-order accurate as the scheme for the simple wave equation \(p_tt = c^2\nabla^2 p\). To avoid an equation system, go for the one-sided difference such that the overall scheme becomes explicit and only of first order in time.

Develop a 3D solver for the specified PDE and introduce a wall layer. Test the solver with the method of manufactured solutions. Make some demonstrations where the wall reflects and absorbs the waves (reflection because of discontinuity in \(b\) and absorption because of growing \(b\)). Experiment with the impact of the \(\tau_s\) parameter. Filename: acoustics.py.

Project 28: Solve a 1D transport equation

We shall study the wave equation

(14)\[ u_t + cu_x = 0,\quad x\in (0,L],\ t\in (0, T],\]

with initial condition

\[u(x,0) = I(x),\quad x\in [0,L],\]

and one periodic boundary condition

\[u(0,t) = u(L,t) {\thinspace .}\]

This boundary condition means that what goes out of the domain at \(x=L\) comes in at \(x=0\). Roughly speaking, we need only one boundary condition because of the spatial derivative is of first order only.

Physical interpretation. The parameter \(c\) can be constant or variable, \(c=c(x)\). The equation (14) arises in transport problems where a quantity \(u\), which could be temperature or concentration of some contaminant, is transported with the velocity \(c\) of a fluid. In addition to the transport imposed by “travelling with the fluid”, \(u\) may also be transported by diffusion (such as heat conduction or Fickian diffusion), but we have in the model \(u_t + cu_x\) assumed that diffusion effects are negligible, which they often are.

A widely used numerical scheme for (14) applies a forward difference in time and a backward difference in space when \(c>0\):

(15)\[ [D_t^+ u + cD_x^-u = 0]_i^n\]\[ {\thinspace .}\]

For \(c<0\) we use a forward difference in space: \([cD_x^+u]_i^n\).

We shall hereafter assume that \(=c(x)>0\).

To compute (20) we need to integrate \(1/c\) to obtain \(C\) and then compute the inverse of \(C\).

The inverse function computation can be easily done if we first think discretely. Say we have some function \(y=g(x)\) and seeks its inverse. Plotting \((x_i,y_i)\), where \(y_i=g(x_i)\) for some mesh points \(x_i\), displays \(g\) as a function of \(x\). The inverse function is simply \(x\) as a function of \(g\), i.e., the curve with points \((y_i,x_i)\). We can therefore quickly compute points at the curve of the inverse function. One way of extending these points to a continuous function is to assume a linear variation (known as linear interpolation) between the points (which actually means to draw straight lines between the points, exactly as done by a plotting program).

The function wrap2callable in scitools.std can take a set of points and return a continuous function that corresponds to linear variation between the points. The computation of the inverse of a function \(g\) on \([0,L]\) can then be done by

def inverse(g, domain, resolution=101):
    x = linspace(domain[0], domain[L], resolution)
    y = g(x)
    from scitools.std import wrap2callable
    g_inverse = wrap2callable((y, x))
    return g_inverse

To compute \(C(x)\) we need to integrate \(1/c\), which can be done by a Trapezoidal rule. Suppose we have computed \(C(x_i)\) and need to compute \(C(x_{i+1})\). Using the Trapezoidal rule with \(m\) subintervals over the integration domain \([x_i,x_{i+1}]\) gives

(16)\[ C(x_{i+1}) = C(x_i) + \int_{x_i}^{x_{i+1}} \frac{dx}{c} \approx h\left( \frac{1}{2}\frac{1}{c(x_i)} + \frac{1}{2}\frac{1}{c(x_{i+1})} + \sum_{j=1}^{m-1} \frac{1}{c(x_i + jh)}\right),\]

where \(h=(x_{i+1}-x_i)/m\) is the length of the subintervals used for the integral over \([x_i,x_{i+1}]\). We observe that (16) is a difference equation which we can solve by repeatedly applying (16) for \(i=0,1,\ldots,N_x-1\) if a mesh \(x_0,x_,\ldots,x_{N_x}\) is prescribed. Note that \(C(0)=0\).

a) Show that under the assumption of \(a=\hbox{const}\),

(17)\[ u(x,t) = I(x - ct)\]

fulfills the PDE as well as the initial and boundary condition (provided \(I(0)=I(L)\)).

b) Set up a computational algorithm and implement it in a function. Assume \(a\) is constant and positive.

c) Test implementation by using the remarkable property that the numerical solution is exact at the mesh points if \(\Delta t = c^{-1}\Delta x\).

d) Make a movie comparing the numerical and exact solution for the following two choices of initial conditions:

(18)\[ I(x) = \left\lbrack\sin\left(\pi\frac{x}{L}\right)\right\rbrack^{2n}\]

where \(n\) is an integer, typically \(n=5\), and

(19)\[ I(x) = \exp{\left( -\frac{(x-L/2)^2}{2\sigma2}\right)} {\thinspace .}\]

Choose \(\Delta t = c^{-1}\Delta x, 0.9c^{-1}\Delta x, 0.5c^{-1}\Delta x\).

e) The performance of the suggested numerical scheme can be investigated by analyzing the numerical dispersion relation. Analytically, we have that the Fourier component

\[u(x,t) = e^{i(kx-\omega t)},\]

is a solution of the PDE if \(\omega = kc\). This is the analytical dispersion relation. A complete solution of the PDE can be built by adding up such Fourier components with different amplitudes, where the initial condition \(I\) determines the amplitudes. The solution \(u\) is then represented by a Fourier series.

A similar discrete Fourier component at \((x_p,t_n)\) is

\[u_p^q = e^{i(kp\Delta x -\tilde\omega n\Delta t)},\]

where in general \(\tilde\omega\) is a function of \(k\), \(\Delta t\), and \(\Delta x\), and differs from the exact \(\omega =kc\).

Insert the discrete Fourier component in the numerical scheme and derive an expression for \(\tilde\omega\), i.e., the discrete dispersion relation. Show in particular that if the \(\Delta t/(c\Delta x)=1\), the discrete solution coincides with the exact solution at the mesh points, regardless of the mesh resolution (!). Show that if the stability condition

\[\frac{\Delta t}{c\Delta x}\leq 1,\]

the discrete Fourier component cannot grow (i.e., \(\tilde\omega\) is real).

f) Write a test for your implementation where you try to use information from the numerical dispersion relation.

g) Set up a computational algorithm for the variable coefficient case and implement it in a function. Make a test that the function works for constant \(a\).

h) It can be shown that for an observer moving with velocity \(c(x)\), \(u\) is constant. This can be used to derive an exact solution when \(a\) varies with \(x\). Show first that

(20)\[ u(x,t) = f(C(x) - t),\]

where

\[C'(x) = \frac{1}{c(x)},\]

is a solution of (14) for any differentiable function \(f\).

Solution. Let \(\xi = C(x) - t\). We have that

\[u_t = f'(\xi)(-1),\]

while

\[u_x = f'(\xi)C'(x) = f'(\xi)\frac{1}{c(x)},\]

implying that \(au_x = f'(\xi)\). Then we have \(u_t + cu_x= -f'(\xi) + f'(\xi) = 0\).

i) Use the initial condition to show that an exact solution is

\[u(x,t) = I(C^{-1}(C(x)-t)),\]

with \(C^{-1}\) being the inverse function of \(C = \int c^{1}dx\). Since \(C(x)\) is an integral \(\int_0^x (1/c)dx\), \(C(x)\) is monotonically increasing and there exists hence an inverse function \(C^{-1}\) with values in \([0,L]\).

Solution. In general we have \(u(x,t) = f(C(x)-t)\) and the solution is of this form with \(f(\xi)=I(C^{-1}(\xi))\). Moreover, at \(t=0\) we have \(I(C^{-1}(C(x)))=I(x)\), which is the required initial condition.

j) Implement a function for computing \(C(x_i)\) and one for computing \(C^{-1}(x)\) for any \(x\). Use these two functions for computing the exact solution \(I(C^{-1}(C(x)-t))\). End up with a function u_exact_variable_c(x, n, c, I) that returns the value of \(I(C^{-1}(C(x)-t_n))\).

k) Make movies showing a comparison of the numerical and exact solutions for the two initial conditions (18) and (19). Choose \(\Delta t = \Delta x /\max_{0,L} c(x)\) and the velocity of the medium as

1. \(c(x) = 1 + \epsilon\sin(k\pi x/L)\), \(\epsilon <1\),
2. \(c(x) = 1 + I(x)\), where \(I\) is given by (18) or (19).

The PDE \(u_t + cu_x=0\) expresses that the initial condition \(I(x)\) is transported with velocity \(c(x)\).

Filename: advec1D.py.

Problem 29: General analytical solution of a 1D damped wave equation

We consider an initial-boundary value problem for the damped wave equation:

\[\begin{split}u_{tt} +bu_t &= c^2 u_{xx}, \quad &x\in (0,L),\ t\in (0,T]\\ u(0,t) &= 0,\\ u(L,t) &=0,\\ u(x, 0) &= I(x),\\ u_t(x, 0) &= V(x){\thinspace .}\end{split}\]

Here, \(b\geq 0\) and \(c\) are given constants. The aim is to derive a general analytical solution of this problem. Familiarity with the method of separation of variables for solving PDEs will be assumed.

a) Seek a solution on the form \(u(x,t)=X(x)T(t)\). Insert this solution in the PDE and show that it leads to two differential equations for \(X\) and \(T\):

\[T'' + bT' + \lambda T = 0,\quad c^2 X'' +\lambda X = 0,\]

with \(X(0)=X(L)=0\) as boundary conditions, and \(\lambda\) as a constant to be determined.

b) Show that \(X(x)\) is on the form

\[X_n(x) = C_n\sin kx,\quad k = \frac{n\pi}{L},\quad n=1,2,\ldots\]

where \(C_n\) is an arbitrary constant.

c) Under the assumption that \((b/2)^2 < k^2\), show that \(T(t)\) is on the form

\[T_n(t) = e^{-{\frac{1}{2}}bt}(a_n\cos\omega t + b_n\sin\omega t), \quad\omega = \sqrt{k^2 - \frac{1}{4}b^2},\quad n=1,2,\ldots\]

The complete solution is then

\[u(x,t) = \sum_{n=1}^\infty \sin kx e^{-{\frac{1}{2}}bt}( A_n\cos\omega t + B_n\sin\omega t),\]

where the constants \(A_n\) and \(B_n\) must be computed from the initial conditions.

d) Derive a formula for \(A_n\) from \(u(x,0)=I(x)\) and developing \(I(x)\) as a sine Fourier series on \([0,L]\).

e) Derive a formula for \(B_n\) from \(u_t(x,0)=V(x)\) and developing \(V(x)\) as a sine Fourier series on \([0,L]\).

f) Calculate \(A_n\) and \(B_n\) from vibrations of a string where \(V(x)=0\) and

(21)\[\begin{split} I(x) = \left\lbrace \begin{array}{ll} ax/x_0, & x < x_0,\\ a(L-x)/(L-x_0), & \hbox{otherwise} \end{array}\right.\end{split}\]

g) Implement the series for \(u(x,t)\) in a function u_series(x, t, tol=1E-10), where tol is a tolerance for truncating the series. Simply sum the terms until \(|a_n|\) and \(|b_b|\) both are less than tol.

h) What will change in the derivation of the analytical solution if we have \(u_x(0,t)=u_x(L,t)=0\) as boundary conditions? And how will you solve the problem with \(u(0,t)=0\) and \(u_x(L,t)=0\)?

Filename: damped_wave1D.pdf.

Problem 30: General analytical solution of a 2D damped wave equation

Carry out Problem 29: General analytical solution of a 1D damped wave equation in the 2D case: \(u_{tt}+bu_t = c^2(u_{xx}+u_{yy})\), where \((x,y)\in (0,L_x)\times (0,L_y)\). Assume a solution on the form \(u(x,y,t)=X(x)Y(y)T(t)\).

Filename: damped_wave2D.pdf.