The ODE is\[ y^{\prime \prime }+2y^{\prime }+\lambda ^{2}y=0 \] With \(y\left ( 0\right ) =y\left ( 1\right ) =0\). The first step is to find 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 find the eigenvalues. Using second order centered difference 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 differential 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 find \(\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 first 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 first three eigenvalues and eigenfunctions found. The main difficulty with using bvp4c for solving the eigenvalue problem is on deciding which guess \(\lambda \) to use for each mode shape to solve for. The first 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.
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 first, second and third. (the problem asked for only the first two mode shapes, but the third one was added for verification).
| Solver | eigenvalue found |
| analytical | \(\lambda _{1}=\sqrt{\pi ^{2}+1}=\) \(3.296\,9\) |
| bvp4c | \(3.2969\) |
| eig | \(3.2960997\) |
| Solver | eigenvalue found |
| analytical | \(\lambda _{1}=\sqrt{\left ( 2\pi \right ) ^{2}+1}=\) \(6.3623\) |
| bvp4c | \(6.3622\) |
| eig | \(6.35738\) |
| 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
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 first 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 first. The following shows the grid used
The grid starts at \(z=1\) and ends at \(z=2\) since this is the domain of the differential 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 difference 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 differential 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 find \(\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 first 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 finding 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 first three eigenvalues and eigenfunctions found.
For the power method, the \(A\) matrix is setup a little different 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.
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 first, second and third. (the problem asked for only the first mode shape, but the second and third were added for verification). For power method, only the lowest eigenvalue and corresponding eigenvector are found.
| 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\) |
| 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\) |
| 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
\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 find \(\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 first kind.
The first 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 first. 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 difference 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 differential 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 find \(\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 find \(\alpha \).
For the power method, the \(A\) matrix is setup a little different 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}\).
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 difference that can be seen between the eig and the power methods since they are both based on the same finite difference 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 different in the middle of the range.
The following shows the result in separate plots
The following is printout of Matlab console running the program