WikiWaves - User contributions [en]

Variable Depth Shallow Water Wave Equation

2009-01-09T04:09:09Z

Slawrence: /* Solution using Separation of Variables */

== 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 [[:Category:Shallow Depth|Shallow Depth]].

== Equations ==

We begin with the shallow depth equation
{{shallow depth one dimension}}

== Waves in a finite basin ==

We consider the problem of waves in a finite basin <math>0<x<1</math>. At the edge of the basin the boundary conditions are
<center>
<math>\left.\partial_x \zeta\right|_{x=0} = \left.\partial_x \zeta\right|_{x=1} =0</math>
</center>.

We solve the equations by expanding in the modes for the basin which satisfy
<center><math>
\partial_x \left(h(x) \partial_x \zeta_n \right) = -\lambda_n,
</math></center>
normalised so that
<center><math>
\int_0^1 \zeta_n \zeta_m = \delta_{mn}.
</math></center>
The solution is then given by
<center><math>
\zeta(x,t) = \sum_{n=0}^{\infty} \left(\int_0^1 \zeta_n(x^\prime) \zeta_0 (x^\prime) dx^\prime \right) \zeta_n(x) \cos(\sqrt{\lambda_n} t )
</math></center>
<center><math>
+ \sum_{n=1} ^{\infty} \left(\int_0^1 \zeta_n(x^\prime) \partial_t\zeta_0 (x^\prime) dx^\prime \right) \zeta_n(x) \frac{\sin(\sqrt{\lambda_n} t )}{\sqrt{\lambda_n}}
</math></center>
where we have assumed that <math>\lambda_0 = 0</math>.

=== Calculation of <math>\zeta_n</math> ===

We can calculate the eigenfunctions <math>\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
<center>
<math>J[\zeta] = \int_0^1 \frac{1}{2}\left\{ \left(h(x) \partial_x \zeta\right)^2 - \lambda \zeta^2 \right\}</math>
</center>
subject to the boundary conditions that the normal derivative vanishes (where <math>\lambda</math> is the eigenvalue).

We expand the displacement in the eigenfunctions for constant depth <math>h=1</math>
<center>
<math>
\zeta = \sum_{n=1}^{N} a_n \psi_n(x)
</math>
</center>
where
<center>
<math>
\psi_n = \sqrt{2} \cos( (n-1) \pi x),\,\,n\ne 1
</math>
</center>
<center>
<math>
\psi_0 = 1,\,
</math>
</center>
and substitute this expansion into the variational equation we obtain
<center>
<math>
\mathbf{M} \vec{a} = \lambda \vec{a}
</math>
</center>
where the elements of the matrix '''M''' are
<center>
<math>
M_{mn} = \int_0^1 \left\{ \left(\partial_x \psi_m h(x) \partial_x \psi_n\right) \right\}
</math>
</center>

=== Matlab code ===

== Waves in an infinite basin ==

We assume that the depth is constant and equal to one outside the region <math>0<x<1</math>.
We can therefore write the wave as

== Solution using Separation of Variables ==

Taking a separable solution gives the eigenvalue problem
<center>
<math>\partial_x \left( h(x) \partial_x\zeta \right) = -\kappa^{2}\zeta \quad (1)</math>
</center>

Given boundary conditions <math>\zeta \mid_0 = a</math> and <math>\zeta \mid_1 = b</math> we can take <math>\zeta = (b-a)x + a + u </math> With <math> u </math> satisfying <math> u |_0 = u |_1 = 0</math>

Substituting this form into (1) gives
<center>
<math> (b-a)\partial_xh(x)+\partial_x(h(x)\partial_xu) = -\kappa^{2}\left((b-a)x+a+u\right)</math>
</center>
Or, on rearranging
<center>
<math> \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>
</center>
Now consider the homogenous Sturm-Liouville problem for u
<center>
<math> \partial_x(h(x)\partial_xu)+\lambda u = 0\quad u|_0=u|_1=0 \quad (3)</math>
</center>
By Sturm-Liouville theory this has an infinite set of eigenvalues <math> \lambda_k </math> with corresponding eigenfunctions <math> u_k </math>. Also since <math> u_k|_0=u_k|_1=0\quad \forall k </math> Each <math> u_k</math> can be expanded as a fourier series in terms of sine functions.
<center>
<math> u_k = \sum_{n=1}^{\infty} a_{n,k}\sin(n\pi x)</math>
</center>
Transforming (3) into the equivalent variational problem gives
<center>
<math> J[u] = \int_{0}^{1}\,hu'^{2}-\lambda u^{2} \, dx \quad (4) </math>
</center>
Substituting the fourier expansion into (4) implies J must be stationary at <math> \frac{\partial J}{\partial a_{n}}=0 \quad \forall n </math>
<center>
<math> \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>
</center>
By defining a vector <math> \textbf{a} = \left(a_{n}\right)</math> and a matrix <math>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> M\textbf{a} = \lambda\textbf{a} </math> which returns the eigenvalues and eigenvectors of equation (3), with eigenvectors <math> \textbf{a} </math> representing coefficient vectors of the fourier expansions of eigenfunctions.

If we now construct <math> u = \sum_{k=1}^{\infty} b_k u_k </math> and substitute this into equation (2) we get
<center>
<math> \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>
</center>
And defining the RHS of equation (5) as <math> f(x) </math>, a known function, we can retrieve the coefficients <math> b_k </math> by integrating against <math> u_k </math>
<center>
<math> b_k = \frac{\int_{0}^{1}\,f u_k\,dx}{(\kappa^{2}-\lambda_k) \int_{0}^{1}\, u_{k}^{2}\,dx} </math>
</center>

Also to find the coefficients <math> c_n </math> of the fourier expansion of u are just <math> \sum_{k=1}^{\infty}a_{n,k}b_{k} </math> with <math> a_{n.k} </math> being the <math>n</math>th coefficient of the <math>k</math>th eigenfunction of the Sturm-Liouville problem.

So <math> \zeta=(b-a)x+a+\sum_{n=1}^{\infty}c_{n} \sin(n\pi x) </math> with <math> \zeta |_0=a \quad \zeta |_1=b </math> and, given <math> a </math> and <math> b </math> explicitly differentiating <math> \zeta </math> gives <math> \partial_x \zeta |_0 </math> and <math> \partial_x |zeta |_1 </math>.

The aim here is to construct a matrix <math> S </math> such that, given <math> a </math> and <math> b </math>
<center>
<math> S \begin{pmatrix} \zeta |_0 \\ \zeta |_1 \end{pmatrix}=\begin{pmatrix} \partial_x \zeta |_0 \\ \partial_x \zeta |_1 \end{pmatrix}</math>
</center>

Taking <math> a=1,\,b=0 </math> to give <math> \zeta_1 </math> shows that the first column of <math> S </math> must be <math> \begin{pmatrix} \partial_x \zeta_1 |_0 \\ \partial_x \zeta_1 |_1 \end{pmatrix} </math> and likewise taking <math> a=0,\,b=1 </math> to give <math> \zeta_2 </math> shows the second column must be <math> \begin{pmatrix} \partial_x \zeta_2 |_0 \\ \partial_x \zeta_2 |_1 \end{pmatrix} </math>. So <math>S</math> is given by
<center>
<math> \begin{pmatrix} \partial_x \zeta_1 |_0 \, \partial_x \zeta_2 |_0 \\ \partial_x \zeta_1 |_1 \, \partial_x \zeta_2 |_1 \end{pmatrix} </math>
</center>

Now for the areas of constant depth there is a potential of the form <math> e^{ikx} </math> which, creates reflected and transmitted potentials from the variable depth area of the form <math> Re^{ikx} </math> and <math> Te^{ikx} </math> respectively where the magnitudes of <math> R </math> and <math> T </math> are unknown. We can calculate that the boundary conditions for <math> \zeta </math> must be
<center>
<math> \zeta |_0 = 1+R, \quad \zeta |_1 = Te^{ik}, \quad \partial_x \zeta |_0 = ik(1-R), \quad \partial_x \zeta |_1 = ikTe^{ik} </math>
</center>

Knowing <math> S </math> these boundary conditions can be solved for <math> R</math> and <math>T</math>, which in turn gives actual numerical boundary conditions to the original problem. Taking a linear combination of the solutions already calculated (<math> \zeta_1 </math> and <math> \zeta_2 </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> f(x) </math> is travelling in from <math> -\infty </math> Taking the fourier transform gives
<center>
<math> \hat{f}(k)=\int_{-\infty}^{\infty}\,f(x)e^{2\pi i k} \,dx </math>
</center>

and then inverting <math> \hat{f}(k) \zeta(x) </math> with respect to <math> t </math> gives
<center>
<math> f(x,t) = \int_{-\infty}^{\infty}\,\hat{f}(k) \zeta(x) e^{-2\pi i k t} \,dk </math>
</center>

which is the time dependent solution for the wave form.

Variable Depth Shallow Water Wave Equation

2009-01-09T03:45:47Z

Slawrence: /* Solution using Separation of Variables */

== 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 [[:Category:Shallow Depth|Shallow Depth]].

== Equations ==

We begin with the shallow depth equation
{{shallow depth one dimension}}

== Waves in a finite basin ==

We consider the problem of waves in a finite basin <math>0<x<1</math>. At the edge of the basin the boundary conditions are
<center>
<math>\left.\partial_x \zeta\right|_{x=0} = \left.\partial_x \zeta\right|_{x=1} =0</math>
</center>.

We solve the equations by expanding in the modes for the basin which satisfy
<center><math>
\partial_x \left(h(x) \partial_x \zeta_n \right) = -\lambda_n,
</math></center>
normalised so that
<center><math>
\int_0^1 \zeta_n \zeta_m = \delta_{mn}.
</math></center>
The solution is then given by
<center><math>
\zeta(x,t) = \sum_{n=0}^{\infty} \left(\int_0^1 \zeta_n(x^\prime) \zeta_0 (x^\prime) dx^\prime \right) \zeta_n(x) \cos(\sqrt{\lambda_n} t )
</math></center>
<center><math>
+ \sum_{n=1} ^{\infty} \left(\int_0^1 \zeta_n(x^\prime) \partial_t\zeta_0 (x^\prime) dx^\prime \right) \zeta_n(x) \frac{\sin(\sqrt{\lambda_n} t )}{\sqrt{\lambda_n}}
</math></center>
where we have assumed that <math>\lambda_0 = 0</math>.

=== Calculation of <math>\zeta_n</math> ===

We can calculate the eigenfunctions <math>\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
<center>
<math>J[\zeta] = \int_0^1 \frac{1}{2}\left\{ \left(h(x) \partial_x \zeta\right)^2 - \lambda \zeta^2 \right\}</math>
</center>
subject to the boundary conditions that the normal derivative vanishes (where <math>\lambda</math> is the eigenvalue).

We expand the displacement in the eigenfunctions for constant depth <math>h=1</math>
<center>
<math>
\zeta = \sum_{n=1}^{N} a_n \psi_n(x)
</math>
</center>
where
<center>
<math>
\psi_n = \sqrt{2} \cos( (n-1) \pi x),\,\,n\ne 1
</math>
</center>
<center>
<math>
\psi_0 = 1,\,
</math>
</center>
and substitute this expansion into the variational equation we obtain
<center>
<math>
\mathbf{M} \vec{a} = \lambda \vec{a}
</math>
</center>
where the elements of the matrix '''M''' are
<center>
<math>
M_{mn} = \int_0^1 \left\{ \left(\partial_x \psi_m h(x) \partial_x \psi_n\right) \right\}
</math>
</center>

=== Matlab code ===

== Waves in an infinite basin ==

We assume that the depth is constant and equal to one outside the region <math>0<x<1</math>.
We can therefore write the wave as

== Solution using Separation of Variables ==

Taking a separable solution gives the eigenvalue problem
<center>
<math>\partial_x \left( h(x) \partial_x\zeta \right) = -\kappa^{2}\zeta \quad (1)</math>
</center>

Given boundary conditions <math>\zeta (0) = a</math> and <math>\zeta (1) = b</math> we can take <math>\zeta = (b-a)x + a + u </math> With <math> u </math> satisfying <math> u (0) = u (1) = 0</math>

Substituting this form into (1) gives
<center>
<math> (b-a)\partial_xh(x)+\partial_x(h(x)\partial_xu) = -\kappa^{2}\left((b-a)x+a+u\right)</math>
</center>
Or, on rearranging
<center>
<math> \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>
</center>
Now consider the homogenous Sturm-Liouville problem for u
<center>
<math> \partial_x(h(x)\partial_xu)+\lambda u = 0\quad u(0)=u(1)=0 \quad (3)</math>
</center>
By Sturm-Liouville theory this has an infinite set of eigenvalues <math> \lambda_k </math> with corresponding eigenfunctions <math> u_k </math>. Also since <math> u_k(0)=u_k(1)=0\quad \forall k </math> Each <math> u_k</math> can be expanded as a fourier series in terms of sine functions.
<center>
<math> u_k = \sum_{n=1}^{\infty} a_{n,k}\sin(n\pi x)</math>
</center>
Transforming (3) into the equivalent variational problem gives
<center>
<math> J[u] = \int_{0}^{1}\,hu'^{2}-\lambda u^{2} \, dx \quad (4) </math>
</center>
Substituting the fourier expansion into (4) implies J must be stationary at <math> \frac{\partial J}{\partial a_{n}}=0 \quad \forall n </math>
<center>
<math> \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>
</center>
By defining a vector <math> \textbf{a} = \left(a_{n}\right)</math> and a matrix <math>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> M\textbf{a} = \lambda\textbf{a} </math> which returns the eigenvalues and eigenvectors of equation (3), with eigenvectors <math> \textbf{a} </math> representing coefficient vectors of the fourier expansions of eigenfunctions.

If we now construct <math> u = \sum_{k=1}^{\infty} b_k u_k </math> and substitute this into equation (2) we get
<center>
<math> \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>
</center>
And defining the RHS of equation (5) as <math> f(x) </math>, a known function, we can retrieve the coefficients <math> b_k </math> by integrating against <math> u_k </math>
<center>
<math> b_k = \frac{\int_{0}^{1}\,f u_k\,dx}{(\kappa^{2}-\lambda_k) \int_{0}^{1}\, u_{k}^{2}\,dx} </math>
</center>

Also to find the coefficients <math> c_n </math> of the fourier expansion of u are just <math> \sum_{k=1}^{\infty}a_{n,k}b_{k} </math> with <math> a_{n.k} </math> being the <math>n</math>th coefficient of the <math>k</math>th eigenfunction of the Sturm-Liouville problem.

So <math> \zeta=(b-a)x+a+\sum_{n=1}^{\infty}c_{n} \sin(n\pi x) </math> with <math> \zeta (0)=a \quad \zeta (1)=b </math> and, given <math> a </math> and <math> b </math> explicitly differentiating <math> \zeta </math> gives <math> \zeta'(0) </math> and <math> \zeta'(1) </math>.

The aim here is to construct a matrix <math> S </math> such that, given <math> a </math> and <math> b </math>
<center>
<math> S \begin{pmatrix} \zeta(0) \\ \zeta(1) \end{pmatrix}=\begin{pmatrix} \zeta'(0) \\ \zeta'(1) \end{pmatrix}</math>
</center>

Taking <math> a=1,\,b=0 </math> shows that the first column of <math> S </math> must be <math> \begin{pmatrix} \zeta'(0) \\ \zeta'(1) \end{pmatrix} \bigg|_{(a=1,b=0)} </math> and likewise taking <math> a=0,\,b=1 </math> shows the second column must be <math> \begin{pmatrix} \zeta'(0) \\ \zeta'(1) \end{pmatrix} \bigg|_{(a=0,b=1)} </math>. So <math>S</math> is given by
<center>
<math> \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>
</center>

Now for the areas of constant depth there is a potential of the form <math> e^{ikx} </math> which, creates reflected and transmitted potentials from the variable depth area of the form <math> Re^{ikx} </math> and <math> Te^{ikx} </math> respectively where the magnitudes of <math> R </math> and <math> T </math> are unknown. We can calculate that the boundary conditions for <math> \zeta </math> must be
<center>
<math> \zeta(0) = 1+R, \quad \zeta(1) = Te^{ik}, \quad \zeta'(0) = ik(1-R), \quad \zeta'(1) = ikTe^{ik} </math>
</center>

Knowing <math> S </math> these boundary conditions can be solved for <math> R</math> and <math>T</math>, which in turn gives actual numerical boundary conditions to the original problem. Taking a linear combination of the solutions already calculated (<math> a=1,\,b=0 </math> and <math> 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> f(x) </math> is travelling in from <math> -\infty </math> Taking the fourier transform gives
<center>
<math> \hat{f}(k)=\int_{-\infty}^{\infty}\,f(x)e^{2\pi i k} \,dx </math>
</center>

and then inverting <math> \hat{f}(k) \zeta(x) </math> with respect to <math> t </math> gives
<center>
<math> f(x,t) = \int_{-\infty}^{\infty}\,\hat{f}(k) \zeta(x) e^{-2\pi i k t} \,dk </math>
</center>

which is the time dependent solution for the wave form.

Variable Depth Shallow Water Wave Equation

2009-01-09T03:42:38Z

Slawrence: /* Solution using Separation of Variables */

== 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 [[:Category:Shallow Depth|Shallow Depth]].

== Equations ==

We begin with the shallow depth equation
{{shallow depth one dimension}}

== Waves in a finite basin ==

We consider the problem of waves in a finite basin <math>0<x<1</math>. At the edge of the basin the boundary conditions are
<center>
<math>\left.\partial_x \zeta\right|_{x=0} = \left.\partial_x \zeta\right|_{x=1} =0</math>
</center>.

We solve the equations by expanding in the modes for the basin which satisfy
<center><math>
\partial_x \left(h(x) \partial_x \zeta_n \right) = -\lambda_n,
</math></center>
normalised so that
<center><math>
\int_0^1 \zeta_n \zeta_m = \delta_{mn}.
</math></center>
The solution is then given by
<center><math>
\zeta(x,t) = \sum_{n=0}^{\infty} \left(\int_0^1 \zeta_n(x^\prime) \zeta_0 (x^\prime) dx^\prime \right) \zeta_n(x) \cos(\sqrt{\lambda_n} t )
</math></center>
<center><math>
+ \sum_{n=1} ^{\infty} \left(\int_0^1 \zeta_n(x^\prime) \partial_t\zeta_0 (x^\prime) dx^\prime \right) \zeta_n(x) \frac{\sin(\sqrt{\lambda_n} t )}{\sqrt{\lambda_n}}
</math></center>
where we have assumed that <math>\lambda_0 = 0</math>.

=== Calculation of <math>\zeta_n</math> ===

We can calculate the eigenfunctions <math>\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
<center>
<math>J[\zeta] = \int_0^1 \frac{1}{2}\left\{ \left(h(x) \partial_x \zeta\right)^2 - \lambda \zeta^2 \right\}</math>
</center>
subject to the boundary conditions that the normal derivative vanishes (where <math>\lambda</math> is the eigenvalue).

We expand the displacement in the eigenfunctions for constant depth <math>h=1</math>
<center>
<math>
\zeta = \sum_{n=1}^{N} a_n \psi_n(x)
</math>
</center>
where
<center>
<math>
\psi_n = \sqrt{2} \cos( (n-1) \pi x),\,\,n\ne 1
</math>
</center>
<center>
<math>
\psi_0 = 1,\,
</math>
</center>
and substitute this expansion into the variational equation we obtain
<center>
<math>
\mathbf{M} \vec{a} = \lambda \vec{a}
</math>
</center>
where the elements of the matrix '''M''' are
<center>
<math>
M_{mn} = \int_0^1 \left\{ \left(\partial_x \psi_m h(x) \partial_x \psi_n\right) \right\}
</math>
</center>

=== Matlab code ===

== Waves in an infinite basin ==

We assume that the depth is constant and equal to one outside the region <math>0<x<1</math>.
We can therefore write the wave as

== Solution using Separation of Variables ==

Taking a separable solution gives the eigenvalue problem
<center>
<math>\partial_x \left( h(x) \partial_x\zeta \right) = -\kappa^{2}\zeta \quad (1)</math>
</center>

Given boundary conditions <math>\zeta (0) = a</math> and <math>\zeta (1) = b</math> we can take <math>\zeta = (b-a)x + a + u </math> With <math> u </math> satisfying <math> u (0) = u (1) = 0</math>

Substituting this form into (1) gives
<center>
<math> (b-a)\partial_xh(x)+\partial_x(h(x)\partial_xu) = -\kappa^{2}\left((b-a)x+a+u\right)</math>
</center>
Or, on rearranging
<center>
<math> \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>
</center>
Now consider the homogenous Sturm-Liouville problem for u
<center>
<math> \partial_x(h(x)\partial_xu)+\lambda u = 0\quad u(0)=u(1)=0 \quad (3)</math>
</center>
By Sturm-Liouville theory this has an infinite set of eigenvalues <math> \lambda_k </math> with corresponding eigenfunctions <math> u_k </math>. Also since <math> u_k(0)=u_k(1)=0\quad \forall k </math> Each <math> u_k</math> can be expanded as a fourier series in terms of sine functions.
<center>
<math> u_k = \sum_{n=1}^{\infty} a_{n,k}\sin(n\pi x)</math>
</center>
Transforming (3) into the equivalent variational problem gives
<center>
<math> J[u] = \int_{0}^{1}\,hu'^{2}-\lambda u^{2} \, dx \quad (4) </math>
</center>
Substituting the fourier expansion into (4) implies J must be stationary at <math> \frac{\partial J}{\partial a_{n}}=0 \quad \forall n </math>
<center>
<math> \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>
</center>
By defining a vector <math> \textbf{a} = \left(a_{n}\right)</math> and a matrix <math>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> M\textbf{a} = \lambda\textbf{a} </math> which returns the eigenvalues and eigenvectors of equation (3), with eigenvectors <math> \textbf{a} </math> representing coefficient vectors of the fourier expansions of eigenfunctions.

If we now construct <math> u = \sum_{k=1}^{\infty} b_k u_k </math> and substitute this into equation (2) we get
<center>
<math> \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>
</center>
And defining the RHS of equation (5) as <math> f(x) </math>, a known function, we can retrieve the coefficients <math> b_k </math> by integrating against <math> u_k </math>
<center>
<math> b_k = \frac{\int_{0}^{1}\,f u_k\,dx}{(\kappa^{2}-\lambda_k) \int_{0}^{1}\, u_{k}^{2}\,dx} </math>
</center>

Also to find the coefficients <math> c_n </math> of the fourier expansion of u are just <math> \sum_{k=1}^{\infty}a_{n,k}b_{k} </math> with <math> a_{n.k} </math> being the <math>n</math>th coefficient of the <math>k</math>th eigenfunction of the Sturm-Liouville problem.

So <math> \zeta=(b-a)x+a+\sum_{n=1}^{\infty}c_{n} \sin(n\pi x) </math> with <math> \zeta (0)=a \quad \zeta (1)=b </math> and, given <math> a </math> and <math> b </math> explicitly differentiating <math> \zeta </math> gives <math> \zeta'(0) </math> and <math> \zeta'(1) </math>.

The aim here is to construct a matrix <math> S </math> such that, given <math> a </math> and <math> b </math>
<center>
<math> S \begin{pmatrix} \zeta(0) \\ \zeta(1) \end{pmatrix}=\begin{pmatrix} \zeta'(0) \\ \zeta'(1) \end{pmatrix}</math>
</center>

Taking <math> a=1,\,b=0 </math> shows that the first column of <math> S </math> must be <math> \begin{pmatrix} \zeta'(0) \\ \zeta'(1) \end{pmatrix} \bigg|_{(a=1,b=0)} </math> and likewise taking <math> a=0,\,b=1 </math> shows the second column must be <math> \begin{pmatrix} \zeta'(0) \\ \zeta'(1) \end{pmatrix} \bigg|_{(a=0,b=1)} </math>. So <math>S</math> is given by
<center>
<math> \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>
</center>

Now for the areas of constant depth there is a potential of the form <math> e^{ikx} </math> which, creates reflected and transmitted potentials from the variable depth area of the form <math> Re^{ikx} </math> and <math> Te^{ikx} </math> respectively where the magnitudes of <math> R </math> and <math> T </math> are unknown. We can calculate that the boundary conditions for <math> \zeta </math> must be
<center>
<math> \zeta(0) = 1+R, \quad \zeta(1) = Te^{ik}, \quad \zeta'(0) = ik(1-R), \quad \zeta'(1) = ikTe^{ik} </math>
</center>

Knowing <math> S </math> these boundary conditions can be solved for <math> R</math> and <math>T</math>, which in turn gives actual numerical boundary conditions to the original problem. Taking a linear combination of the solutions already calculated (<math> a=1,\,b=0 </math> and <math> 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> f(x) </math> is travelling in from <math> -\infty </math> Taking the fourier transform gives a potential
<center>
<math> \hat{f}(k)=\int_{-\infty}^{\infty}\,f(x)e^{2\pi i k} \,dx </math>
</center>

And then inverting <math> \hat{f}(k) \zeta(x) </math> with respect to <math> t </math> gives
<center>
<math> f(x,t) = \int_{-\infty}^{\infty}\,\hat{f}(k) \zeta(x) e^{-2\pi i k t} \,dk </math>
</center>

which is the time dependent solution for the wave form.

Variable Depth Shallow Water Wave Equation

2009-01-09T03:26:19Z

Slawrence: /* Calculation of <math>\zeta_n</math> */

== 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 [[:Category:Shallow Depth|Shallow Depth]].

== Equations ==

We begin with the shallow depth equation
{{shallow depth one dimension}}

== Waves in a finite basin ==

We consider the problem of waves in a finite basin <math>0<x<1</math>. At the edge of the basin the boundary conditions are
<center>
<math>\left.\partial_x \zeta\right|_{x=0} = \left.\partial_x \zeta\right|_{x=1} =0</math>
</center>.

We solve the equations by expanding in the modes for the basin which satisfy
<center><math>
\partial_x \left(h(x) \partial_x \zeta_n \right) = -\lambda_n,
</math></center>
normalised so that
<center><math>
\int_0^1 \zeta_n \zeta_m = \delta_{mn}.
</math></center>
The solution is then given by
<center><math>
\zeta(x,t) = \sum_{n=0}^{\infty} \left(\int_0^1 \zeta_n(x^\prime) \zeta_0 (x^\prime) dx^\prime \right) \zeta_n(x) \cos(\sqrt{\lambda_n} t )
</math></center>
<center><math>
+ \sum_{n=1} ^{\infty} \left(\int_0^1 \zeta_n(x^\prime) \partial_t\zeta_0 (x^\prime) dx^\prime \right) \zeta_n(x) \frac{\sin(\sqrt{\lambda_n} t )}{\sqrt{\lambda_n}}
</math></center>
where we have assumed that <math>\lambda_0 = 0</math>.

=== Calculation of <math>\zeta_n</math> ===

We can calculate the eigenfunctions <math>\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
<center>
<math>J[\zeta] = \int_0^1 \frac{1}{2}\left\{ \left(h(x) \partial_x \zeta\right)^2 - \lambda \zeta^2 \right\}</math>
</center>
subject to the boundary conditions that the normal derivative vanishes (where <math>\lambda</math> is the eigenvalue).

We expand the displacement in the eigenfunctions for constant depth <math>h=1</math>
<center>
<math>
\zeta = \sum_{n=1}^{N} a_n \psi_n(x)
</math>
</center>
where
<center>
<math>
\psi_n = \sqrt{2} \cos( (n-1) \pi x),\,\,n\ne 1
</math>
</center>
<center>
<math>
\psi_0 = 1,\,
</math>
</center>
and substitute this expansion into the variational equation we obtain
<center>
<math>
\mathbf{M} \vec{a} = \lambda \vec{a}
</math>
</center>
where the elements of the matrix '''M''' are
<center>
<math>
M_{mn} = \int_0^1 \left\{ \left(\partial_x \psi_m h(x) \partial_x \psi_n\right) \right\}
</math>
</center>

=== Matlab code ===

== Waves in an infinite basin ==

We assume that the depth is constant and equal to one outside the region <math>0<x<1</math>.
We can therefore write the wave as

== Solution using Separation of Variables ==

Taking a separable solution gives the eigenvalue problem
<center>
<math>\partial_x \left( h(x) \partial_x\zeta \right) = -\kappa^{2}\zeta \quad (1)</math>
</center>

Given boundary conditions <math>\zeta (0) = a</math> and <math>\zeta (1) = b</math> we can take <math>\zeta = (b-a)x + a + u </math> With <math> u </math> satisfying <math> u (0) = u (1) = 0</math>

Substituting this form into (1) gives
<center>
<math> (b-a)\partial_xh(x)+\partial_x(h(x)\partial_xu) = -\kappa^{2}\left((b-a)x+a+u\right)</math>
</center>
Or, on rearranging
<center>
<math> \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>
</center>
Now consider the homogenous Sturm-Liouville problem for u
<center>
<math> \partial_x(h(x)\partial_xu)+\lambda u = 0\quad u(0)=u(1)=0 \quad (3)</math>
</center>
By Sturm-Liouville theory this has an infinite set of eigenvalues <math> \lambda_k </math> with corresponding eigenfunctions <math> u_k </math>. Also since <math> u_k(0)=u_k(1)=0\quad \forall k </math> Each <math> u_k</math> can be expanded as a fourier series in terms of sine functions.
<center>
<math> u_k = \sum_{n=1}^{\infty} a_{n,k}\sin(n\pi x)</math>
</center>
Transforming (3) into the equivalent variational problem gives
<center>
<math> J[u] = \int_{0}^{1}\,hu'^{2}-\lambda u^{2} \, dx \quad (4) </math>
</center>
Substituting the fourier expansion into (4) implies J must be stationary at <math> \frac{\partial J}{\partial a_{n}}=0 \quad \forall n </math>
<center>
<math> \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>
</center>
By defining a vector <math> \textbf{a} = \left(a_{n}\right)</math> and a matrix <math>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> M\textbf{a} = \lambda\textbf{a} </math> which returns the eigenvalues and eigenvectors of equation (3), with eigenvectors <math> \textbf{a} </math> representing coefficient vectors of the fourier expansions of eigenfunctions.

If we now construct <math> u = \sum_{k=1}^{\infty} b_k u_k </math> and substitute this into equation (2) we get
<center>
<math> \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>
</center>
And defining the RHS of equation (5) as <math> f(x) </math>, a known function, we can retrieve the coefficients <math> b_k </math> by integrating against <math> u_k </math>
<center>
<math> b_k = \frac{\int_{0}^{1}\,f u_k\,dx}{(\kappa^{2}-\lambda_k) \int_{0}^{1}\, u_{k}^{2}\,dx} </math>
</center>

Also to find the coefficients <math> c_n </math> of the fourier expansion of u are just <math> \sum_{k=1}^{\infty}a_{n,k}b_{k} </math> with <math> a_{n.k} </math> being the <math>n</math>th coefficient of the <math>k</math>th eigenfunction of the Sturm-Liouville problem.

So <math> \zeta=(b-a)x+a+\sum_{n=1}^{\infty}c_{n} \sin(n\pi x) </math> with <math> \zeta (0)=a \quad \zeta (1)=b </math> and, given <math> a </math> and <math> b </math> explicitly differentiating <math> \zeta </math> gives <math> \zeta'(0) </math> and <math> \zeta'(1) </math>.

The aim here is to construct a matrix <math> S </math> such that, given <math> a </math> and <math> b </math>
<center>
<math> S \begin{pmatrix} \zeta(0) \\ \zeta(1) \end{pmatrix}=\begin{pmatrix} \zeta'(0) \\ \zeta'(1) \end{pmatrix}</math>
</center>

Taking <math> a=1,\,b=0 </math> shows that the first column of <math> S </math> must be <math> \begin{pmatrix} \zeta'(0) \\ \zeta'(1) \end{pmatrix} \bigg|_{(a=1,b=0)} </math> and likewise taking <math> a=0,\,b=1 </math> shows the second column must be <math> \begin{pmatrix} \zeta'(0) \\ \zeta'(1) \end{pmatrix} \bigg|_{(a=0,b=1)} </math>. So <math>S</math> is given by
<center>
<math> \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>
</center>

Now with a potential of the form <math> e^{ikx} </math> which, creates reflected and transmitted potentials from the variable depth area of the form <math> Re^{ikx} </math> and <math> Te^{ikx} </math> respectively where the magnitudes of R and T are unknown. We can calculate that the boundary conditions for <math> \zeta </math> must be
<center>
<math> \zeta(0) = 1+R, \quad \zeta(1) = Te^{ik}, \quad \zeta'(0) = ik(1-R), \quad \zeta'(1) = ikTe^{ik} </math>
</center>

Knowing <math> S </math> these boundary conditions can be solved for <math> R</math> and <math>T</math>, which in turn gives actual numerical boundary conditions to the original problem. Taking a linear combination of the solutions already calculated (<math> a=1,\,b=0 </math> and <math> 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> f(x) </math> is travelling in from <math> -\infty </math> Taking the fourier transform gives a potential
<center>
<math> \hat{f}(k)=\int_{-\infty}^{\infty}\,f(x)e^{2\pi i k} \,dx </math>
</center>

And then inverting <math> \hat{f}(k) \zeta(x) </math> with respect to <math> t </math> gives
<center>
<math> f(x,t) = \int_{-\infty}^{\infty}\,\hat{f}(k) \zeta(x) e^{-2\pi i k t} \,dk </math>
</center>

which is the time dependent solution.

Variable Depth Shallow Water Wave Equation

2009-01-09T03:24:45Z

Slawrence: /* Solution using Separation of Variables */

== 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 [[:Category:Shallow Depth|Shallow Depth]].

== Equations ==

We begin with the shallow depth equation
{{shallow depth one dimension}}

== Waves in a finite basin ==

We consider the problem of waves in a finite basin <math>0<x<1</math>. At the edge of the basin the boundary conditions are
<center>
<math>\left.\partial_x \zeta\right|_{x=0} = \left.\partial_x \zeta\right|_{x=1} =0</math>
</center>.

We solve the equations by expanding in the modes for the basin which satisfy
<center><math>
\partial_x \left(h(x) \partial_x \zeta_n \right) = -\lambda_n,
</math></center>
normalised so that
<center><math>
\int_0^1 \zeta_n \zeta_m = \delta_{mn}.
</math></center>
The solution is then given by
<center><math>
\zeta(x,t) = \sum_{n=0}^{\infty} \left(\int_0^1 \zeta_n(x^\prime) \zeta_0 (x^\prime) dx^\prime \right) \zeta_n(x) \cos(\sqrt{\lambda_n} t )
</math></center>
<center><math>
+ \sum_{n=1} ^{\infty} \left(\int_0^1 \zeta_n(x^\prime) \partial_t\zeta_0 (x^\prime) dx^\prime \right) \zeta_n(x) \frac{\sin(\sqrt{\lambda_n} t )}{\sqrt{\lambda_n}}
</math></center>
where we have assumed that <math>\lambda_0 = 0</math>.

=== Calculation of <math>\zeta_n</math> ===

We can calculate the eigenfunctions <math>\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
<center>
<math>J[\zeta] = \int_0^1 \frac{1}{2}\left\{ \left(h(x) \partial_x \zeta\right)^2 - \lambda \zeta^2 \right\}</math>
</center>
subject to the boundary conditions that the normal derivative vanishes (where <math>\lambda</math> is the eigenvalue).

We expand the displacement in the eigenfunctions for constant depth <math>h=1</math>
<center>
<math>
\zeta = \sum_{n=1}^{N} a_n \psi_n(x)
</math>
</center>
where
<center>
<math>
\psi_n = \sqrt{2} \cos( (n-1) \pi x),\,\,n\ne 1
</math>
</center>
<center>
<math>
\psi_0 = 1,\,
</math>
</center>
and substitute this expansion into the variational equation we obtain
<center>
<math>
\mathbf{M} \vec{a} = \lambda \vec{a}
</math>
</center>
where the elements of the matrix '''M''' are
<center>
<math>
m_{mn} = \int_0^1 \left\{ \left(\partial_x \psi_m h(x) \partial_x \psi_n\right) \right\}
</math>
</center>

=== Matlab code ===

== Waves in an infinite basin ==

We assume that the depth is constant and equal to one outside the region <math>0<x<1</math>.
We can therefore write the wave as

== Solution using Separation of Variables ==

Taking a separable solution gives the eigenvalue problem
<center>
<math>\partial_x \left( h(x) \partial_x\zeta \right) = -\kappa^{2}\zeta \quad (1)</math>
</center>

Given boundary conditions <math>\zeta (0) = a</math> and <math>\zeta (1) = b</math> we can take <math>\zeta = (b-a)x + a + u </math> With <math> u </math> satisfying <math> u (0) = u (1) = 0</math>

Substituting this form into (1) gives
<center>
<math> (b-a)\partial_xh(x)+\partial_x(h(x)\partial_xu) = -\kappa^{2}\left((b-a)x+a+u\right)</math>
</center>
Or, on rearranging
<center>
<math> \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>
</center>
Now consider the homogenous Sturm-Liouville problem for u
<center>
<math> \partial_x(h(x)\partial_xu)+\lambda u = 0\quad u(0)=u(1)=0 \quad (3)</math>
</center>
By Sturm-Liouville theory this has an infinite set of eigenvalues <math> \lambda_k </math> with corresponding eigenfunctions <math> u_k </math>. Also since <math> u_k(0)=u_k(1)=0\quad \forall k </math> Each <math> u_k</math> can be expanded as a fourier series in terms of sine functions.
<center>
<math> u_k = \sum_{n=1}^{\infty} a_{n,k}\sin(n\pi x)</math>
</center>
Transforming (3) into the equivalent variational problem gives
<center>
<math> J[u] = \int_{0}^{1}\,hu'^{2}-\lambda u^{2} \, dx \quad (4) </math>
</center>
Substituting the fourier expansion into (4) implies J must be stationary at <math> \frac{\partial J}{\partial a_{n}}=0 \quad \forall n </math>
<center>
<math> \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>
</center>
By defining a vector <math> \textbf{a} = \left(a_{n}\right)</math> and a matrix <math>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> M\textbf{a} = \lambda\textbf{a} </math> which returns the eigenvalues and eigenvectors of equation (3), with eigenvectors <math> \textbf{a} </math> representing coefficient vectors of the fourier expansions of eigenfunctions.

If we now construct <math> u = \sum_{k=1}^{\infty} b_k u_k </math> and substitute this into equation (2) we get
<center>
<math> \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>
</center>
And defining the RHS of equation (5) as <math> f(x) </math>, a known function, we can retrieve the coefficients <math> b_k </math> by integrating against <math> u_k </math>
<center>
<math> b_k = \frac{\int_{0}^{1}\,f u_k\,dx}{(\kappa^{2}-\lambda_k) \int_{0}^{1}\, u_{k}^{2}\,dx} </math>
</center>

Also to find the coefficients <math> c_n </math> of the fourier expansion of u are just <math> \sum_{k=1}^{\infty}a_{n,k}b_{k} </math> with <math> a_{n.k} </math> being the <math>n</math>th coefficient of the <math>k</math>th eigenfunction of the Sturm-Liouville problem.

So <math> \zeta=(b-a)x+a+\sum_{n=1}^{\infty}c_{n} \sin(n\pi x) </math> with <math> \zeta (0)=a \quad \zeta (1)=b </math> and, given <math> a </math> and <math> b </math> explicitly differentiating <math> \zeta </math> gives <math> \zeta'(0) </math> and <math> \zeta'(1) </math>.

The aim here is to construct a matrix <math> S </math> such that, given <math> a </math> and <math> b </math>
<center>
<math> S \begin{pmatrix} \zeta(0) \\ \zeta(1) \end{pmatrix}=\begin{pmatrix} \zeta'(0) \\ \zeta'(1) \end{pmatrix}</math>
</center>

Taking <math> a=1,\,b=0 </math> shows that the first column of <math> S </math> must be <math> \begin{pmatrix} \zeta'(0) \\ \zeta'(1) \end{pmatrix} \bigg|_{(a=1,b=0)} </math> and likewise taking <math> a=0,\,b=1 </math> shows the second column must be <math> \begin{pmatrix} \zeta'(0) \\ \zeta'(1) \end{pmatrix} \bigg|_{(a=0,b=1)} </math>. So <math>S</math> is given by
<center>
<math> \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>
</center>

Now with a potential of the form <math> e^{ikx} </math> which, creates reflected and transmitted potentials from the variable depth area of the form <math> Re^{ikx} </math> and <math> Te^{ikx} </math> respectively where the magnitudes of R and T are unknown. We can calculate that the boundary conditions for <math> \zeta </math> must be
<center>
<math> \zeta(0) = 1+R, \quad \zeta(1) = Te^{ik}, \quad \zeta'(0) = ik(1-R), \quad \zeta'(1) = ikTe^{ik} </math>
</center>

Knowing <math> S </math> these boundary conditions can be solved for <math> R</math> and <math>T</math>, which in turn gives actual numerical boundary conditions to the original problem. Taking a linear combination of the solutions already calculated (<math> a=1,\,b=0 </math> and <math> 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> f(x) </math> is travelling in from <math> -\infty </math> Taking the fourier transform gives a potential
<center>
<math> \hat{f}(k)=\int_{-\infty}^{\infty}\,f(x)e^{2\pi i k} \,dx </math>
</center>

And then inverting <math> \hat{f}(k) \zeta(x) </math> with respect to <math> t </math> gives
<center>
<math> f(x,t) = \int_{-\infty}^{\infty}\,\hat{f}(k) \zeta(x) e^{-2\pi i k t} \,dk </math>
</center>

which is the time dependent solution.

Variable Depth Shallow Water Wave Equation

2009-01-09T01:50:56Z

Slawrence: /* Solution using Separation of Variables */

== 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 [[:Category:Shallow Depth|Shallow Depth]].

== Equations ==

We begin with the shallow depth equation
{{shallow depth one dimension}}

== Waves in a finite basin ==

We consider the problem of waves in a finite basin <math>0<x<1</math>. At the edge of the basin the boundary conditions are
<center>
<math>\left.\partial_x \zeta\right|_{x=0} = \left.\partial_x \zeta\right|_{x=1} =0</math>
</center>.

We solve the equations by expanding in the modes for the basin which satisfy
<center><math>
\partial_x \left(h(x) \partial_x \zeta_n \right) = -\lambda_n,
</math></center>
normalised so that
<center><math>
\int_0^1 \zeta_n \zeta_m = \delta_{mn}.
</math></center>
The solution is then given by
<center><math>
\zeta(x,t) = \sum_{n=0}^{\infty} \left(\int_0^1 \zeta_n(x^\prime) \zeta_0 (x^\prime) dx^\prime \right) \zeta_n(x) \cos(\sqrt{\lambda_n} t )
</math></center>
<center><math>
+ \sum_{n=1} ^{\infty} \left(\int_0^1 \zeta_n(x^\prime) \partial_t\zeta_0 (x^\prime) dx^\prime \right) \zeta_n(x) \frac{\sin(\sqrt{\lambda_n} t )}{\sqrt{\lambda_n}}
</math></center>
where we have assumed that <math>\lambda_0 = 0</math>.

=== Calculation of <math>\zeta_n</math> ===

We can calculate the eigenfunctions <math>\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
<center>
<math>J[\zeta] = \int_0^1 \frac{1}{2}\left\{ \left(h(x) \partial_x \zeta\right)^2 - \lambda \zeta^2 \right\}</math>
</center>
subject to the boundary conditions that the normal derivative vanishes (where <math>\lambda</math> is the eigenvalue).

We expand the displacement in the eigenfunctions for constant depth <math>h=1</math>
<center>
<math>
\zeta = \sum_{n=1}^{N} a_n \psi_n(x)
</math>
</center>
where
<center>
<math>
\psi_n = \sqrt{2} \cos( (n-1) \pi x),\,\,n\ne 1
</math>
</center>
<center>
<math>
\psi_0 = 1,\,
</math>
</center>
and substitute this expansion into the variational equation we obtain
<center>
<math>
\mathbf{M} \vec{a} = \lambda \vec{a}
</math>
</center>
where the elements of the matrix '''M''' are
<center>
<math>
m_{mn} = \int_0^1 \left\{ \left(\partial_x \psi_m h(x) \partial_x \psi_n\right) \right\}
</math>
</center>

=== Matlab code ===

== Waves in an infinite basin ==

We assume that the depth is constant and equal to one outside the region <math>0<x<1</math>.
We can therefore write the wave as

== Solution using Separation of Variables ==

Taking a separable solution <math>\ w (x,t) = \Tau (t) \hat{w} (x)</math> gives the eigenvalue problem
<center>
<math>\partial_x \left( h(x) \partial_x\hat{w} \right) = -\kappa^{2}\hat{w} \quad (1)</math>
</center>

Given boundary conditions <math>\hat{w} (0) = a</math> and <math>\hat{w} (1) = b</math> we can take <math>\hat{w} = (b-a)x + a + u </math> With <math> u </math> satisfying <math> u (0) = u (1) = 0</math>

Substituting this form into (1) gives
<center>
<math> (b-a)\partial_xh(x)+\partial_x(h(x)\partial_xu) = -\kappa^{2}\left((b-a)x+a+u\right)</math>
</center>
Or, on rearranging
<center>
<math> \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>
</center>
Now consider the homogenous Sturm-Liouville problem for u
<center>
<math> \partial_x(h(x)\partial_xu)+\lambda u = 0\quad u(0)=u(1)=0 \quad (3)</math>
</center>
By Sturm-Liouville theory this has an infinite set of eigenvalues <math> \lambda_k </math> with corresponding eigenfunctions <math> u_k </math>. Also since <math> u_k(0)=u_k(1)=0\quad \forall k </math> Each <math> u_k</math> can be expanded as a fourier series in terms of sine functions.
<center>
<math> u_k = \sum_{n=1}^{\infty} a_{n,k}\sin(n\pi x)</math>
</center>
Transforming (3) into the equivalent variational problem gives
<center>
<math> J[u] = \int_{0}^{1}\,hu'^{2}-\lambda u^{2} \, dx \quad (4) </math>
</center>
Substituting the fourier expansion into (4) implies J must be stationary at <math> \frac{\partial J}{\partial a_{n}}=0 \quad \forall n </math>
<center>
<math> \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>
</center>
By defining a vector <math> \textbf{a} = \left(a_{n}\right)</math> and a matrix <math>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> M\textbf{a} = \lambda\textbf{a} </math> which returns the eigenvalues and eigenvectors of equation (3), with eigenvectors <math> \textbf{a} </math> representing coefficient vectors of the fourier expansions of eigenfunctions.

If we now construct <math> u = \sum_{k=1}^{\infty} b_k u_k </math> and substitute this into equation (2) we get
<center>
<math> \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>
</center>
And defining the RHS of equation (5) as <math> f(x) </math>, a known function, we can retrieve the coefficients <math> b_k </math> by integrating against <math> u_k </math>
<center>
<math> b_k = \frac{\int_{0}^{1}\,f u_k\,dx}{(\kappa^{2}-\lambda_k) \int_{0}^{1}\, u_{k}^{2}\,dx} </math>
</center>

Also to find the coefficients <math> c_n </math> of the fourier expansion of u are just <math> \sum_{k=1}^{\infty}a_{n,k}b_{k} </math> with <math> a_{n.k} </math> being the <math>n</math>th coefficient of the <math>k</math>th eigenfunction of the Sturm-Liouville problem.

So <math> w=(b-a)x+a+\sum_{n=1}^{\infty}c_{n} \sin(n\pi x) </math> with <math> w(0)=a \quad w(1)=b </math> and, given <math> a </math> and <math> b </math> explicitly differentiating <math> w </math> gives <math> w'(0) </math> and <math> w'(1) </math>.

The aim here is to construct a matrix <math> S </math> such that, given <math> a </math> and <math> b </math>
<center>
<math> S \begin{pmatrix} w(0) \\ w(1) \end{pmatrix}=\begin{pmatrix} w'(0) \\ w'(1) \end{pmatrix}</math>
</center>

Taking <math> a=1,\,b=0 </math> shows that the first column of <math> S </math> must be <math> \begin{pmatrix} w'(0) \\ w'(1) \end{pmatrix} \bigg|_{(a=1,b=0)} </math> and likewise taking <math> a=0,\,b=1 </math> shows the second column must be <math> \begin{pmatrix} w'(0) \\ w'(1) \end{pmatrix} \bigg|_{(a=0,b=1)} </math>. So <math>S</math> is given by
<center>
<math> \begin{pmatrix} \begin{matrix} w'(0) \\ w'(1) \end{matrix} \bigg|_{(a=1,b=0)} \begin{matrix} w'(0) \\ w'(1) \end{matrix} \bigg|_{(a=0,b=1)} \end{pmatrix} </math>
</center>

Now with a potential of the form <math> e^{ikx} </math> which, creates reflected and transmitted potentials from the variable depth area of the form <math> Re^{ikx} </math> and <math> Te^{ikx} </math> respectively where the magnitudes of R and T are unknown. We can calculate that the boundary conditions for <math> w </math> must be
<center>
<math> w(0) = 1+R, \quad w(1) = Te^{ik}, \quad w'(0) = ik(1-R), \quad w'(1) = ikTe^{ik} </math>
</center>

Knowing <math> S </math> these boundary conditions can be solved for <math> R</math> and <math>T</math>, which in turn gives actual numerical boundary conditions to the original problem. Taking a linear combination of the solutions already calculated (<math> a=1,\,b=0 </math> and <math> 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 line

If a waveform <math> f(x) </math> is travelling in from <math> -\infty </math> Taking the fourier transform gives
<center>
<math> \hat{f}(k)=\int_{-\infty}^{\infty}\,f(x)e^{2\pi i k} \,dx </math>
</center>

And then
<center>
<math> f(x,t) = \int_{-\infty}^{\infty}\,\hat{f}(k) w(x) e^{-2\pi i k t} \,dk </math>

Variable Depth Shallow Water Wave Equation

2009-01-08T02:36:45Z

Slawrence: /* Solution using Separation of Variables */

== 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 [[:Category:Shallow Depth|Shallow Depth]].

== Equations ==

We begin with the shallow depth equation
{{shallow depth one dimension}}

== Waves in a finite basin ==

We consider the problem of waves in a finite basin <math>0<x<1</math>. At the edge of the basin the boundary conditions are
<center>
<math>\left.\partial_x \zeta\right|_{x=0} = \left.\partial_x \zeta\right|_{x=1} =0</math>
</center>.

We solve the equations by expanding in the modes for the basin which satisfy
<center><math>
\partial_x \left(h(x) \partial_x \zeta_n \right) = -\lambda_n,
</math></center>
normalised so that
<center><math>
\int_0^1 \zeta_n \zeta_m = \delta_{mn}.
</math></center>
The solution is then given by
<center><math>
\zeta(x,t) = \sum_{n=0}^{\infty} \left(\int_0^1 \zeta_n(x^\prime) \zeta_0 (x^\prime) dx^\prime \right) \zeta_n(x) \cos(\sqrt{\lambda_n} t )
</math></center>
<center><math>
+ \sum_{n=1} ^{\infty} \left(\int_0^1 \zeta_n(x^\prime) \partial_t\zeta_0 (x^\prime) dx^\prime \right) \zeta_n(x) \frac{\sin(\sqrt{\lambda_n} t )}{\sqrt{\lambda_n}}
</math></center>
where we have assumed that <math>\lambda_0 = 0</math>.

=== Calculation of <math>\zeta_n</math> ===

We can calculate the eigenfunctions <math>\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
<center>
<math>J[\zeta] = \int_0^1 \frac{1}{2}\left\{ \left(h(x) \partial_x \zeta\right)^2 - \lambda \zeta^2 \right\}</math>
</center>
subject to the boundary conditions that the normal derivative vanishes (where <math>\lambda</math> is the eigenvalue).

We expand the displacement in the eigenfunctions for constant depth <math>h=1</math>
<center>
<math>
\zeta = \sum_{n=1}^{N} a_n \psi_n(x)
</math>
</center>
where
<center>
<math>
\psi_n = \sqrt{2} \cos( (n-1) \pi x),\,\,n\ne 1
</math>
</center>
<center>
<math>
\psi_0 = 1,\,
</math>
</center>
and substitute this expansion into the variational equation we obtain
<center>
<math>
\mathbf{M} \vec{a} = \lambda \vec{a}
</math>
</center>
where the elements of the matrix '''M''' are
<center>
<math>
m_{mn} = \int_0^1 \left\{ \left(\partial_x \psi_m h(x) \partial_x \psi_n\right) \right\}
</math>
</center>

=== Matlab code ===

== Waves in an infinite basin ==

We assume that the depth is constant and equal to one outside the region <math>0<x<1</math>.
We can therefore write the wave as

== Solution using Separation of Variables ==

Taking a separable solution <math>\ w (x,t) = \Tau (t) \hat{w} (x)</math> gives the eigenvalue problem
<center>
<math>\partial_x \left( h(x) \partial_x\hat{w} \right) = -\kappa^{2}\hat{w} \quad (1)</math>
</center>

Given boundary conditions <math>\hat{w} (0) = a</math> and <math>\hat{w} (1) = b</math> we can take <math>\hat{w} = (b-a)x + a + u </math> With <math> u </math> satisfying <math> u (0) = u (1) = 0</math>

Substituting this form into (1) gives
<center>
<math> (b-a)\partial_xh(x)+\partial_x(h(x)\partial_xu) = -\kappa^{2}\left((b-a)x+a+u\right)</math>
</center>
Or, on rearranging
<center>
<math> \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>
</center>
Now consider the homogenous Sturm-Liouville problem for u
<center>
<math> \partial_x(h(x)\partial_xu)+\lambda u = 0\quad u(0)=u(1)=0 \quad (3)</math>
</center>
By Sturm-Liouville theory this has an infinite set of eigenvalues <math> \lambda_k </math> with corresponding eigenfunctions <math> u_k </math>. Also since <math> u_k(0)=u_k(1)=0\quad \forall k </math> Each <math> u_k</math> can be expanded as a fourier series in terms of sine functions.
<center>
<math> u_k = \sum_{n=1}^{\infty} a_{n,k}\sin(n\pi x)</math>
</center>
Transforming (3) into the equivalent variational problem gives
<center>
<math> J[u] = \int_{0}^{1}\,hu'^{2}-\lambda u^{2} \, dx \quad (4) </math>
</center>
Substituting the fourier expansion into (4) implies J must be stationary at <math> \frac{\partial J}{\partial a_{n}}=0 \quad \forall n </math>
<center>
<math> \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>
</center>
By defining a vector <math> \textbf{a} = \left(a_{n}\right)</math> and a matrix <math>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> M\textbf{a} = \lambda\textbf{a} </math> which returns the eigenvalues and eigenvectors of equation (3), with eigenvectors <math> \textbf{a} </math> representing coefficient vectors of the fourier expansions of eigenfunctions.

If we now construct <math> u = \sum_{k=1}^{\infty} b_k u_k </math> and substitute this into equation (2) we get
<center>
<math> \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>
</center>
And defining the RHS of equation (5) as <math> f(x) </math>, a known function, we can retrieve the coefficients <math> b_k </math> by integrating against <math> u_k </math>
<center>
<math> b_k = \frac{\int_{0}^{1}\,f u_k\,dx}{(\kappa^{2}-\lambda_k) \int_{0}^{1}\, u_{k}^{2}\,dx} </math>
</center>

Also to find the coefficients <math> c_n </math> of the fourier expansion of u are just <math> \sum_{k=1}^{\infty}a_{n,k}b_{k} </math> with <math> a_{n.k} </math> being the <math>n</math>th coefficient of the <math>k</math>th eigenfunction of the Sturm-Liouville problem.

So <math> w=(b-a)x+a+\sum_{n=1}^{\infty}c_{n} \sin(n\pi x) </math> with <math> w(0)=a \quad w(1)=b </math> and, given <math> a </math> and <math> b </math> explicitly differentiating <math> w </math> gives <math> w'(0) </math> and <math> w'(1) </math>.

The aim here is to construct a matrix <math> S </math> such that, given <math> a </math> and <math> b </math>
<center>
<math> S \begin{pmatrix} w(0) \\ w(1) \end{pmatrix}=\begin{pmatrix} w'(0) \\ w'(1) \end{pmatrix}</math>
</center>

Taking <math> a=1,\,b=0 </math> shows that the first column of <math> S </math> must be <math> \begin{pmatrix} w'(0) \\ w'(1) \end{pmatrix} \bigg|_{(a=1,b=0)} </math> and likewise taking <math> a=0,\,b=1 </math> shows the second column must be <math> \begin{pmatrix} w'(0) \\ w'(1) \end{pmatrix} \bigg|_{(a=0,b=1)} </math>. So <math>S</math> is given by
<center>
<math> \begin{pmatrix} \begin{matrix} w'(0) \\ w'(1) \end{matrix} \bigg|_{(a=1,b=0)} \begin{matrix} w'(0) \\ w'(1) \end{matrix} \bigg|_{(a=0,b=1)} \end{pmatrix} </math>
</center>

Now considering the case were there is an wave of the form <math> e^{ikx} </math> which, creates reflected and transmitted waves from the variable depth area of the form <math> Re^{ikx} </math> and <math> Te^{ikx} </math> respectively where the magnitudes of R and T are unknown. We can calculate that the boundary conditions for <math> w </math> must be
<center>
<math> w(0) = 1+R, \quad w(1) = Te^{ik}, \quad w'(0) = ik(1-R), \quad w'(1) = ikTe^{ik} </math>
</center>

Now, knowing <math> S </math> these boundary conditions can be solved for <math> R</math> and <math>T</math>, which in turn gives actual numerical boundary conditions ) to the original problem. Taking a linear combination of the solutions already calculated (<math> a=1,\,b=0 </math> and <math> a=0,\,b=1 </math>) will provide the solution for these new boundary conditions.

Variable Depth Shallow Water Wave Equation

2009-01-08T01:08:12Z

Slawrence:

Variable Depth Shallow Water Wave Equation

2009-01-07T22:49:40Z

Slawrence: /* Solution using Separation of Variables */

Variable Depth Shallow Water Wave Equation

2008-12-16T01:57:45Z

Slawrence:

The wave equation can be written as
<center>
<math>\partial_t^2 w- \partial_x \left( c(x)^2 \partial_x w \right) \quad (1)</math>
</center>

Taking a seperable solution <math>\ w (x,t) = \Tau (t) \hat{w} (x)</math> gives the eigenvalue problem
<center>
<math>\partial_x \left( c(x)^2 \partial_x\hat{w} \right) = \lambda\hat{w} \quad (2)</math>
</center>

Given boundary conditions <math>\hat{w} (0) = a</math> and <math>\hat{w} (1) = b</math> we can take
<center>
<math>\hat{w} = (b-a)x + a + u \quad (3)</math>
</center>

With <math> u = \sum_{n=1}^{N} a_n \sin (n \pi x)</math> a series solution satisfying <math> u (0) = u (1) = 0</math>

We wish to solve for <math> a_n </math> given <math> a </math> and <math> b </math>. Equation (2) can be transformed into the Sturm-Liouville problem of minimizing the functional
<center>
<math> J [u] = \int_{0}^{1} c^2 \left(\frac{d \hat{w}}{d x}\right)^2 + \lambda \hat{w}^2 \,dx</math>
</center>
which is minimal exactly when <math>\partial_{a_n} J = 0 \quad \forall a_n </math>

Substituting in <math>\hat{w} = (b-a)x + a + \sum_{n=1}^{N} a_n \sin (n \pi x)</math> and <math>\frac {d \hat{w}}{d x} = (b-a) + \sum_{n=1}^{N} a_n n \pi \cos (n \pi x)</math>

gives

<center>
<math> \partial_{a_n} J = \int_{0}^{1} c^2 \left( n \pi \cos (n \pi x) \left((b-a) + \sum_{m=1}^{N} a_m m \pi \cos (m \pi x) \right) \right) \, dx + \int_{0}^{1} \lambda \left( \sin (n \pi x) \left((b-a)x + a + \sum_{m=1}^{N} a_m \sin (m \pi x) \right) \right) \, dx = 0 \quad (4) </math>
</center>

Since <math> \sin(n \pi x) </math> and <math> \sin(m \pi x) </math> are orthogonal the second integral in (4) can be calculated.

<math> \int_{0}^{1} \lambda \left( \sin (n \pi x) \left((b-a)x + a + \sum_{m=1}^{N} a_m \sin (m \pi x) \right) \right) \, dx =\lambda \int_{0}^{1} \left( (b-a)x+a \right) \sin (n \pi x) \, dx + \frac{\lambda a_n}{2}</math>

and

<math> \int_{0}^{1} \lambda \left( \sin (n \pi x) \left((b-a)x + a + \sum_{m=1}^{N} a_m \sin (m \pi x) \right) \right) \, dx </math>
<math>= \lambda \left((b-a) \int_{0}^{1} x \sin(n \pi x) \, dx + a \int_{0}^{1} \sin(n \pi x) \, dx \right) </math>
<math>= \frac{\lambda (b-a)}{(n \pi)^2} \Big[ \sin(n \pi x) -n \pi x \cos( n \pi x) \Big]_{0}^{1} - \frac{a \lambda}{n \pi} \Big[\cos(n \pi x) \Big]_{0}^{1} </math>
<math> = \frac{\lambda (b-a)}{(n \pi)^2} \left(n \pi (-1)^{n+1} \right) + \frac{a \lambda}{n \pi} \left(1 - (-1)^n \right) </math>
<math> = \frac{\lambda}{n \pi} \left( b (-1)^{n+1}+a \right) </math>

So equation (4) can be written as
<math>\int_{0}^{1} c^2 \left( n \pi \cos (n \pi x) \left((b-a) + \sum_{m=1}^{N} a_m m \pi \cos (m \pi x) \right) \right) \, dx + \frac{\lambda}{n \pi} \left( b (-1)^{n+1}+a \right) + \frac{\lambda a_n}{2}= 0 </math>

The remaining integral can be split into two parts and with the sum taken outside the integral we obtain
<math>\int_{0}^{1} c^2 (b-a) n \pi \cos (n \pi x) + \sum_{m=1}^{N} \left( \int_{0}^{1} c^2 n m \pi^2 \cos(n \pi x) \cos (m \pi x) \, dx \right) a_m + \frac{\lambda}{n \pi} \left( b (-1)^{n+1}+a \right) + \frac{\lambda a_n}{2}= 0 </math>

or, on rearranging
<math> \sum_{m=1}^{N} \left( \int_{0}^{1} c^2 n m \pi^2 \cos(n \pi x) \cos (m \pi x) \, dx \right) a_m+ \frac{\lambda a_n}{2}= -\int_{0}^{1} c^2 (b-a) n \pi \cos (n \pi x) - \frac{\lambda}{n \pi} \left( b (-1)^{n+1}+a \right) </math>

Now if we form the vector <math> \mathbf{a} = \begin{pmatrix} a_1 \\ \vdots \\ a_N \end{pmatrix} </math> and recalling that we have the above expression for all n, we can write the above as a matrix multiplication of <math> \mathbf{a} </math>
<center>
<math> \mathrm{M} \mathbf{a} = \mathbf{f} </math>
</center>

with
<math> \mathrm{M}_{(n,m)} = \int_{0}^{1} c^2 n m \pi^2 \cos(n \pi x) \cos (m \pi x) \, dx + \delta_{nm} \frac{\lambda}{2} \quad \delta_{nm} = \left\{ \begin{matrix} 1 & \mathrm{if} \quad n = m \\ 0 & \mathrm{if} \quad n \neq m \end{matrix} \right. </math>

and
<math> \mathbf{f} = \begin{pmatrix} f_1 \\ \vdots \\ f_N \end{pmatrix} </math> with <math> f_n = -\int_{0}^{1} c^2 (b-a) n \pi \cos (n \pi x) + \frac{\lambda}{n \pi} \left( b (-1)^n-a \right) </math>

So, given a suitable <math> c(x) </math> and boundary conditions <math>\hat{w} (0) = a</math> and <math>\hat{w} (1) = b</math> we have a system of linear equations that can be solved to give the coefficients <math> a_n </math> which in turn define the function <math>\hat{w} </math>.

Variable Depth Shallow Water Wave Equation

2008-12-16T01:42:04Z

Slawrence:

The wave equation can be written as
<center>
<math>\partial_t^2 w- \partial_x \left( c(x)^2 \partial_x w \right) \quad (1)</math>
</center>

Taking a seperable solution <math>\ w (x,t) = \Tau (t) \hat{w} (x)</math> gives the eigenvalue problem
<center>
<math>\partial_x \left( c(x)^2 \partial_x\hat{w} \right) = \lambda\hat{w} \quad (2)</math>
</center>

Given boundary conditions <math>\hat{w} (0) = a</math> and <math>\hat{w} (1) = b</math> we can take
<center>
<math>\hat{w} = (b-a)x + a + u \quad (3)</math>
</center>

With <math> u = \sum_{n=1}^{N} a_n \sin (n \pi x)</math> a series solution satisfying <math> u (0) = u (1) = 0</math>

We wish to solve for <math> a_n </math> given <math> a </math> and <math> b </math>. Equation (2) can be transformed into the Sturm-Liouville problem of minimizing the functional
<center>
<math> J [u] = \int_{0}^{1} c^2 \left(\frac{d \hat{w}}{d x}\right)^2 + \lambda \hat{w}^2 \,dx</math>
</center>
which is minimal exactly when <math>\partial_{a_n} J = 0 \quad \forall a_n </math>

Substituting in <math>\hat{w} = (b-a)x + a + \sum_{n=1}^{N} a_n \sin (n \pi x)</math> and <math>\frac {d \hat{w}}{d x} = (b-a) + \sum_{n=1}^{N} a_n n \pi \cos (n \pi x)</math>

gives

<center>
<math> \partial_{a_n} J = \int_{0}^{1} c^2 \left( n \pi \cos (n \pi x) \left((b-a) + \sum_{m=1}^{N} a_m m \pi \cos (m \pi x) \right) \right) \, dx + \int_{0}^{1} \lambda \left( \sin (n \pi x) \left((b-a)x + a + \sum_{m=1}^{N} a_m \sin (m \pi x) \right) \right) \, dx = 0 \quad (4) </math>
</center>

Since <math> \sin(n \pi x) </math> and <math> \sin(m \pi x) </math> are orthogonal the second integral in (4) can be calculated.

<math> \int_{0}^{1} \lambda \left( \sin (n \pi x) \left((b-a)x + a + \sum_{m=1}^{N} a_m \sin (m \pi x) \right) \right) \, dx =\lambda \int_{0}^{1} \left( (b-a)x+a \right) \sin (n \pi x) \, dx + \frac{\lambda a_n}{2}</math>

and

<math> \int_{0}^{1} \lambda \left( \sin (n \pi x) \left((b-a)x + a + \sum_{m=1}^{N} a_m \sin (m \pi x) \right) \right) \, dx </math>
<math>= \lambda \left((b-a) \int_{0}^{1} x \sin(n \pi x) \, dx + a \int_{0}^{1} \sin(n \pi x) \, dx \right) </math>
<math>= \frac{\lambda (b-a)}{(n \pi)^2} \Big[ \sin(n \pi x) -n \pi x \cos( n \pi x) \Big]_{0}^{1} - \frac{a \lambda}{n \pi} \Big[\cos(n \pi x) \Big]_{0}^{1} </math>
<math> = \frac{\lambda (b-a)}{(n \pi)^2} \left(n \pi (-1)^{n+1} \right) + \frac{a \lambda}{n \pi} \left(1 - (-1)^n \right) </math>
<math> = \frac{\lambda}{n \pi} \left( b (-1)^{n+1}+a \right) </math>

So equation (4) can be written as
<math>\int_{0}^{1} c^2 \left( n \pi \cos (n \pi x) \left((b-a) + \sum_{m=1}^{N} a_m m \pi \cos (m \pi x) \right) \right) \, dx + \frac{\lambda}{n \pi} \left( b (-1)^{n+1}+a \right) + \frac{\lambda a_n}{2}= 0 </math>

The remaining integral can be split into two parts and with the sum taken outside the integral we obtain
<math>\int_{0}^{1} c^2 (b-a) n \pi \cos (n \pi x) + \sum_{m=1}^{N} \left( \int_{0}^{1} c^2 n m \pi^2 \cos(n \pi x) \cos (m \pi x) \, dx \right) a_m + \frac{\lambda}{n \pi} \left( b (-1)^{n+1}+a \right) + \frac{\lambda a_n}{2}= 0 </math>

or, on rearranging
<math> \sum_{m=1}^{N} \left( \int_{0}^{1} c^2 n m \pi^2 \cos(n \pi x) \cos (m \pi x) \, dx \right) a_m+ \frac{\lambda a_n}{2}= -\int_{0}^{1} c^2 (b-a) n \pi \cos (n \pi x) - \frac{\lambda}{n \pi} \left( b (-1)^{n+1}+a \right) </math>

Now if we form the vector <math> \mathbf{a} = \begin{pmatrix} a_1 \\ \vdots \\ a_N \end{pmatrix} </math> and recalling that we have the above expression for all n, we can write the above as a matrix multiplication of <math> \mathbf{a} </math>
<center>
<math> \mathrm{M} \mathbf{a} = \mathbf{f} </math>
</center>

with
<math> \mathrm{M}_{(n,m)} = \int_{0}^{1} c^2 n m \pi^2 \cos(n \pi x) \cos (m \pi x) \, dx + \delta_{nm} \frac{\lambda}{2} \quad \delta_{nm} = \left\{ \begin{matrix} 1 & \mathrm{if} \quad n = m \\ 0 & \mathrm{if} \quad n \neq m \end{matrix} \right. </math>

and
<math> \mathbf{f} = \begin{pmatrix} f_1 \\ \vdots \\ f_N \end{pmatrix} </math> with <math> f_n = -\int_{0}^{1} c^2 (b-a) n \pi \cos (n \pi x) - \frac{\lambda}{n \pi} \left( b (-1)^{n+1}+a \right) </math>

So, given a suitable <math> c(x) </math> and boundary conditions <math>\hat{w} (0) = a</math> and <math>\hat{w} (1) = b</math> we have a system of linear equations that can be solved to give the coefficients <math> a_n </math> which in turn define the function <math>\hat{w} </math>.