Variable Depth Shallow Water Wave Equation
Introduction
We consider here the problem of waves reflected by a region of variable depth in an otherwise uniform depth region assuming the equations of Shallow Depth.
Equations
We begin with the shallow depth equation
subject to the initial conditions
where [math]\displaystyle{ \zeta }[/math] is the displacement, [math]\displaystyle{ \rho }[/math] is the string density and [math]\displaystyle{ h(x) }[/math] is the variable depth (note that we are unifying the variable density string and the wave equation in variable depth because the mathematical treatment is identical).
Waves in a finite basin
We consider the problem of waves in a finite basin [math]\displaystyle{ 0\lt x\lt 1 }[/math]. At the edge of the basin the boundary conditions are
[math]\displaystyle{ \left.\partial_x \zeta\right|_{x=0} = \left.\partial_x \zeta\right|_{x=1} =0 }[/math]
.
We solve the equations by expanding in the modes for the basin which satisfy
normalised so that
The solution is then given by
where we have assumed that [math]\displaystyle{ \lambda_0 = 0 }[/math].
Calculation of [math]\displaystyle{ \zeta_n }[/math]
We can calculate the eigenfunctions [math]\displaystyle{ \zeta_n }[/math] by an expansion in the modes for the case of uniform depth. We use the Rayleigh-Ritz method. The eigenfunctions are local minimums of
[math]\displaystyle{ J[\zeta] = \int_0^1 \frac{1}{2}\left\{ \left(h(x) \partial_x \zeta\right)^2 - \lambda \zeta^2 \right\} }[/math]
subject to the boundary conditions that the normal derivative vanishes (where [math]\displaystyle{ \lambda }[/math] is the eigenvalue).
We expand the displacement in the eigenfunctions for constant depth [math]\displaystyle{ h=1 }[/math]
[math]\displaystyle{ \zeta = \sum_{n=1}^{N} a_n \psi_n(x) }[/math]
where
[math]\displaystyle{ \psi_n = \sqrt{2} \cos( (n-1) \pi x),\,\,n\ne 1 }[/math]
[math]\displaystyle{ \psi_0 = 1,\, }[/math]
and substitute this expansion into the variational equation we obtain
[math]\displaystyle{ \mathbf{M} \vec{a} = \lambda \vec{a} }[/math]
where the elements of the matrix M are
[math]\displaystyle{ M_{mn} = \int_0^1 \left\{ \left(\partial_x \psi_m h(x) \partial_x \psi_n\right) \right\} }[/math]
Matlab code
Waves in an infinite basin
We assume that the depth is constant and equal to one outside the region [math]\displaystyle{ 0\lt x\lt 1 }[/math]. We can therefore write the wave as
Solution using Separation of Variables
Taking a separable solution gives the eigenvalue problem
[math]\displaystyle{ \partial_x \left( h(x) \partial_x\zeta \right) = -\kappa^{2}\zeta \quad (1) }[/math]
Given boundary conditions [math]\displaystyle{ \zeta (0) = a }[/math] and [math]\displaystyle{ \zeta (1) = b }[/math] we can take [math]\displaystyle{ \zeta = (b-a)x + a + u }[/math] With [math]\displaystyle{ u }[/math] satisfying [math]\displaystyle{ u (0) = u (1) = 0 }[/math]
Substituting this form into (1) gives
[math]\displaystyle{ (b-a)\partial_xh(x)+\partial_x(h(x)\partial_xu) = -\kappa^{2}\left((b-a)x+a+u\right) }[/math]
Or, on rearranging
[math]\displaystyle{ \partial_x(h(x)\partial_xu)+\kappa^{2}u = -(b-a)\partial_xh(x)-\kappa^{2}\left((b-a)x+a\right) \quad (2) }[/math]
Now consider the homogenous Sturm-Liouville problem for u
[math]\displaystyle{ \partial_x(h(x)\partial_xu)+\lambda u = 0\quad u(0)=u(1)=0 \quad (3) }[/math]
By Sturm-Liouville theory this has an infinite set of eigenvalues [math]\displaystyle{ \lambda_k }[/math] with corresponding eigenfunctions [math]\displaystyle{ u_k }[/math]. Also since [math]\displaystyle{ u_k(0)=u_k(1)=0\quad \forall k }[/math] Each [math]\displaystyle{ u_k }[/math] can be expanded as a fourier series in terms of sine functions.
[math]\displaystyle{ u_k = \sum_{n=1}^{\infty} a_{n,k}\sin(n\pi x) }[/math]
Transforming (3) into the equivalent variational problem gives
[math]\displaystyle{ J[u] = \int_{0}^{1}\,hu'^{2}-\lambda u^{2} \, dx \quad (4) }[/math]
Substituting the fourier expansion into (4) implies J must be stationary at [math]\displaystyle{ \frac{\partial J}{\partial a_{n}}=0 \quad \forall n }[/math]
[math]\displaystyle{ \frac{\partial J}{\partial a_{n}}=\int_{0}^{1}\,hn\pi \cos(n\pi x)\sum_{m=1}^{\infty} a_{m}\cos(m\pi x)\,dx-\frac{\lambda} {2}a_{n}=0 }[/math]
By defining a vector [math]\displaystyle{ \textbf{a} = \left(a_{n}\right) }[/math] and a matrix [math]\displaystyle{ M_{(n,m)} = 2\int_{0}^{1}\,hnm\pi^{2} \cos(n\pi x)\cos(m\pi x)\,dx }[/math] we have the linear system [math]\displaystyle{ M\textbf{a} = \lambda\textbf{a} }[/math] which returns the eigenvalues and eigenvectors of equation (3), with eigenvectors [math]\displaystyle{ \textbf{a} }[/math] representing coefficient vectors of the fourier expansions of eigenfunctions.
If we now construct [math]\displaystyle{ u = \sum_{k=1}^{\infty} b_k u_k }[/math] and substitute this into equation (2) we get
[math]\displaystyle{ \sum_{k=1}^{\infty} (\kappa^{2}-\lambda_k) b_k u_k = -(b-a)\partial_xh(x)-\kappa^{2}\left((b-a)x+a\right) \quad (5) }[/math]
And defining the RHS of equation (5) as [math]\displaystyle{ f(x) }[/math], a known function, we can retrieve the coefficients [math]\displaystyle{ b_k }[/math] by integrating against [math]\displaystyle{ u_k }[/math]
[math]\displaystyle{ b_k = \frac{\int_{0}^{1}\,f u_k\,dx}{(\kappa^{2}-\lambda_k) \int_{0}^{1}\, u_{k}^{2}\,dx} }[/math]
Also to find the coefficients [math]\displaystyle{ c_n }[/math] of the fourier expansion of u are just [math]\displaystyle{ \sum_{k=1}^{\infty}a_{n,k}b_{k} }[/math] with [math]\displaystyle{ a_{n.k} }[/math] being the [math]\displaystyle{ n }[/math]th coefficient of the [math]\displaystyle{ k }[/math]th eigenfunction of the Sturm-Liouville problem.
So [math]\displaystyle{ \zeta=(b-a)x+a+\sum_{n=1}^{\infty}c_{n} \sin(n\pi x) }[/math] with [math]\displaystyle{ \zeta (0)=a \quad \zeta (1)=b }[/math] and, given [math]\displaystyle{ a }[/math] and [math]\displaystyle{ b }[/math] explicitly differentiating [math]\displaystyle{ \zeta }[/math] gives [math]\displaystyle{ \zeta'(0) }[/math] and [math]\displaystyle{ \zeta'(1) }[/math].
The aim here is to construct a matrix [math]\displaystyle{ S }[/math] such that, given [math]\displaystyle{ a }[/math] and [math]\displaystyle{ b }[/math]
[math]\displaystyle{ S \begin{pmatrix} \zeta(0) \\ \zeta(1) \end{pmatrix}=\begin{pmatrix} \zeta'(0) \\ \zeta'(1) \end{pmatrix} }[/math]
Taking [math]\displaystyle{ a=1,\,b=0 }[/math] shows that the first column of [math]\displaystyle{ S }[/math] must be [math]\displaystyle{ \begin{pmatrix} \zeta'(0) \\ \zeta'(1) \end{pmatrix} \bigg|_{(a=1,b=0)} }[/math] and likewise taking [math]\displaystyle{ a=0,\,b=1 }[/math] shows the second column must be [math]\displaystyle{ \begin{pmatrix} \zeta'(0) \\ \zeta'(1) \end{pmatrix} \bigg|_{(a=0,b=1)} }[/math]. So [math]\displaystyle{ S }[/math] is given by
[math]\displaystyle{ \begin{pmatrix} \begin{matrix} \zeta'(0) \\ \zeta'(1) \end{matrix} \bigg|_{(a=1,b=0)} \begin{matrix} \zeta'(0) \\ \zeta'(1) \end{matrix} \bigg|_{(a=0,b=1)} \end{pmatrix} }[/math]
Now for the areas of constant depth there is a potential of the form [math]\displaystyle{ e^{ikx} }[/math] which, creates reflected and transmitted potentials from the variable depth area of the form [math]\displaystyle{ Re^{ikx} }[/math] and [math]\displaystyle{ Te^{ikx} }[/math] respectively where the magnitudes of [math]\displaystyle{ R }[/math] and [math]\displaystyle{ T }[/math] are unknown. We can calculate that the boundary conditions for [math]\displaystyle{ \zeta }[/math] must be
[math]\displaystyle{ \zeta(0) = 1+R, \quad \zeta(1) = Te^{ik}, \quad \zeta'(0) = ik(1-R), \quad \zeta'(1) = ikTe^{ik} }[/math]
Knowing [math]\displaystyle{ S }[/math] these boundary conditions can be solved for [math]\displaystyle{ R }[/math] and [math]\displaystyle{ T }[/math], which in turn gives actual numerical boundary conditions to the original problem. Taking a linear combination of the solutions already calculated ([math]\displaystyle{ a=1,\,b=0 }[/math] and [math]\displaystyle{ a=0,\,b=1 }[/math]) will provide the solution for these new boundary conditions. This solution, along with the potentials outside this region gives a potential for the whole real axis.
If a waveform [math]\displaystyle{ f(x) }[/math] is travelling in from [math]\displaystyle{ -\infty }[/math] Taking the fourier transform gives a potential
[math]\displaystyle{ \hat{f}(k)=\int_{-\infty}^{\infty}\,f(x)e^{2\pi i k} \,dx }[/math]
And then inverting [math]\displaystyle{ \hat{f}(k) \zeta(x) }[/math] with respect to [math]\displaystyle{ t }[/math] gives
[math]\displaystyle{ f(x,t) = \int_{-\infty}^{\infty}\,\hat{f}(k) \zeta(x) e^{-2\pi i k t} \,dk }[/math]
which is the time dependent solution for the wave form.
