### 2.3 HW 3

2.3.1 Problem 1
2.3.2 Problem 2
2.3.3 Problem 3

1.
This HW in one PDF (letter size)
2.
This HW in one PDF (legal size)

#### 2.3.1 Problem 1

2.3.1.1 Results
2.3.1.2 Source code

The ODE is$y^{\prime \prime }+2y^{\prime }+\lambda ^{2}y=0$ With $$y\left ( 0\right ) =y\left ( 1\right ) =0$$. The ﬁrst step is to ﬁnd the state space representation. Let $$x_{1}=y,x_{2}=y^{\prime }.$$ Taking derivatives gives\begin{align*} \dot{x}_{1} & =x_{2}\\ \dot{x}_{2} & =-2x_{2}-\lambda ^{2}x_{1} \end{align*}

The above is used with bvp4c as shown in the source code. To use eig, the problem is converted to the form $$Ay=\alpha By$$ and then Matlab $$eig(A,B)$$ is used to ﬁnd the eigenvalues. Using second order centered diﬀerence gives\begin{align*} \left . \frac{dy}{dx}\right \vert _{i} & =\frac{y_{i+1}-y_{i-1}}{2h}\\ \left . \frac{d^{2}y}{dx^{2}}\right \vert _{i} & =\frac{y_{i+1}-2y_{i}+y_{i-1}}{h^{2}} \end{align*}

Therefore, the approximation to the diﬀerential equation at grid $$i$$ (on the internal nodes as shown in the above diagram) is$\left . y^{\prime \prime }+2y^{\prime }+\lambda ^{2}y\right \vert _{i}\approx \frac{y_{i+1}-2y_{i}+y_{i-1}}{h^{2}}+2\frac{y_{i+1}-y_{i-1}}{2h}+\lambda ^{2}y_{i}$ Hence\begin{align*} \frac{y_{i+1}-2y_{i}+y_{i-1}}{h^{2}}+2\frac{y_{i+1}-y_{i-1}}{2h}+\lambda ^{2}y_{i} & =0\\ y_{i+1}-2y_{i}+y_{i-1}+h\left ( y_{i+1}-y_{i-1}\right ) +h^{2}\lambda ^{2}y_{i} & =0\\ y_{i-1}\left ( 1-h\right ) +y_{i}\left ( h^{2}\lambda ^{2}-2\right ) +y_{i+1}\left ( 1+h\right ) & =0 \end{align*}

At node $$i=1$$,$y_{0}\left ( 1-h\right ) +y_{1}\left ( h^{2}\lambda ^{2}-2\right ) +y_{2}\left ( 1+h\right ) =0$ Moving the known quantities and any quantity with $$\lambda$$ to the right side$-2y_{1}+y_{2}\left ( 1+h\right ) =y_{0}\left ( h-1\right ) -y_{1}\left ( h^{2}\lambda ^{2}\right )$ At node $$i=2$$\begin{align*} y_{1}\left ( 1-h\right ) +y_{2}\left ( h^{2}\lambda ^{2}-2\right ) +y_{3}\left ( 1+h\right ) & =0\\ y_{1}\left ( 1-h\right ) -2y_{2}+y_{3}\left ( 1+h\right ) & =-y_{2}\left ( h^{2}\lambda ^{2}\right ) \end{align*}

And so on. At the last node, $$i=N$$\begin{align*} y_{N-1}\left ( 1-h\right ) +y_{N}\left ( h^{2}\lambda ^{2}-2\right ) +y_{N+1}\left ( 1+h\right ) & =0\\ y_{N-1}\left ( 1-h\right ) -2y_{N} & =-y_{N}\left ( h^{2}\lambda ^{2}\right ) -y_{N+1}\left ( 1+h\right ) \end{align*}

At $$i=N-1$$\begin{align*} y_{N-2}\left ( 1-h\right ) +y_{N-1}\left ( h^{2}\lambda ^{2}-2\right ) +y_{N}\left ( 1+h\right ) & =0\\ y_{N-2}\left ( 1-h\right ) -2y_{N-1}+y_{N}\left ( 1+h\right ) & =-y_{N-1}\left ( h^{2}\lambda ^{2}\right ) \end{align*}

Hence the structure is$\begin{bmatrix} -2 & 1+h & 0 & 0 & 0 & \cdots & 0\\ 1-h & -2 & 1+h & 0 & 0 & \cdots & \vdots \\ 0 & 1-h & -2 & 1+h & 0 & \cdots & \vdots \\ 0 & 0 & \cdots & \ddots & \cdots & \cdots & \vdots \\ \vdots & \cdots & \cdots & \cdots & \ddots & \cdots & 0\\ \vdots & \cdots & \cdots & \cdots & 1-h & -2 & 1+h\\ 0 & \cdots & \cdots & \cdots & 0 & 1-h & -2 \end{bmatrix}\begin{bmatrix} y_{1}\\ y_{2}\\ y_{3}\\ \vdots \\ y_{N-2}\\ y_{N-1}\\ y_{N}\end{bmatrix} =\begin{bmatrix} y_{0}\left ( h-1\right ) \\ 0\\ 0\\ \vdots \\ 0\\ 0\\ -y_{N+1}\left ( 1+h\right ) \end{bmatrix} -h^{2}\lambda ^{2}\begin{bmatrix} 1 & 0 & 0 & 0 & 0 & \cdots & 0\\ 0 & 1 & 0 & 0 & 0 & \cdots & \vdots \\ 0 & 0 & 1 & 0 & 0 & \cdots & \vdots \\ 0 & 0 & 0 & 1 & 0 & \cdots & \vdots \\ \vdots & \cdots & \cdots & \cdots & \ddots & \cdots & \vdots \\ \vdots & \vdots & \vdots & \vdots & \vdots & \ddots & 0\\ 0 & 0 & 0 & 0 & \cdots & 0 & 1 \end{bmatrix}\begin{bmatrix} y_{1}\\ y_{2}\\ y_{3}\\ \vdots \\ y_{N-2}\\ y_{N-1}\\ y_{N}\end{bmatrix}$ Since $$y_{0}=y_{N+1}=0$$ the above reduces to\begin{align*} \begin{bmatrix} -2 & 1+h & 0 & 0 & 0 & \cdots & 0\\ 1-h & -2 & 1+h & 0 & 0 & \cdots & \vdots \\ 0 & 1-h & -2 & 1+h & 0 & \cdots & \vdots \\ 0 & 0 & \cdots & \ddots & \cdots & \cdots & \vdots \\ \vdots & \cdots & \cdots & \cdots & \ddots & \cdots & 0\\ \vdots & \cdots & \cdots & \cdots & 1-h & -2 & 1+h\\ 0 & \cdots & \cdots & \cdots & 0 & 1-h & -2 \end{bmatrix}\begin{bmatrix} y_{1}\\ y_{2}\\ y_{3}\\ \vdots \\ y_{N-2}\\ y_{N-1}\\ y_{N}\end{bmatrix} & =-h^{2}\lambda ^{2}\begin{bmatrix} 1 & 0 & 0 & 0 & 0 & \cdots & 0\\ 0 & 1 & 0 & 0 & 0 & \cdots & \vdots \\ 0 & 0 & 1 & 0 & 0 & \cdots & \vdots \\ 0 & 0 & 0 & 1 & 0 & \cdots & \vdots \\ \vdots & \cdots & \cdots & \cdots & \ddots & \cdots & \vdots \\ \vdots & \vdots & \vdots & \vdots & \vdots & \ddots & 0\\ 0 & 0 & 0 & 0 & \cdots & 0 & 1 \end{bmatrix}\begin{bmatrix} y_{1}\\ y_{2}\\ y_{3}\\ \vdots \\ y_{N-2}\\ y_{N-1}\\ y_{N}\end{bmatrix} \\ Ay & =\alpha By \end{align*}

Where $$\alpha =\lambda ^{2}$$ and $$B=-h^{2}\mathbf{I}$$. The above is now implemented in Matlab and eig is used to ﬁnd $$\alpha$$.

The analytical value of the eigenvalue is given from\begin{align*} \sqrt{\lambda _{n}^{2}-1} & =n\pi \\ \lambda _{n} & =\sqrt{\left ( n\pi \right ) ^{2}+1} \end{align*}

Hence the ﬁrst three eigenvalues are \begin{align*} \lambda _{1} & =\sqrt{\pi ^{2}+1}=3.2969\\ \lambda _{2} & =\sqrt{\left ( 2\pi \right ) ^{2}+1}=6.3623\\ \lambda _{3} & =\sqrt{\left ( 3\pi \right ) ^{2}+1}=9.4777 \end{align*}

And the corresponding analytical mode shapes, using $$A_{n}=1$$ when normalized is\begin{align*} y_{1}\left ( x\right ) & =e^{-x}\sin \left ( \pi x\right ) \\ y_{2}\left ( x\right ) & =e^{-x}\sin \left ( 2\pi x\right ) \end{align*}

These are used to compare the numerical solutions from bvp4c and from eig against. The following plots show the result for the ﬁrst three eigenvalues and eigenfunctions found. The main diﬃculty with using bvp4c for solving the eigenvalue problem is on deciding which guess $$\lambda$$ to use for each mode shape to solve for. The ﬁrst three mode shapes are solved for, and also a plot of the initial mode shape guess passed to bvp4c is plotted. Using large grid size, the solution by eig and bvp4c matched very well as can be seen from the plots below. The eigenvalue produced by bvp4c was little closer to the analytical one than the eigenvalue produced by eig() command.

##### 2.3.1.1 Results

Each mode shape plot is given, showing the eigenvalue produced by each solver and the initial mode shape guess used. There are 3 plots, one for each mode shape. The ﬁrst, second and third. (the problem asked for only the ﬁrst two mode shapes, but the third one was added for veriﬁcation).

1.
First mode shape
 Solver eigenvalue found analytical $$\lambda _{1}=\sqrt{\pi ^{2}+1}=$$ $$3.296\,9$$ bvp4c $$3.2969$$ eig $$3.2960997$$

2.
Second mode shape
 Solver eigenvalue found analytical $$\lambda _{1}=\sqrt{\left ( 2\pi \right ) ^{2}+1}=$$ $$6.3623$$ bvp4c $$6.3622$$ eig $$6.35738$$

3.
Third mode shape
 Solver eigenvalue found analytical $$\lambda _{1}=\sqrt{\left ( 3\pi \right ) ^{2}+1}=$$ $$9.4777$$ bvp4c $$9.4777$$ eig $$9.4623$$

Printout of Matlab console running the program

>>nma_HW3_EMA_471_problem_1
*********************
running mode 1. Eigenvalue, obtained with bvp4c, is 3.2968962.
eigenvalue from eig is 3.2960997
*********************
running mode 2. Eigenvalue, obtained with bvp4c, is 6.3622025.
eigenvalue from eig is 6.3573774
*********************
running mode 3. Eigenvalue, obtained with bvp4c, is 9.4777434.
eigenvalue from eig is 9.4622719

#### 2.3.2 Problem 2

2.3.2.1 Power method
2.3.2.2 Results
2.3.2.3 Source code

The geometry of the problem is as follows

Using the normalized ODE

$z^{4}y^{\prime \prime }+\lambda ^{2}y=0$

With BC \begin{align*} y\left ( \frac{a}{L}\right ) & =y\left ( \frac{3}{3}\right ) =y\left ( 1\right ) =0\\ y\left ( \frac{b}{L}\right ) & =y\left ( \frac{6}{3}\right ) =y\left ( 2\right ) =0 \end{align*}

And $\lambda ^{2}=\frac{Pb^{4}}{EL^{2}I_{0}}$

For domain $$1\leq z\leq 2$$. The analytical solution is $$P_{n}=n^{2}\pi ^{2}\left ( \frac{a}{b}\right ) ^{2}\frac{EI_{0}}{L^{2}}$$ and $$y_{n}\left ( z\right ) =A_{n}z\sin \left ( n\pi \frac{b}{L}\left ( 1-\frac{a}{Lz}\right ) \right )$$. The ﬁrst step is to convert the ODE into state space for use with bvp4c. Let $$x_{1}=y,x_{2}=y^{\prime }$$. Taking derivatives gives\begin{align*} \dot{x}_{1} & =x_{2}\\ \dot{x}_{2} & =-\frac{\lambda ^{2}}{z^{4}}x_{1} \end{align*}

For using eig, the problem needs to discretized ﬁrst. The following shows the grid used

The grid starts at $$z=1$$ and ends at $$z=2$$ since this is the domain of the diﬀerential equation being solved. Therefore $$i=0$$ corresponds to the left boundary conditions which is $$y(z=1)=0$$ and $$i=(N+2)h$$ corresponds to the right boundary conditions which is $$y(z=2)=0$$.

Using second order centered diﬀerence gives$\left . \frac{d^{2}y}{dx^{2}}\right \vert _{i}=\frac{y_{i+1}-2y_{i}+y_{i-1}}{h^{2}}$ Therefore, the approximation to the diﬀerential equation at grid $$i$$ (on the internal nodes as shown in the above diagram) is as follows. Notice we needed to add $$1$$ to the grid spacing since the left boundary starts at $$z=1$$ in this case and not at $$z=0$$ as normally the case in other problems.$\left . z^{4}y^{\prime \prime }+\lambda ^{2}y\right \vert _{i}\approx \left ( 1+ih\right ) ^{4}\frac{y_{i+1}-2y_{i}+y_{i-1}}{h^{2}}+\lambda ^{2}y_{i}$ Hence\begin{align*} \left ( 1+ih\right ) ^{4}\left ( y_{i+1}-2y_{i}+y_{i-1}\right ) +h^{2}\lambda ^{2}y_{i} & =0\\ \left ( 1+ih\right ) ^{4}y_{i+1}-2\left ( 1+ih\right ) ^{4}y_{i}+\left ( 1+ih\right ) ^{4}y_{i-1} & =-h^{2}\lambda ^{2}y_{i} \end{align*}

At node $$i=1$$,$\left ( 1+h\right ) ^{4}y_{2}-2\left ( 1+h\right ) ^{4}y_{1}+\left ( 1+h\right ) ^{4}y_{0}=-h^{2}\lambda ^{2}y_{1}$ Moving the known quantities to the right side$\left ( 1+h\right ) ^{4}y_{2}-2\left ( 1+h\right ) ^{4}y_{1}=-\left ( 1+h\right ) ^{4}y_{0}-h^{2}\lambda ^{2}y_{1}$ At node $$i=2$$$\left ( 1+2h\right ) ^{4}y_{3}-2\left ( 1+2h\right ) ^{4}y_{2}+\left ( 1+2h\right ) ^{4}y_{1}=-h^{2}\lambda ^{2}y_{2}$ And so on. At the last node, $$i=N$$\begin{align*} \left ( 1+Nh\right ) ^{4}y_{N+1}-2\left ( 1+Nh\right ) ^{4}y_{N}+\left ( 1+Nh\right ) ^{4}y_{N-1} & =-h^{2}\lambda ^{2}y_{N}\\ -2\left ( Nh\right ) ^{4}y_{N}+\left ( Nh\right ) ^{4}y_{N-1} & =-\left ( 1+Nh\right ) ^{4}y_{N+1}-h^{2}\lambda ^{2}y_{N} \end{align*}

At $$i=N-1$$$\left ( 1+\left ( N-1\right ) h\right ) ^{4}y_{N}-2\left ( 1+\left ( N-1\right ) h\right ) ^{4}y_{N-1}+\left ( 1+\left ( N-1\right ) h\right ) ^{4}y_{N-2}=-h^{2}\lambda ^{2}y_{\left ( N-1\right ) }$ Hence the structure is$\begin{bmatrix} -2\left ( 1+h\right ) ^{4} & \left ( 1+h\right ) ^{4} & 0 & 0 & 0 & \cdots & 0\\ \left ( 1+2h\right ) ^{4} & -2\left ( 1+2h\right ) ^{4} & \left ( 1+2h\right ) ^{4} & 0 & 0 & \cdots & \vdots \\ 0 & \left ( 1+3h\right ) ^{4} & -2\left ( 1+3h\right ) ^{4} & \left ( 1+3h\right ) ^{4} & 0 & \cdots & \vdots \\ 0 & 0 & \cdots & \ddots & \cdots & \cdots & \vdots \\ \vdots & \cdots & \cdots & \cdots & \ddots & \cdots & 0\\ \vdots & \cdots & \cdots & \cdots & \left ( 1+\left ( N-1\right ) h\right ) ^{4} & -2\left ( 1+\left ( N-1\right ) h\right ) ^{4} & \left ( 1+\left ( N-1\right ) h\right ) ^{4}\\ 0 & \cdots & \cdots & \cdots & 0 & \left ( 1+Nh\right ) ^{4} & -2\left ( 1+Nh\right ) ^{4}\end{bmatrix}\begin{bmatrix} y_{1}\\ y_{2}\\ y_{3}\\ \vdots \\ y_{N-2}\\ y_{N-1}\\ y_{N}\end{bmatrix} =$ $= \begin{bmatrix} -\left ( 1+h\right ) ^{4}y_{0}\\ 0\\ 0\\ \vdots \\ 0\\ 0\\ -\left ( 1+Nh\right ) ^{4}y_{N+1}\end{bmatrix} -h^{2}\lambda ^{2}\begin{bmatrix} 1 & 0 & 0 & 0 & 0 & \cdots & 0\\ 0 & 1 & 0 & 0 & 0 & \cdots & \vdots \\ 0 & 0 & 1 & 0 & 0 & \cdots & \vdots \\ 0 & 0 & 0 & 1 & 0 & \cdots & \vdots \\ \vdots & \cdots & \cdots & \cdots & \ddots & \cdots & \vdots \\ \vdots & \vdots & \vdots & \vdots & \vdots & \ddots & 0\\ 0 & 0 & 0 & 0 & \cdots & 0 & 1 \end{bmatrix}\begin{bmatrix} y_{1}\\ y_{2}\\ y_{3}\\ \vdots \\ y_{N-2}\\ y_{N-1}\\ y_{N}\end{bmatrix}$ Since $$y_{0}=y_{N+1}=0$$ the above reduces to$\begin{bmatrix} -2\left ( 1+h\right ) ^{4} & \left ( 1+h\right ) ^{4} & 0 & 0 & 0 & \cdots & 0\\ \left ( 1+2h\right ) ^{4} & -2\left ( 1+2h\right ) ^{4} & \left ( 1+2h\right ) ^{4} & 0 & 0 & \cdots & \vdots \\ 0 & \left ( 1+3h\right ) ^{4} & -2\left ( 1+3h\right ) ^{4} & \left ( 1+3h\right ) ^{4} & 0 & \cdots & \vdots \\ 0 & 0 & \cdots & \ddots & \cdots & \cdots & \vdots \\ \vdots & \cdots & \cdots & \cdots & \ddots & \cdots & 0\\ \vdots & \cdots & \cdots & \cdots & \left ( 1+\left ( N-1\right ) h\right ) ^{4} & -2\left ( 1+\left ( N-1\right ) h\right ) ^{4} & \left ( 1+\left ( N-1\right ) h\right ) ^{4}\\ 0 & \cdots & \cdots & \cdots & 0 & \left ( 1+Nh\right ) ^{4} & -2\left ( 1+Nh\right ) ^{4}\end{bmatrix}\begin{bmatrix} y_{1}\\ y_{2}\\ y_{3}\\ \vdots \\ y_{N-2}\\ y_{N-1}\\ y_{N}\end{bmatrix}$ \begin{align*} &= -h^{2}\lambda ^{2}\begin{bmatrix} 1 & 0 & 0 & 0 & 0 & \cdots & 0\\ 0 & 1 & 0 & 0 & 0 & \cdots & \vdots \\ 0 & 0 & 1 & 0 & 0 & \cdots & \vdots \\ 0 & 0 & 0 & 1 & 0 & \cdots & \vdots \\ \vdots & \cdots & \cdots & \cdots & \ddots & \cdots & \vdots \\ \vdots & \vdots & \vdots & \vdots & \vdots & \ddots & 0\\ 0 & 0 & 0 & 0 & \cdots & 0 & 1 \end{bmatrix}\begin{bmatrix} y_{1}\\ y_{2}\\ y_{3}\\ \vdots \\ y_{N-2}\\ y_{N-1}\\ y_{N}\end{bmatrix}\\ Ay &=\alpha By \end{align*}

Where $$\alpha =\lambda ^{2}$$ and $$B=-h^{2}\mathbf{I}$$. The above is implemented in Matlab and eig is used to ﬁnd $$\alpha$$.

The analytical value of the eigenvalue is given from$\lambda _{n}=\sqrt{\frac{P_{n}b^{4}}{EL^{2}I_{0}}}$ Where$P_{n}=n^{2}\pi ^{2}\left ( \frac{a}{b}\right ) ^{2}\frac{EI_{0}}{L^{2}}$ Hence$\lambda _{n}=\sqrt{\frac{n^{2}\pi ^{2}\left ( \frac{a}{b}\right ) ^{2}\frac{EI_{0}}{L^{2}}b^{4}}{EL^{2}I_{0}}}=\sqrt{\frac{n^{2}\pi ^{2}a^{2}b^{2}}{L^{4}}}$ Using $$a=3,b=6,L=3,$$the ﬁrst three eigenvalues are \begin{align*} \lambda _{1} & =\sqrt{\frac{\pi ^{2}\left ( 3^{2}\right ) \left ( 6^{2}\right ) }{\left ( 3^{4}\right ) }}=6.2832\\ \lambda _{2} & =\sqrt{\frac{2^{2}\pi ^{2}\left ( 3^{2}\right ) \left ( 6^{2}\right ) }{\left ( 3^{4}\right ) }}=12.566\\ \lambda _{3} & =\sqrt{\frac{3^{2}\pi ^{2}\left ( 3^{2}\right ) \left ( 6^{2}\right ) }{\left ( 3^{4}\right ) }}=18.850 \end{align*}

And the corresponding analytical mode shapes, using $$A_{n}=1$$ when normalized is\begin{align*} y_{1}\left ( z\right ) & =z\sin \left ( \pi \frac{b}{L}\left ( 1-\frac{a}{Lz}\right ) \right ) =z\sin \left ( 2\pi \left ( 1-\frac{1}{z}\right ) \right ) \\ y_{2}\left ( x\right ) & =z\sin \left ( 2\pi \frac{b}{L}\left ( 1-\frac{a}{Lz}\right ) \right ) =z\sin \left ( 4\pi \left ( 1-\frac{1}{z}\right ) \right ) \\ y_{3}\left ( x\right ) & =z\sin \left ( 3\pi \frac{b}{L}\left ( 1-\frac{a}{Lz}\right ) \right ) =z\sin \left ( 6\pi \left ( 1-\frac{1}{z}\right ) \right ) \end{align*}

And the corresponding buckling loads at each mode shape are, using $$E=10^{9}$$Pa, $$I_{0}=\frac{1}{4}\pi \left ( r_{b}\right ) ^{4}$$ where $$r_{b}=0.2$$ meter are$P_{n}=n^{2}\pi ^{2}\left ( \frac{a}{b}\right ) ^{2}\frac{EI_{0}}{L^{2}}=n^{2}\pi ^{2}\left ( \frac{3}{6}\right ) ^{2}\frac{\left ( 10^{9}\right ) \frac{1}{4}\pi \left ( 0.2\right ) ^{4}}{\left ( 3^{2}\right ) }=n^{2}\left ( 3.445\,1\times 10^{5}\right ) \text{ N}$ Hence\begin{align*} P_{1} & =3.445\,1\times 10^{5}\text{ N}\\ P_{2} & =4\left ( 3.445\,1\times 10^{5}\right ) =1.378\,1\times 10^{6}\text{ N}\\ P_{3} & =9\left ( 3.445\,1\times 10^{5}\right ) =3.100\,6\times 10^{6}\text{ N} \end{align*}

These are used to compare the numerical solutions from bvp4c, eig and power method against. (for power method, only the lowest eigenvalue is obtained). For the numerical computation of $$P_{n}$$, after ﬁnding the numerical eigenvalue $$\lambda _{n}$$, then $$P_{n}$$ is found from$P_{n}=\frac{\lambda _{n}^{2}EL^{2}I_{0}}{b^{4}}$ And the values obtained are compared to the analytical $$P_{n}$$. The following plots show the result for the ﬁrst three eigenvalues and eigenfunctions found.

##### 2.3.2.1 Power method

For the power method, the $$A$$ matrix is setup a little diﬀerent than with the above eig method. Starting from

$\left . z^{4}y^{\prime \prime }+\lambda ^{2}y\right \vert _{i}\approx \left ( 1+ih\right ) ^{4}\frac{y_{i+1}-2y_{i}+y_{i-1}}{h^{2}}+\lambda ^{2}y_{i}$

Hence\begin{align*} \left ( 1+ih\right ) ^{4}\left ( y_{i+1}-2y_{i}+y_{i-1}\right ) +h^{2}\lambda ^{2}y_{i} & =0\\ \frac{-\left ( 1+ih\right ) ^{4}y_{i+1}+2\left ( 1+ih\right ) ^{4}y_{i}-\left ( 1+ih\right ) ^{4}y_{i-1}}{h^{2}} & =\lambda ^{2}y_{i} \end{align*}

At node $$i=1$$,$\frac{-\left ( 1+h\right ) ^{4}y_{2}+2\left ( 1+h\right ) ^{4}y_{1}-\left ( 1+h\right ) ^{4}y_{0}}{h^{2}}=\lambda ^{2}y_{1}$

Since $$y_{0}=0$$ then

$\frac{-\left ( 1+h\right ) ^{4}y_{2}+2\left ( 1+h\right ) ^{4}y_{1}}{h^{2}}=\lambda ^{2}y_{1}$

At node $$i=2$$$\frac{-\left ( 1+2h\right ) ^{4}y_{3}+2\left ( 1+2h\right ) ^{4}y_{2}-\left ( 1+2h\right ) ^{4}y_{1}}{h^{2}}=\lambda ^{2}y_{2}$ And so on. At the last node, $$i=N$$$\frac{-\left ( 1+Nh\right ) ^{4}y_{N+1}+2\left ( 1+Nh\right ) ^{4}y_{N}-\left ( 1+Nh\right ) ^{4}y_{N-1}}{h^{2}}=\lambda ^{2}y_{N}$

Since $$y_{N+1}=0$$ then

$\frac{2\left ( 1+Nh\right ) ^{4}y_{N}-\left ( 1+Nh\right ) ^{4}y_{N-1}}{h^{2}}=\lambda ^{2}y_{N}$

At $$i=N-1$$$\frac{-\left ( 1+\left ( N-1\right ) h\right ) ^{4}y_{N}+2\left ( 1+\left ( N-1\right ) h\right ) ^{4}y_{\left ( N-1\right ) }-\left ( 1+\left ( N-1\right ) h\right ) ^{4}y_{N-2}}{h^{2}}=\lambda ^{2}y_{\left ( N-1\right ) }$ Hence the structure is\begin{align*} \begin{bmatrix} \frac{2\left ( 1+h\right ) ^{4}}{h^{2}} & \frac{-\left ( 1+h\right ) ^{4}}{h^{2}} & 0 & 0 & 0 & \cdots & 0\\ \frac{-\left ( 1+2h\right ) ^{4}}{h^{2}} & \frac{2\left ( 1+2h\right ) ^{4}}{h^{2}} & \frac{-\left ( 1+2h\right ) ^{4}}{h^{2}} & 0 & 0 & \cdots & \vdots \\ 0 & \frac{-\left ( 1+3h\right ) ^{4}}{h^{2}} & \frac{2\left ( 1+3h\right ) ^{4}}{h^{2}} & \frac{-\left ( 1+3h\right ) ^{4}}{h^{2}} & 0 & \cdots & \vdots \\ 0 & 0 & \cdots & \ddots & \cdots & \cdots & \vdots \\ \vdots & \cdots & \cdots & \cdots & \ddots & \cdots & 0\\ \vdots & \cdots & \cdots & \cdots & \frac{-\left ( 1+\left ( N-1\right ) h\right ) ^{4}}{h^{2}} & \frac{2\left ( 1+\left ( N-1\right ) h\right ) ^{4}}{h^{2}} & \frac{-\left ( 1+\left ( N-1\right ) h\right ) ^{4}}{h^{2}}\\ 0 & \cdots & \cdots & \cdots & 0 & \frac{-\left ( 1+Nh\right ) ^{4}}{h^{2}} & \frac{2\left ( 1+Nh\right ) ^{4}}{h^{2}}\end{bmatrix}\begin{bmatrix} y_{1}\\ y_{2}\\ y_{3}\\ \vdots \\ y_{N-2}\\ y_{N-1}\\ y_{N}\end{bmatrix} & =\lambda ^{2}\begin{bmatrix} 1 & 0 & 0 & 0 & 0 & \cdots & 0\\ 0 & 1 & 0 & 0 & 0 & \cdots & \vdots \\ 0 & 0 & 1 & 0 & 0 & \cdots & \vdots \\ 0 & 0 & 0 & 1 & 0 & \cdots & \vdots \\ \vdots & \cdots & \cdots & \cdots & \ddots & \cdots & \vdots \\ \vdots & \vdots & \vdots & \vdots & \vdots & \ddots & 0\\ 0 & 0 & 0 & 0 & \cdots & 0 & 1 \end{bmatrix}\begin{bmatrix} y_{1}\\ y_{2}\\ y_{3}\\ \vdots \\ y_{N-2}\\ y_{N-1}\\ y_{N}\end{bmatrix} \end{align*}

$Ay =\lambda ^{2}y$

The above structure is now used to solve for lowest eigenvalue and corresponding eigenvector.

##### 2.3.2.2 Results

Each mode shape plot is given, showing the eigenvalue produced by each solver and the initial mode shape guess used. There are 3 plots, one for each mode shape. The ﬁrst, second and third. (the problem asked for only the ﬁrst mode shape, but the second and third were added for veriﬁcation). For power method, only the lowest eigenvalue and corresponding eigenvector are found.

1.
First mode shape
 Solver eigenvalue found $$\lambda _{n}$$ Corresponding Critical load $$P_{n}$$ (N) analytical $$6.2817063$$ $$344352.012$$ bvp4c $$6.2821629$$ $$344402.076$$ Matlab eig $$6.2817063$$ $$344352.012$$ Power method $$6.2817055$$ $$344351.929$$

2.
Second mode shape
 Solver eigenvalue found $$\lambda _{n}$$ Corresponding Critical load $$P_{n}$$ (N) analytical $$12.5663706$$ $$1378056.741$$ bvp4c $$12.5663983$$ $$1378062.820$$ eig $$12.5534143$$ $$1375216.578$$

3.
Third mode shape
 Solver eigenvalue found $$\lambda _{n}$$ Corresponding Critical load $$P_{n}$$ (N) analytical $$18.8495559$$ $$3100627.668$$ bvp4c $$18.8499237$$ $$3100748.676$$ eig $$18.80506$$ $$3086006.365$$

Printout of Matlab console running the program

>>nma_HW3_EMA_471_problem_2
*********************
running mode 1
Eigenvalue obtained with bvp4c, is 6.2821629
eigenvalue from eig is 6.2817063
eigenvalue from analytical is 6.2831853
critical load from analytical is 344514.185
eigenvalue obtained with the power iteration method 6.2817055
*********************
running mode 2
Eigenvalue obtained with bvp4c, is 12.5663983
eigenvalue from eig is 12.5534143
eigenvalue from analytical is 12.5663706
critical load from analytical is 1378056.741
*********************
running mode 3
Eigenvalue obtained with bvp4c, is 18.8499237
eigenvalue from eig is 18.8050600
eigenvalue from analytical is 18.8495559
critical load from analytical is 3100627.668

#### 2.3.3 Problem 3

2.3.3.1 Power method
2.3.3.2 Results
2.3.3.3 Source code

\begin{align*} \frac{d^{2}\theta }{dz^{2}}+\lambda ^{2}z\theta & =0\\ \theta ^{\prime }\left ( 1\right ) & =0\\ \theta \left ( 0\right ) & =0 \end{align*}

For domain $$0\leq z\leq 1$$. By numerically solving for the lowest eigenvalue $$\lambda _{1}$$, the limiting height $$L$$ can next be found from solving for $$L$$ in $$\lambda ^{2}=\frac{\rho gAL^{3}}{EI}$$. Three methods are used to ﬁnd $$\lambda _{1}$$: Power method, bpv4c and Matlab eig. The buckled shape (eigen shapes) found from the numerical method is compared to the analytical shape given $\theta _{1}\left ( z\right ) =A_{1}\sqrt{z}J_{\left ( -\frac{1}{3}\right ) }\left ( \frac{2}{3}\lambda _{1}z^{\frac{3}{2}}\right )$ $$A_{1}$$ is taken as $$1$$ due to the normalization used and $$J$$ is the Bessel function of ﬁrst kind.

The ﬁrst step is to convert the ODE into state space for use with bvp4c. Let $$x_{1}=\theta ,x_{2}=\theta ^{\prime }$$. Taking derivatives gives\begin{align*} \dot{x}_{1} & =x_{2}\\ \dot{x}_{2} & =-\lambda ^{2}zx_{1} \end{align*}

For using eig, the problem needs to discretized ﬁrst. The following shows the grid used

The grid starts at $$i=0$$ which corresponds to $$z=0$$ and ends at $$i=N+1$$ which corresponds to $$z=1$$. Since $$\theta$$ is not known at $$z=0$$, then in this problem $$i=0$$ is included in the internal grid points, hence the $$A$$ matrix will have size $$(N+1)\times (N+1)$$. Using second order centered diﬀerence gives$\left . \frac{d^{2}\theta }{dz^{2}}\right \vert _{i}=\frac{\theta _{i+1}-2\theta _{i}+\theta _{i-1}}{h^{2}}$ Therefore, the approximation to the diﬀerential equation at grid $$i$$ (on the internal nodes as shown in the above diagram) is as follows.$\left . \frac{1}{z}\frac{d^{2}\theta }{dz^{2}}+\lambda ^{2}\theta =0\right \vert _{i}\approx \frac{1}{ih}\frac{\theta _{i+1}-2\theta _{i}+\theta _{i-1}}{h^{2}}+\lambda ^{2}\theta _{i}$ Hence\begin{align*} \frac{1}{ih}\frac{\theta _{i+1}-2\theta _{i}+\theta _{i-1}}{h^{2}}+\lambda ^{2}\theta _{i} & =0\\ \frac{1}{ih}\left ( \theta _{i+1}-2\theta _{i}+\theta _{i-1}\right ) & =-h^{2}\lambda ^{2}\theta _{i} \end{align*}

At node $$i=0$$\begin{align*} \frac{\theta _{1}-2\theta _{0}+\theta _{-1}}{ih+\varepsilon } & =-h^{2}\lambda ^{2}\theta _{0}\\ \frac{\theta _{1}-2\theta _{0}+\theta _{-1}}{\varepsilon } & =-h^{2}\lambda ^{2}\theta _{0} \end{align*}

Where $$\varepsilon$$ is small value $$10^{-6}$$ in order to handle the condition at $$z=0$$.

To ﬁnd $$\theta _{i=-1}$$, the condition $$\theta ^{\prime }\left ( 0\right ) =0$$ is used. Since $$\theta ^{\prime }\left ( 0\right ) =\frac{\theta _{1}-\theta _{-1}}{2h}=0$$ then $$\theta _{-1}=\theta _{1}$$ and the above becomes$\fbox{\frac{2\theta _1-2\theta _0}{\varepsilon }=-h^2\lambda ^2\theta _0}$ At $$i=1$$$\fbox{\frac{\theta _2-2\theta _1+\theta _0}{h}=-h^2\lambda ^2\theta _1}$ At node $$i=2$$$\fbox{\frac{\theta _3-2\theta _2+\theta _1}{2h}=-h^2\lambda ^2\theta _2}$ And so on. At the last internal node, $$i=N$$$\frac{\theta _{N+1}-2\theta _{N}+\theta _{N-1}}{Nh}=-h^{2}\lambda ^{2}\theta _{N}$ But $$\theta _{N+1}=0$$ from boundary conditions, hence$\fbox{\frac{-2\theta _N+\theta _N-1}{Nh}=-h^2\lambda ^2\theta _N}$ At $$i=N-1$$$\fbox{\frac{\theta _N-2\theta _N-1+\theta _N-2}{\left ( N-1\right ) h}=-h^2\lambda ^2\theta _N-1}$ Hence the structure is\begin{align*} \begin{bmatrix} -\frac{2}{\varepsilon } & \frac{2}{\varepsilon } & 0 & 0 & 0 & \cdots & 0\\ \frac{1}{h} & -\frac{2}{h} & \frac{1}{h} & 0 & 0 & \cdots & \vdots \\ 0 & \frac{1}{2h} & -\frac{2}{2h} & \frac{1}{2h} & 0 & \cdots & \vdots \\ 0 & 0 & \cdots & \ddots & \cdots & \cdots & \vdots \\ \vdots & \cdots & \cdots & \cdots & \ddots & \cdots & 0\\ \vdots & \cdots & \cdots & \cdots & \frac{1}{\left ( N-1\right ) h} & -\frac{2}{\left ( N-1\right ) h} & \frac{1}{\left ( N-1\right ) h}\\ 0 & \cdots & \cdots & \cdots & 0 & \frac{1}{Nh} & -\frac{2}{Nh}\end{bmatrix}\begin{bmatrix} \theta _{0}\\ \theta _{1}\\ \theta _{2}\\ \vdots \\ \theta _{N-2}\\ \theta _{N-1}\\ \theta _{N}\end{bmatrix} & =-h^{2}\lambda ^{2}\begin{bmatrix} 1 & 0 & 0 & 0 & 0 & \cdots & 0\\ 0 & 1 & 0 & 0 & 0 & \cdots & \vdots \\ 0 & 0 & 1 & 0 & 0 & \cdots & \vdots \\ 0 & 0 & 0 & 1 & 0 & \cdots & \vdots \\ \vdots & \cdots & \cdots & \cdots & \ddots & \cdots & \vdots \\ \vdots & \vdots & \vdots & \vdots & \vdots & 1 & 0\\ 0 & 0 & 0 & 0 & \cdots & 0 & 1 \end{bmatrix}\begin{bmatrix} \theta _{0}\\ \theta _{1}\\ \theta _{2}\\ \vdots \\ \theta _{N-2}\\ \theta _{N-1}\\ \theta _{N}\end{bmatrix} \\ A\theta & =\alpha B\theta \end{align*}

Where $$\alpha =\lambda ^{2}$$ and $$B=-h^{2}\mathbf{I}$$. The above is implemented in Matlab and eig is used to ﬁnd $$\alpha$$.

##### 2.3.3.1 Power method

For the power method, the $$A$$ matrix is setup a little diﬀerent than with the above eig method which results in\begin{align*} \frac{-1}{h^{2}}\begin{bmatrix} -\frac{2}{\varepsilon } & \frac{2}{\varepsilon } & 0 & 0 & 0 & \cdots & 0\\ \frac{1}{h} & -\frac{2}{h} & \frac{1}{h} & 0 & 0 & \cdots & \vdots \\ 0 & \frac{1}{2h} & -\frac{2}{2h} & \frac{1}{2h} & 0 & \cdots & \vdots \\ 0 & 0 & \cdots & \ddots & \cdots & \cdots & \vdots \\ \vdots & \cdots & \cdots & \cdots & \ddots & \cdots & 0\\ \vdots & \cdots & \cdots & \cdots & \frac{1}{\left ( N-1\right ) h} & -\frac{2}{\left ( N-1\right ) h} & \frac{1}{\left ( N-1\right ) h}\\ 0 & \cdots & \cdots & \cdots & 0 & \frac{1}{Nh} & -\frac{2}{Nh}\end{bmatrix}\begin{bmatrix} \theta _{0}\\ \theta _{1}\\ \theta _{2}\\ \vdots \\ \theta _{N-2}\\ \theta _{N-1}\\ \theta _{N}\end{bmatrix} & =\lambda ^{2}\begin{bmatrix} 1 & 0 & 0 & 0 & 0 & \cdots & 0\\ 0 & 1 & 0 & 0 & 0 & \cdots & \vdots \\ 0 & 0 & 1 & 0 & 0 & \cdots & \vdots \\ 0 & 0 & 0 & 1 & 0 & \cdots & \vdots \\ \vdots & \cdots & \cdots & \cdots & \ddots & \cdots & \vdots \\ \vdots & \vdots & \vdots & \vdots & \vdots & 1 & 0\\ 0 & 0 & 0 & 0 & \cdots & 0 & 1 \end{bmatrix}\begin{bmatrix} \theta _{0}\\ \theta _{1}\\ \theta _{2}\\ \vdots \\ \theta _{N-2}\\ \theta _{N-1}\\ \theta _{N}\end{bmatrix} \\ A\theta & =\lambda ^{2}\theta \end{align*}

The above structure is now used to solve for lowest eigenvalue and corresponding eigenvector.

One the system is solved for the lowest eigenvalue, the critical length of the column is found by solving for $$L$$ from $$\lambda ^{2}=\frac{\rho gAL^{3}}{EI}$$.

##### 2.3.3.2 Results

The following table shows the lowest eigenvalue found by each method, and the corresponding $$L$$ found.

 method $$\lambda$$ $$L_{critical}$$ meter bvp4c $$2.7995616$$ $$18.81235953$$ eig $$2.71870438$$ $$18.44836607$$ power $$2.71870430$$ $$18.44836567$$

The following are the plots of the mode shape by each method. There is little diﬀerence that can be seen between the eig and the power methods since they are both based on the same ﬁnite diﬀerence scheme. The bvp4c is the most similar to the analytical solution. In order to evaluate and plot the analytical solution given in the problem, the eigenvalue found from bvp4c was used.

The following plot shows the result on one plot for all the methods. As can be seen, they are very similar to each others.

Below is a zoomed version, showing the bvp4c is in very good agreement with the analytical plot. The power method and the eig methods are almost exactly the same. All methods become very close to each others at the boundaries and they are most diﬀerent in the middle of the range.

The following shows the result in separate plots

The following is printout of Matlab console running the program

>>nma_HW3_EMA_471_problem_3
*********************
Eigenvalue obtained with bvp4c is
2.79956162718772

Critical length is
18.8123595369211
*********************
eigenvalue from eig is
2.71870438941484

Critical length is
18.4483660724537
*********************
eigenvalue obtained with the power iteration method
2.7187043018523

Critical length is
18.4483656763372