View source for Waves on a Variable Beam

== Introduction ==
We consider here the equations for a non-uniform free beam. We begin by presenting the theory for a uniform beam.<br\>

An example of the motion for a non-uniform beam is demonstrated below:
[[Image:variable_beam_ii.gif|right|Vibration of variable beam]]

== Equations for a beam ==
{{equations for a beam}}

== Solution for a uniform beam in [[Eigenfunctions for a Uniform Free Beam|eigenfunctions]] with no pressure ==

{{solution for a uniform beam in eigenfunctions}}

== Variational Techniques ==
As an alternative to solving the eigenvalue problem, we can equivalently minimise the [http://en.wikipedia.org/wiki/Rayleigh_quotient Rayleigh Quotient] ([[Linton and McIver 2001]]):
<center><math>  R[\zeta]=\frac{\varepsilon[\zeta]}{H[\zeta]}   </math></center>
where <math> \varepsilon[\zeta] \,\!</math> is known as the energy functional, and <math> H[\zeta] \,\!</math> is some constraint.
<br /><br />
In the case of a uniform beam with zero  transverse load, our energy functional is ([[Lanczos 1949]]):
<center><math>  \varepsilon[\zeta]=\frac{\beta}{2}\int_{-L}^{L}\left(\frac{\partial^2 v}{\partial x^2} \right)^2 \mathrm{d}x    </math></center>
If we choose
<center><math>  H[\zeta]=\int_{-L}^{L}\gamma \zeta^2 \mathrm{d}x=1   </math></center>
Then minimizing the Rayleigh Quotient can be expressed as a variational problem <br />
<center>
min <math>  R[\zeta]=\frac{\beta}{2}\int_{-L}^{L}\left(\frac{\partial^2 \zeta}{\partial x^2} \right)^2 \mathrm{d}x   </math>
subject to <math> \int_{-L}^{L}\gamma \zeta^2 \mathrm{d}x=1   </math></center>
which will in turn, also solve our eigenvalue problem.

== Rayleigh-Ritz method ==
We can essentially replace the variational problem of finding a <math> y(x) \,\!</math> that extremises <math> J[y] \,\!</math> to finding a set of constants <math> a_1, a_2, ..., a_N \,\!</math> which extremise <math> J[a_1, a_2, ..., a_N] \,\!</math>, through using the [http://en.wikipedia.org/wiki/Rayleigh-Ritz_method Rayleigh-Ritz method].<br /><br />
We solve
<center><math> \frac{\partial}{\partial a_k}J[a_1, a_2, ..., a_N]=0 \,\!</math></center>
for all <math> k=1,2,...,N \,\!</math>.<br /><br />
In our example of non-uniform beams, the assumption is made that we are able to approximate the eigenfunctions for a non-uniform beam as a linear combination of the eigenfunctions for a linear beam:
<center><math> \widehat{X}_n=\sum_{i=1}^{N}a_{i}X_{i}(x) \,\!</math></center>

== Solution for a Non-uniform free beam in eigenfunctions ==
In the case of a non-uniform free beam, the Euler-Bernoulli beam equation becomes:
<center><math> \frac{\partial^2}{\partial x^2}\left( \beta(x) \frac{\partial^{2}\zeta}{\partial x^{2}}\right) +\gamma(x)\frac{\partial^{2}\zeta}{\partial t^{2}}=0 \,\!</math></center>
Using separation of variables, where <math> \zeta(x,t)=\widehat{X}(x)\widehat{T}(t) \,\!</math>, we obtain a corresponding eigenfunction problem (with eigenvalues denoted by <math> \mu_n=\widehat{k}_n^4 \,\!</math>):
<center><math> \frac{\partial^2}{\partial x^2}\left( \beta(x) \frac{\partial^{2} \widehat{X}_n}{\partial x^{2}}\right)=\gamma(x)\mu_n \widehat{X}_n \,\!</math></center>
which leads us to the following variational expression:
<center>
min <math>  J[\widehat{X}_n]=\frac{1}{2}\int_{-L}^{L} \bigg\{\beta(x)(\widehat{X}_n^{''})^2 -\mu_n \gamma(x) \widehat{X}_n^2 \bigg\}\mathrm{d}x    </math>
</center>
We can approximate this using Rayleigh-Ritz and obtain:
<center><math> J[\vec{a}]=\frac{1}{2}\int_{-L}^{L} \bigg\{\beta(x)[\sum_{i=1}^{N}a_{i}{X}_{i}^{''}]^2 -\mu_i \gamma(x) [\sum_{i=1}^{N}a_{i}{X}_{i}]^2 \bigg\}\mathrm{d}x  \,\!</math></center>
Then extremising <math> J[\vec{a}] \,\!</math> as follows
<center><math> \frac{\partial}{\partial a_k}J[\vec{a}]=0 \,\!</math></center>
for all <math> k=1,2,...,N \,\!</math>, allows us to obtain:
<center><math> \frac{\partial}{\partial a_k}J[\vec{a}]=\sum_{i=1}^{N} \left\{ \int_{-L}^{L}\beta(x) {X}_{i}^{''} {X}_{k}^{''}\mathrm{d}x -\mu_i \int_{-L}^{L} \gamma(x){X}_{i}{X}_{k}\mathrm{d}x\right\}a_i=0 \,\!</math></center>
for all <math> k=1,2,...N \,\!</math> (here <math> \mu_i \,\!</math> denotes both the [http://en.wikipedia.org/wiki/Lagrange_multipliers Lagrange multiplier] and eigenvalues of a non-uniform beam)

== Converting to Matrix Form ==
If we define
<center>
:<math> K_{ik}=\int_{-L}^{L}D(x) {X}_{i}^{''} {X}_{k}^{''}\mathrm{d}x \,\!</math>
:<math> M_{ik}=\int_{-L}^{L} m(x){X}_{i}{X}_{k}\mathrm{d}x \,\!</math>
</center>
then for all <math> k=1,2,...N \,\!</math>:
<center><math> \frac{\partial}{\partial a_k}J[\mathbf{a}]=\sum_{i=1}^{N}(K_{ik}-\mu_i M_{ik})a_i=0 \,\!</math></center>
Let us create the matrix <math> K \,\!</math>, with elements <math> K_{ik} \,\!</math>, the matrix <math> M \,\!</math>, with elements <math> M_{ik} \,\!</math>, the sparse matrix <math> \Lambda \,\!</math> with <math> \mu_i \,\!</math> terms on the diagonal, and the vector <math> \mathbf{a} \,\!</math>.  We can consequently express the above sum in the following way:
<center><math> (K-\Lambda M)\mathbf{a} = 0 \,\!</math></center>
So we can solve for <math> \mathbf{a} \,\!</math> in any of the problems below to obtain coefficients <math> a_i \,\!</math>:
<center>
:<math> (K-\Lambda M)\mathbf{a} = 0  \,\!</math>
:<math> K\mathbf{a}=\Lambda M\mathbf{a}  \,\!</math>
:<math> (M^{-1}K)\mathbf{a}=\Lambda \mathbf{a} \,\!</math>
</center>

== Quadrature ==
Rather than tediously evaluating an integral for each element of the matrices <math> K \,\!</math> and <math> M \,\!</math>, we can break up the integral into subintervals with weights <math> w_h \,\!</math> using [http://en.wikipedia.org/wiki/Composite_simpson%27s_rule#Composite_Simpson.27s_rule Composite Simpson's Rule]:
<center>
:<math> \int_a^b f(x) dx \approx \frac{h}{3} \left[ f(x_0)+2 \sum_{j=1}^{n/2-1}f(x_{2j})+4\sum_{j=1}^{n/2}f(x_{2j-1})+f(x_n)\right] \,\!</math>
:<math> \approx \frac{h}{3}\left[ f(x_0)+4f(x_1)+2f(x_2)+4f(x_3)+2f(x_4)+...+4f(x_{n-1})+f(x_n)\right] \,\!</math>
</center>
So we can express <math> K_{ik} \,\!</math> as follows:
<center>
:<math> K_{ik}=\int_{-L}^{L}D(x) {X}_{i}^{''} {X}_{k}^{''}\mathrm{d}x \,\!</math>
:<math> \approx \sum D(x_h){X}_{i}^{''}(x_h) {X}_{k}^{''}(x_h)w_h \,\!</math>
:<math> \approx \vec{X}_{i}^{''} H \vec{X}_{k}^{''T}  \,\!</math>
</center>
where
<center><math> \vec{X}_i^{''}=[X_i^{''}(x_1),X_i^{''}(x_2),...,X_i^{''}(x_h)] \,\!</math>
<math> H =
\begin{bmatrix}
  D(x_1)w_1 & 0 & ... &0 \\
  0 & D(x_2)w_2 & ... &0 \\
  \vdots &  & \ddots&\vdots \\
  0 & ... & ...&D(x_h)w_h
\end{bmatrix}
\,\!</math>
</center>
and weights <math> w_h \,\!</math> are defined as: <math> w_1=h/3, w_2=4h/3, ... , w_h=h/3 \,\!</math>.<br /><br />
We can extend this concept to the the full matrix <math> K \,\!</math> if we form:
<center><math> K=X^{''}_{mat}H X^{''T}_{mat}  \,\!</math></center>
where
<center>
<math> X^{''}_{mat} =
\begin{bmatrix}
{X}_{1}^{''}(x_1) & {X}_{1}^{''}(x_2) & ... & {X}^{''}_{1}(x_h) \\
{X}_{2}^{''}(x_1) & {X}_{2}^{''}(x_2) & ... & {X}^{''}_{2}(x_h) \\
\vdots & \vdots & & \vdots \\
{X}_{N}^{''}(x_1) & {X}_{N}^{''}(x_2) & ... & {X}^{''}_{N}(x_h)
\end{bmatrix}
\,\!</math>
</center>
For <math> M \,\!</math> we use an identical approach:
<center><math> M=X_{mat}J X^{T}_{mat} \,\!</math></center>
where
<center>
<math> J =
\begin{bmatrix}
m(x_1)w_1 & 0 & ... &0 \\
0 & m(x_2)w_2 & ... &0 \\
\vdots &  & \ddots&\vdots \\
0 & ... & ...&m(x_h)w_h
\end{bmatrix}
\,\!</math>
</center>
<center>
<math> X_{mat} =
\begin{bmatrix}
{X}_{1}(x_1) & {X}_{1}(x_2) & ... & {X}_{1}(x_h) \\
{X}_{2}(x_1) & {X}_{2}(x_2) & ... & {X}_{2}(x_h) \\
\vdots & \vdots & & \vdots \\
{X}_{N}(x_1) & {X}_{N}(x_2) & ... & {X}_{N}(x_h)
\end{bmatrix}
\,\!</math>
</center>

== Evaluation ==
Using the MATLAB program eig.m, we can solve the eigenvalue problem
<center><math> K\vec{a}=\Lambda M\vec{a} \,\!</math></center>
using the either of the commands <br /><br />
[Ai_values,nonuniform_eigvals] = eig(M^(-1)*K) <br />
[Ai_values,nonuniform_eigvals] = eig(K,M)
<br /><br />
Where the diagonal of nonuniform_eigvals denotes the eigenvalues for the nonuniform beam (<math> \mu_n \,\!</math>), and each column of the matrix Ai_values represents the coefficients <math> a_1, a_2,... \,\!</math> corresponding to <math> \widehat{X}_n \,\!</math>:
<center><math> \widehat{X}_n=\sum_{i=1}^{N}a_{i}X_{i}(x) \,\!</math></center>
We can obtain the matrix <math> \widehat{X}_{mat} \,\!</math> by taking Xhat_mat= Ai_values<sup>T</sup>X_mat.

== Non-uniform beam revisited ==
We have already solved the eigenfunction problem. We now turn our attention to the time component arising from separation of variables:
<center><math> \widehat{T}_n^{''}+\mu_n\widehat{T}_n=0 \,\!</math></center>
which has solutions of the form
<center><math> \widehat{T}_n(t)=A_n \cos(\widehat{k}_n^2 t)+B_n \sin(\widehat{k}_n^2 t) \,\!</math></center>
Introducing the initial conditions
<center>
:<math> v(x,0)=f(x) \,\!</math>
:<math>  \frac{\partial v(x,0)}{\partial t}=g(x)\,\!</math>
</center>
Then using the the first of these initial conditions we obtain:
<center><math> \sum_{n=0}^{\infty}A_n\widehat{X}_n(x)=f(x) \,\!</math></center>
multiplying by <math> m(x)\widehat{X}_n(x) \,\!</math> and integrating allows us to obtain:
<center><math> A_n=\frac{\int_{-L}^{L}f(x)m(x)\widehat{X}_n(x)\mathrm{d}x}{\int_{-L}^{L}m(x)\widehat{X}_n(x)\widehat{X}_n(x)\mathrm{d}x} \,\!</math></center>
Then using the second initial condition, and applying the same technique yields
<center><math> B_n=\frac{1}{\widehat{k}_n^2}\frac{\int_{-L}^{L}g(x)m(x)\widehat{X}_n(x)\mathrm{d}x}{\int_{-L}^{L}m(x)\widehat{X}_n(x)\widehat{X}_n(x)\mathrm{d}x} \,\!</math></center>
Consequently,
<center><math> v(x,t)=A_0\widehat{X}_0(x)+A_1\widehat{X}_1(x)+\sum_{n=2}^{\infty} A_n\widehat{X}_n(x)\cos(\widehat{k}_n^2 t)+\sum_{n=2}^{\infty} B_n\widehat{X}_n(x) \sin(\widehat{k}_n^2 t)  \,\!</math></center>
where <math> \widehat{X}_0(x) \,\!</math> and <math> \widehat{X}_2(x) \,\!</math> are the two rigid modes of the non-uniform beam.  Note that <math> B_0 \,\!</math> and <math> B_1 \,\!</math> do not exist.

== Matlab Code ==

A program to solve for a variable beam can be found here {{variable beam}}

=== Additional code ===

This program requires
* {{eigenvalues beam}}
* {{eigenvectors beam}}
* {{double derivative eigenvectors beam}}

[[Category:Simple Linear Waves]]

==References==
[[Lanczos 1949]],
[[Linton and McIver 2001]],
[[Rao 1986]]

[[Category:Complete Pages]]