From 509b7cf4564756bb37ea9dd2ee0dcebeccaa715a Mon Sep 17 00:00:00 2001 From: Kyle Niemeyer Date: Wed, 26 Feb 2020 17:54:20 -0800 Subject: [PATCH] more on eigenvalue module --- _build/bvps/eigenvalue.html | 121 +++++++++++++++++++++++++++++++++- content/bvps/eigenvalue.ipynb | 105 +++++++++++++++++++++++++++++ 2 files changed, 225 insertions(+), 1 deletion(-) diff --git a/_build/bvps/eigenvalue.html b/_build/bvps/eigenvalue.html index b30aa4d..71cea06 100644 --- a/_build/bvps/eigenvalue.html +++ b/_build/bvps/eigenvalue.html @@ -11,7 +11,7 @@ next_page: url: /quizzes/quiz2-IVPs.html suffix: .ipynb -search: y lambda l frac x n align begin p end beam pi equation b cos sin quad boundary ei rightarrow ldots infty different solution lets conditions mz modes gather eigenvalue problems eigenvalues example buckling consider e neq load means characteristic value where obtain certain values us deflection supported get governing considering sum general because text eigenfunctions associated corresponding buckle properties cr case types show areas involving not able analytical identify tell important information system simply static vertical start moments around upper pin m also d dx governs stability small deflections simplify things define gives ode coefficients apply starting otherwise trivial +search: y lambda begin align end x l k frac equation n pi eigenvalues p beam delta b boundary quad rightarrow bmatrix solution conditions ei cos sin lets case our values get ldots infty left right different equations mathbf det eigenvalue problems where buckling consider deflection mz modes gather load yi obtain us system example supported also e ode general neq because instead need associated represent corresponding cr finite above using matrix given means characteristic value not analytical certain simply governing considering sum simplify trivial text problem infinite eigenfunctions three buckle recall properties sqrt critical slightly same differences points into modify comment: "***PROGRAMMATICALLY GENERATED, DO NOT EDIT. SEE ORIGINAL FILES IN /content***" --- @@ -156,12 +156,131 @@

Example: beam \text{so} \quad \lambda L &= \frac{(2n-1) \pi}{2} \quad n = 1,2,3,\ldots, \infty \\ \lambda &= \frac{(2n-1) \pi}{2 L} \quad n = 1,2,3,\ldots, \infty \end{align}

+

Then, the critical buckling load, again corresponding to $n=1$, is +\begin{equation} +P_{cr} = \frac{\pi^2 EI}{4 L^2} +\end{equation}

+ + + + + + +
+ +
+
+

Getting eigenvalues numerically

We can only get the eigenvalues analytically if we can obtain an analytical solution of the ODE, but we might want to get eigenvalues for more complex problems too. In that case, we can use an approach based on finite differences to find the eigenvalues.

+

Consider the same problem as above, for deflection of a simply supported beam: +\begin{equation} +y'' + \lambda^2 y = 0 +\end{equation} +with boundary conditions $y(0) = 0$ and $y(L) = 0$. Let's represent this using finite differences, for a case where $L=3$ and $\Delta x = 1$, so we have four points in our solution grid.

+

The finite difference representation of the ODE is: +\begin{align} +\frac{y_{i-1} - 2y_i + y_{i+1}}{\Delta x^2} + \lambda^2 y_i &= 0 \\ +y_{i-1} + \left( \lambda^2 \Delta x^2 - 2 \right) y_i + y_{i+1} &= 0 +\end{align} +However, in this case, we are not solving for the values of deflection ($y_i$), but instead the eigenvalues $\lambda$.

+

Then, we can write the system of equations using the above recursion formula and our two boundary conditions: +\begin{align} +y_1 &= 0 \\ +y_1 + y_2 \left( \lambda^2 \Delta x^2 - 2 \right) + y_3 &= 0 \\ +y_2 + y_3 \left( \lambda^2 \Delta x^2 - 2 \right) + y_4 &= 0 \\ +y_4 &= 0 +\end{align} +which we can simplify down to two equations by incorporating the boundary conditions into the equations for the two middle points, and also letting $k = \lambda^2 \Delta x^2$: +\begin{align} +y_2 (k-2) + y_3 &= 0 \\ +y_2 + y_3 (k-2) &= 0 +\end{align} +Let's modify this once more by multiplying both equations by $-1$: +\begin{align} +y_2 (2-k) - y_3 &= 0 \\ +-y_2 + y_3 (2-k) &= 0 +\end{align}

+

Now we can represent this system of equations as a matrix equation $A \mathbf{y} = \mathbf{b} = \mathbf{0}$: +\begin{equation} +\begin{bmatrix} 2-k & -1 \\ -1 & 2-k \end{bmatrix} +\begin{bmatrix} y_2 \\ y_3 \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \end{bmatrix} +\end{equation} +$\mathbf{y} = \mathbf{0}$ is a trivial solution to this, so instead $\det(A) = 0$ satisfies this equation. +For our $2\times 2$ matrix, that looks like: +\begin{align} +\det(A) = \begin{vmatrix} 2-k & -1 \\ -1 & 2-k \end{vmatrix} = (2-k)^2 - 1 &= 0 \\ +k^2 - 4k + 3 &= 0 \\ +(k-3)(k-1) &= 0 +\end{align} +so the roots of this equation are $k_1 = 1$ and $k_2 = 3$. Recall that $k$ is directly related to the eigenvalue: $k = \lambda^2 \Delta x^2$, and $\Delta x = 1$ for this case, so we can calculate the two associated eigenvalues: +\begin{align} +k_1 &= \lambda_1^2 \Delta x^2 = 1 \rightarrow \lambda_1 = 1 \\ +k_2 &= \lambda_2^2 \Delta x^2 = 3 \rightarrow \lambda_2 = \sqrt{3} = 1.732 +\end{align}

+

Our work has given us approximations for the first two eigenvalues. We can compare these against the exact values, given in general by $\lambda = n \pi / L$ (which we determined above): +\begin{align} +n=1: \quad \lambda_1 &= \frac{\pi}{L} = \frac{\pi}{3} = 1.0472 \\ +n=2: \quad \lambda_2 &= \frac{2\pi}{L} = \frac{2\pi}{3} = 2.0944 +\end{align} +So, our approximations are close, but with some obvious error. This is because we used a fairly crude step size of $\Delta x = 1$, dividing the domain into just three segments. By using a finer resolution, we can get more-accurate eigenvalues and also more of them (remember, there are actually an infinite number!).

+

To do that, we will need to use Matlab, which offers the eig() function for calculating eigenvalues---essentially it is finding the roots to the polynomial given by $\det(A) = 0$. We need to modify this slightly, though, to use the function: +\begin{align} +\det(A) &= 0 \\ +\det \left( A^* - k I \right) = 0 +\end{align} +where the new matrix is +\begin{equation} +A^* = \begin{bmatrix} 2 & -1 \\ -1 & 2 \end{bmatrix} +\end{equation} +Then, eig(A*) will provide the values of $k$, which we can use to find the $\lambda$s:

+
+ +
+
+ +
+
+
clear all; clc
+
+dx = 1.0;
+L = 3.0;
+
+Astar = [2 -1; -1 2];
+k = eig(Astar);
+
+lambda = sqrt(k) / dx^2;
+
+fprintf('lambda 1: %6.3f\n', lambda(1));
+fprintf('lambda 2: %6.3f', lambda(2));
+
+ +
+
+
+ +
+
+ +
+
+ +
+
lambda 1:  1.000
+lambda 2:  1.732
+
+
+
+
+
+ +
+
+ diff --git a/content/bvps/eigenvalue.ipynb b/content/bvps/eigenvalue.ipynb index 34f749c..23d77ac 100644 --- a/content/bvps/eigenvalue.ipynb +++ b/content/bvps/eigenvalue.ipynb @@ -134,6 +134,111 @@ "\\end{equation}" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Getting eigenvalues numerically\n", + "\n", + "We can only get the eigenvalues analytically if we can obtain an analytical solution of the ODE, but we might want to get eigenvalues for more complex problems too. In that case, we can use an approach based on *finite differences* to find the eigenvalues.\n", + "\n", + "Consider the same problem as above, for deflection of a simply supported beam:\n", + "\\begin{equation}\n", + "y'' + \\lambda^2 y = 0 \n", + "\\end{equation}\n", + "with boundary conditions $y(0) = 0$ and $y(L) = 0$. Let's represent this using finite differences, for a case where $L=3$ and $\\Delta x = 1$, so we have four points in our solution grid.\n", + "\n", + "The finite difference representation of the ODE is:\n", + "\\begin{align}\n", + "\\frac{y_{i-1} - 2y_i + y_{i+1}}{\\Delta x^2} + \\lambda^2 y_i &= 0 \\\\\n", + "y_{i-1} + \\left( \\lambda^2 \\Delta x^2 - 2 \\right) y_i + y_{i+1} &= 0\n", + "\\end{align}\n", + "However, in this case, we are not solving for the values of deflection ($y_i$), but instead the **eigenvalues** $\\lambda$.\n", + "\n", + "Then, we can write the system of equations using the above recursion formula and our two boundary conditions:\n", + "\\begin{align}\n", + "y_1 &= 0 \\\\\n", + "y_1 + y_2 \\left( \\lambda^2 \\Delta x^2 - 2 \\right) + y_3 &= 0 \\\\\n", + "y_2 + y_3 \\left( \\lambda^2 \\Delta x^2 - 2 \\right) + y_4 &= 0 \\\\\n", + "y_4 &= 0\n", + "\\end{align}\n", + "which we can simplify down to two equations by incorporating the boundary conditions into the equations for the two middle points, and also letting $k = \\lambda^2 \\Delta x^2$:\n", + "\\begin{align}\n", + "y_2 (k-2) + y_3 &= 0 \\\\\n", + "y_2 + y_3 (k-2) &= 0\n", + "\\end{align}\n", + "Let's modify this once more by multiplying both equations by $-1$:\n", + "\\begin{align}\n", + "y_2 (2-k) - y_3 &= 0 \\\\\n", + "-y_2 + y_3 (2-k) &= 0\n", + "\\end{align}\n", + "\n", + "Now we can represent this system of equations as a matrix equation $A \\mathbf{y} = \\mathbf{b} = \\mathbf{0}$:\n", + "\\begin{equation}\n", + "\\begin{bmatrix} 2-k & -1 \\\\ -1 & 2-k \\end{bmatrix}\n", + "\\begin{bmatrix} y_2 \\\\ y_3 \\end{bmatrix} = \\begin{bmatrix} 0 \\\\ 0 \\end{bmatrix}\n", + "\\end{equation}\n", + "$\\mathbf{y} = \\mathbf{0}$ is a trivial solution to this, so instead $\\det(A) = 0$ satisfies this equation.\n", + "For our $2\\times 2$ matrix, that looks like:\n", + "\\begin{align}\n", + "\\det(A) = \\begin{vmatrix} 2-k & -1 \\\\ -1 & 2-k \\end{vmatrix} = (2-k)^2 - 1 &= 0 \\\\\n", + "k^2 - 4k + 3 &= 0 \\\\\n", + "(k-3)(k-1) &= 0\n", + "\\end{align}\n", + "so the roots of this equation are $k_1 = 1$ and $k_2 = 3$. Recall that $k$ is directly related to the eigenvalue: $k = \\lambda^2 \\Delta x^2$, and $\\Delta x = 1$ for this case, so we can calculate the two associated eigenvalues:\n", + "\\begin{align}\n", + "k_1 &= \\lambda_1^2 \\Delta x^2 = 1 \\rightarrow \\lambda_1 = 1 \\\\\n", + "k_2 &= \\lambda_2^2 \\Delta x^2 = 3 \\rightarrow \\lambda_2 = \\sqrt{3} = 1.732\n", + "\\end{align}\n", + "\n", + "Our work has given us approximations for the first two eigenvalues. We can compare these against the exact values, given in general by $\\lambda = n \\pi / L$ (which we determined above):\n", + "\\begin{align}\n", + "n=1: \\quad \\lambda_1 &= \\frac{\\pi}{L} = \\frac{\\pi}{3} = 1.0472 \\\\\n", + "n=2: \\quad \\lambda_2 &= \\frac{2\\pi}{L} = \\frac{2\\pi}{3} = 2.0944\n", + "\\end{align}\n", + "So, our approximations are close, but with some obvious error. This is because we used a fairly crude step size of $\\Delta x = 1$, dividing the domain into just three segments. By using a finer resolution, we can get more-accurate eigenvalues and also more of them (remember, there are actually an infinite number!). \n", + "\n", + "To do that, we will need to use Matlab, which offers the `eig()` function for calculating eigenvalues---essentially it is finding the roots to the polynomial given by $\\det(A) = 0$. We need to modify this slightly, though, to use the function:\n", + "\\begin{align}\n", + "\\det(A) &= 0 \\\\\n", + "\\det \\left( A^* - k I \\right) = 0\n", + "\\end{align}\n", + "where the new matrix is\n", + "\\begin{equation}\n", + "A^* = \\begin{bmatrix} 2 & -1 \\\\ -1 & 2 \\end{bmatrix}\n", + "\\end{equation}\n", + "Then, `eig(A*)` will provide the values of $k$, which we can use to find the $\\lambda$s:" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "lambda 1: 1.000\n", + "lambda 2: 1.732" + ] + } + ], + "source": [ + "clear all; clc\n", + "\n", + "dx = 1.0;\n", + "L = 3.0;\n", + "\n", + "Astar = [2 -1; -1 2];\n", + "k = eig(Astar);\n", + "\n", + "lambda = sqrt(k) / dx^2;\n", + "\n", + "fprintf('lambda 1: %6.3f\\n', lambda(1));\n", + "fprintf('lambda 2: %6.3f', lambda(2));" + ] + }, { "cell_type": "code", "execution_count": null,