{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Lecture 12 - applications of eigenvalue problems\n", "\n", "This notebook presents two applications of eigenvalue problems. The examples were covered in the lectures, and are presented in more depth here. The two examples are \n", "\n", "1. The PageRank algorithm for racking the importance of web pages\n", "1. Vibration of a multi degree-of-freedom mass-spring system" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import sys\n", "if 'conda' in sys.version: # Install using conda if we're using Anaconda Python\n", " !conda install -y pygraphviz" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## PageRank\n", "\n", "PageRank is a algorithm used by Google to rank the 'importance' of a web page. The concept is that every web page has an importance score, and it 'endows' its importance evenly amongst the pages to which it links. The will lead to a problem where the ranking of all pages is an eigenvalue problem. We will demonstrate the process for the small web of four pages.\n", "\n", "This is a larger version of the the example from the lectures." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Model problem\n", "\n", "Consider a web of five pages, $p_{0}$, $p_{1}$, $p_{2}$, $p_{3}$ and $p_{4}$. Consider the scenario:\n", "\n", "- $p_{0}$ links to: $p_{1}$, $p_{2}$, $p_{3}$ and $p_{4}$\n", "- $p_{1}$ links to: $p_{2}$, $p_{3}$ and $p_{4}$\n", "- $p_{2}$ links to: $p_{0}$ and $p_{1}$\n", "- $p_{3}$ links to: $p_{0}$ and $p_{2}$\n", "- $p_{4}$ links to: $p_{1}$ and $p_{2}$ and $p_{3}$\n", "\n", "We can build a directed graph to describe the connections, and for this we use the package `networkx` (https://networkx.github.io/). To visualise the graph we use Matplotlib." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "\n", "\n", "\n", "\n", "\n", "p0\n", "\n", "p0\n", "\n", "\n", "\n", "p1\n", "\n", "p1\n", "\n", "\n", "\n", "p0->p1\n", "\n", "\n", "1/4\n", "\n", "\n", "\n", "p2\n", "\n", "p2\n", "\n", "\n", "\n", "p0->p2\n", "\n", "\n", "1/4\n", "\n", "\n", "\n", "p3\n", "\n", "p3\n", "\n", "\n", "\n", "p0->p3\n", "\n", "\n", "1/4\n", "\n", "\n", "\n", "p4\n", "\n", "p4\n", "\n", "\n", "\n", "p0->p4\n", "\n", "\n", "1/4\n", "\n", "\n", "\n", "p1->p2\n", "\n", "\n", "1/3\n", "\n", "\n", "\n", "p1->p3\n", "\n", "\n", "1/3\n", "\n", "\n", "\n", "p1->p4\n", "\n", "\n", "1/3\n", "\n", "\n", "\n", "p2->p0\n", "\n", "\n", "1/2\n", "\n", "\n", "\n", "p2->p1\n", "\n", "\n", "1/2\n", "\n", "\n", "\n", "p3->p0\n", "\n", "\n", "1/2\n", "\n", "\n", "\n", "p3->p2\n", "\n", "\n", "1/2\n", "\n", "\n", "\n", "p4->p1\n", "\n", "\n", "1/3\n", "\n", "\n", "\n", "p4->p2\n", "\n", "\n", "1/3\n", "\n", "\n", "\n", "p4->p3\n", "\n", "\n", "1/3\n", "\n", "\n", "" ], "text/plain": [ "" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "import networkx as nx\n", "\n", "# Create a directed networkx graph\n", "G = nx.DiGraph()\n", "\n", "# Add outgoing web links (directed graph edges)\n", "G.add_edge(\"p0\", \"p1\", label=\"0.2\")\n", "G.add_edge(\"p0\", \"p2\")\n", "G.add_edge(\"p0\", \"p3\")\n", "G.add_edge(\"p0\", \"p4\")\n", "\n", "G.add_edge(\"p1\", \"p2\")\n", "G.add_edge(\"p1\", \"p3\")\n", "G.add_edge(\"p1\", \"p4\")\n", "\n", "G.add_edge(\"p2\", \"p0\")\n", "G.add_edge(\"p2\", \"p1\")\n", "\n", "G.add_edge(\"p3\", \"p0\")\n", "G.add_edge(\"p3\", \"p2\")\n", "\n", "G.add_edge(\"p4\", \"p1\")\n", "G.add_edge(\"p4\", \"p2\")\n", "G.add_edge(\"p4\", \"p3\")\n", "\n", "# Add weights to graph edges\n", "for node in G.nodes():\n", " # Get number of outgoing edges\n", " num_edges = len(G.edges(node))\n", "\n", " # For each outgoing edge, weight = 1/num_edges \n", " for edge in G.edges(node):\n", " G[edge[0]][edge[1]]['weight'] = 1/num_edges\n", " G[edge[0]][edge[1]]['label'] = \"{}/{}\".format(1, num_edges)\n", "\n", "# To plot graph, convert to a PyGraphviz graph for drawing\n", "Ag = nx.nx_agraph.to_agraph(G)\n", "Ag.layout(prog='dot')\n", "Ag.draw('web.svg')\n", "from IPython.display import SVG\n", "SVG('web.svg')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that edges have been given a weight which is $1/n$, where $n$ is the number of outgoing link/edges from a page/node. The is the fraction of page's importance/rank that it can give to the pages to which it links. \n", "\n", "If we denote the rank/importance of page $i$ by $x_{i}$, we can express the importance of each page:\n", "\n", "\\begin{align}\n", "x_{0} =& \\tfrac{1}{2} x_{2} + \\tfrac{1}{2} x_{3} \n", "\\\\\n", "x_{1} =& \\tfrac{1}{4}x_{0} + \\tfrac{1}{2}x_{2} + \\tfrac{1}{3}x_{4} \n", "\\\\\n", "x_{2} =& \\tfrac{1}{4}x_{0} + \\tfrac{1}{3}x_{1} + \\tfrac{1}{2}x_{3} + \\tfrac{1}{3}x_{4} \n", "\\\\\n", "x_{3} =& \\tfrac{1}{4}x_{0} + \\tfrac{1}{3}x_{1} + \\tfrac{1}{3}x_{4} \n", "\\\\\n", "x_{4} =& \\tfrac{1}{4}x_{0} + \\tfrac{1}{3}x_{1}\n", "\\end{align}\n", "\n", "We can express this as a system of equation:\n", "\n", "$$\n", "\\underbrace{\n", "\\begin{bmatrix}\n", "0 & 0 & \\tfrac{1}{2} & \\tfrac{1}{2} & 0\n", "\\\\\n", "\\tfrac{1}{4} & 0 & \\tfrac{1}{2} & 0 & \\tfrac{1}{3}\n", "\\\\\n", "\\tfrac{1}{4} & \\tfrac{1}{3} & 0 & \\tfrac{1}{2} & \\tfrac{1}{3}\n", "\\\\\n", "\\tfrac{1}{4} & \\tfrac{1}{3} & 0 & 0 & \\tfrac{1}{3}\n", "\\\\\n", "\\tfrac{1}{4} & \\tfrac{1}{3} & 0 & 0 & 0\n", "\\end{bmatrix}}_{\\boldsymbol{A}}\n", "\\begin{bmatrix}\n", "x_{0} \\\\ x_{1} \\\\ x_{2} \\\\ x_{3} \\\\ x_{4}\n", "\\end{bmatrix}\n", "=\n", "\\begin{bmatrix}\n", "x_{0} \\\\ x_{1} \\\\ x_{2} \\\\ x_{3} \\\\ x_{4}\n", "\\end{bmatrix}\n", "$$\n", "\n", "This is an eigenvalue problem with eigenvalue $\\lambda = 1$. To solve the problem we need to find the corresponding eigenvector." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "First, we can create the matrix $\\boldsymbol{A}$ directly from the graph (note that we added weights to the graph edges when building the graph):" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[0. 0. 0.5 0.5 0. ]\n", " [0.25 0. 0.5 0. 0.33333333]\n", " [0.25 0.33333333 0. 0.5 0.33333333]\n", " [0.25 0.33333333 0. 0. 0.33333333]\n", " [0.25 0.33333333 0. 0. 0. ]]\n" ] } ], "source": [ "# Get the matrix containing the graph weights. We provide node nodelist \n", "# so that the rows/columns are in the right order \n", "A = (nx.adjacency_matrix(G, nodelist=sorted(G.nodes())).T).toarray()\n", "print(A)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that the columns of $\\boldsymbol{A}$ sum to one. Such a matrix is known as a *stochastic matrix*. In our context, what is important is that the largest eigenvalue (in absolute terms) for such a matrix is one. Therefore, we are looking for the eigenvector associated with the largest eigenvalue." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Direct computation of the PageRank vector\n", "\n", "To find the solution, we can compute the eigenvectors of $\\boldsymbol{A}$, and pick the eigenvectors that corresponds to $\\lambda = 1$:" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The maxiumum eigenvalue is: (0.9999999999999991+0j)\n", "The PageRank vector (eigenvector) is: \n", " [0.215054+0.j 0.225806+0.j 0.258065+0.j 0.172043+0.j 0.129032+0.j]\n" ] } ], "source": [ "# Compute the eigenvalues and eigenvectors\n", "import numpy as np\n", "np.set_printoptions(precision=6)\n", "evalues, evectors = np.linalg.eig(A)\n", "\n", "# Print largest eigenvalue and corresponding eigenvector\n", "np.set_printoptions(precision=6)\n", "print(\"The maxiumum eigenvalue is: {}\".format(np.max(evalues)))\n", "evector = evectors[:, np.argmax(evalues)]\n", "evector /= np.linalg.norm(evector, 1)\n", "print(\"The PageRank vector (eigenvector) is: \\n {}\".format(evector))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The PageRank vector has been normalised using the $l_{1}$ norm such the the entries sum to one." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Approximation of the PageRank vector via power iteration\n", "\n", "The direct computation of all eigenvalues and eigenvectors is an expensive operation, and can only be reasonably performed for small matrices. However, we are interested only in the eigenvector corresponding to the largest eigenvalue. Recall that, in almost all cases, the repeated multiplication of a vector by a matrix $\\boldsymbol{A}$ will yield a vector that tends to the direction of the eigenvector corresponding to the largest eigenvalue. We can use this to approximate the PageRank vector using only matrix-vector multiplication:" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Estimated PageRank after 5 iterations: [0.216721 0.226431 0.256868 0.170802 0.129177]\n", "Estimated PageRank after 10 iterations: [0.215095 0.22581 0.258071 0.172022 0.129001]\n" ] } ], "source": [ "# Create random starting vector\n", "x0 = np.random.rand(A.shape[0])\n", "\n", "# Perform 5 iterations\n", "for i in range(5):\n", " x0 = A.dot(x0)\n", "x0 = x0/np.linalg.norm(x0 ,1)\n", "print(\"Estimated PageRank after 5 iterations: {}\".format(x0))\n", "\n", "# Perform another 5 iterations\n", "for i in range(5):\n", " x0 = A.dot(x0)\n", "x0 = x0/np.linalg.norm(x0 ,1)\n", "print(\"Estimated PageRank after 10 iterations: {}\".format(x0))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "It is clear that the approximate solution is very close to the exact solution with relatively few iterations. The major advantage of this method is that it is computationally inexpensive, which makes it tractable for very large matrices." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Multiple degree of freedom mass-spring system\n", "\n", "Eigenvalues and eigenvectors are used to understand the dynamic response of mass-spring systems.\n", "\n", "\n", "### System of equations\n", "\n", "The equations of motion for a multi degree-of-freedom mass-spring system are:\n", "\n", "$$\n", "\\boldsymbol{M} \\ddot{\\boldsymbol{x}} + \\boldsymbol{K} \\boldsymbol{x} = \\boldsymbol{0}\n", "$$\n", "\n", "where $\\boldsymbol{x}$ are the displacements of the masses. For the problem in the lecture with two masses (of equal mass) connected by springs of equal stiffness, the 'mass matrix' $\\boldsymbol{M}$ is \n", "\n", "$$\n", "\\boldsymbol{M} = \n", "\\begin{bmatrix}\n", "m_{1} & 0\\\\ 0 & m_{2} \n", "\\end{bmatrix}\n", "$$\n", "\n", "and the 'stiffness matrix' $\\boldsymbol{K}$ is \n", "\n", "$$\n", "\\boldsymbol{K} = \n", "\\begin{bmatrix}\n", "k_{1} + k_{2} & -k_{2} \\\\ -k_{2} & k_{2} \n", "\\end{bmatrix}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Eigenvectors and eigenvalues\n", "\n", "We now take a brief excursion into the properties of eiegenvectors related to the matrices $\\boldsymbol{K}$ and $\\boldsymbol{M}$. The purpose will soon be clear.\n", "\n", "Consider the eigenpairs for the 'generalised' eigenvalue problem\n", "\n", "$$\n", "\\boldsymbol{K} \\boldsymbol{u}_{i} = \\lambda_{i} \\boldsymbol{M} \\boldsymbol{u}_{i} \n", "$$\n", "\n", "where $\\lambda_{i}$ is an eigenvalue and $\\boldsymbol{u}_{i}$ is an eigenvector. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Orthogonality of the eigenvectors\n", "\n", "Multiplying both sides of the above equation by $\\boldsymbol{u}_{j}^{T}$:\n", "\n", "$$\n", "\\boldsymbol{u}_{j}^{T} \\boldsymbol{K} \\boldsymbol{u}_{i} = \\lambda_{i} \\boldsymbol{u}_{j}^{T} \\boldsymbol{u}_{i},\n", "$$\n", "\n", "Since the eigenvectors are orthogonal with respect to $\\boldsymbol{K}$ and $\\boldsymbol{M}$ (see lecture notes), and if they are normalised such that $\\boldsymbol{u}^{T}_{i} \\boldsymbol{M} \\boldsymbol{u}_{i} = 1$, it follows that\n", "\n", "$$\n", "\\boldsymbol{u}_{j}^{T} \\boldsymbol{M} \\boldsymbol{u}_{i} = \n", "\\begin{cases}\n", "1 & \\text{if} \\ i = j,\n", "\\\\\n", "0 & \\text{if} \\ i \\ne j.\n", "\\end{cases}\n", "$$\n", "\n", "From $\\boldsymbol{K} \\boldsymbol{u}_{i} = \\lambda_{i} \\boldsymbol{M} \\boldsymbol{u}_{i}$ and the above orthogonality property,\n", "it then also follows that\n", "\n", "$$\n", "\\boldsymbol{u}_{j}^{T} \\boldsymbol{K} \\boldsymbol{u}_{i}\n", "\\begin{cases}\n", "\\lambda_{i} & \\text{if} \\ i = j,\n", "\\\\\n", "0 & \\text{if} \\ i \\ne j.\n", "\\end{cases}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Decoupling the modes\n", "\n", "Since the eigenvectors are linearly independent, we can use them as a basis to express the displacement vector $\\boldsymbol{x}$\n", "\n", "$$\n", "\\boldsymbol{x} = \\sum_{i} \\alpha_{i} \\boldsymbol{u}_{i} = \\alpha_{1} \\boldsymbol{u}_{1} + \\alpha_{2} \\boldsymbol{u}_{2},\n", "$$\n", "\n", "where $\\alpha_{i}$ are unknowns. Taking two derivatives with respect to time:\n", "\n", "$$\n", "\\ddot{\\boldsymbol{x}} = \\sum_{i} \\ddot{\\alpha}_{i} \\boldsymbol{u}_{i} = \\ddot{\\alpha}_{1} \\boldsymbol{u}_{1} + \\ddot{\\alpha}_{2} \\boldsymbol{u}_{2}\n", "$$\n", "\n", "Inserting these expressions into $\\boldsymbol{M} \\ddot{\\boldsymbol{x}} + \\boldsymbol{K} \\boldsymbol{x} = \\boldsymbol{0}$,\n", "\n", "$$\n", "\\boldsymbol{M} \\sum_{i} \\ddot{\\alpha}_{i} \\boldsymbol{u}_{i} + \\boldsymbol{K} \\sum_{i} \\alpha_{i} \\boldsymbol{u}_{i} = \\boldsymbol{0}\n", "$$\n", "\n", "Multiplying both sides of the above equation by $\\boldsymbol{u}^{T}_{j}$ \n", "\n", "$$\n", "\\boldsymbol{u}^{T}_{j} \\boldsymbol{M} \\sum_{i} \\ddot{\\alpha}_{i} \\boldsymbol{u}_{i} + \\boldsymbol{u}^{T}_{j} \\boldsymbol{K} \\sum_{i} \\alpha_{i} \\boldsymbol{u}_{i} = \\boldsymbol{0}\n", "$$\n", "\n", "From the orthogonality properties of the eigenvectors, for the case $i \\ne j$ the above becomes the trivial '$0 =0$', and when $i = j$ we get the ordinary differential equation\n", "\n", "$$\n", " \\ddot{\\alpha}_{i} + \\lambda_{i} \\alpha_{i} = 0\n", "$$\n", "\n", "for every $i$. This is a second-order, constant coefficient equation which is easy to solve:\n", "\n", "$$\n", " \\alpha_{i} = A_{i}\\sin(\\sqrt{\\lambda_{i}}t) + B_{i}\\cos(\\sqrt{\\lambda_{i}}t)\n", "$$\n", "\n", "Hence the displacement vector is \n", "\n", "$$\n", "\\boldsymbol{x} = \\sum_{i} (A_{i}\\sin(\\sqrt{\\lambda_{i}}t) + B_{i}\\cos(\\sqrt{\\lambda_{i}}t)) \\boldsymbol{u}_{i}\n", "$$\n", "\n", "where $\\sqrt{\\lambda_{i}}$ is the $i$th natural frequency, and $\\boldsymbol{u}_{i}$ is known as \n", "the $i$th vibration mode. The constants $A_{i}$ and $B_{i}$ are determined from the initial conditions." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Solution for the $2 \\times 2$ system\n", "\n", "To get the complete solution describing the motion of the masses we need to compute the eigenvalues and eigenvectors of of $\\boldsymbol{K} \\boldsymbol{x} = \\lambda \\boldsymbol{M} \\boldsymbol{x}$. We first create the matrices $\\boldsymbol{K}$ and $\\boldsymbol{M}$: " ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[ 2. -1.]\n", " [-1. 1.]]\n" ] } ], "source": [ "import numpy as np\n", "k1, k2 = 1.0, 1.0\n", "K = np.array([[k1 + k2, -k2], [-k2, k2]] )\n", "print(K)" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[1. 0.]\n", " [0. 1.]]\n" ] } ], "source": [ "m1, m2 = 1.0, 1.0\n", "M = np.array([[m1, 0], [0, m2]] )\n", "print(M)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We now compute the eigenvalues and eigenvectors of the generalised problem:" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[0.381966 2.618034]\n", "[[-0.525731 -0.850651]\n", " [-0.850651 0.525731]]\n" ] } ], "source": [ "from scipy.linalg import eigh # Use eigh, which is specialised for symmetric matrices\n", "\n", "evalues, evectors = eigh(K, M)\n", "print(evalues)\n", "print(evectors)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The motion of the masses is therefore\n", "\n", "$$\n", "\\boldsymbol{x} = \n", "(A_{1}\\sin(\\sqrt{0.381966} t) + B_{1}\\cos(\\sqrt{0.381966} t)) \\begin{bmatrix}-0.525731 \\\\ -0.850651\\end{bmatrix}\n", "+\n", "(A_{2}\\sin(\\sqrt{2.618034} t) + B_{2}\\cos(\\sqrt{2.618034} t)) \\begin{bmatrix} -0.850651 \\\\ 0.525731\\end{bmatrix}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Exercise\n", "\n", "Give the two springs different stiffness and examine the influence on the natural frequencies." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Motion for given initial conditions\n", "\n", "If at $t = 0$ the conditions are $\\boldsymbol{x} = \\begin{bmatrix} 1 & 1 \\end{bmatrix}^{T}$ and $\\dot{\\boldsymbol{x}} = \\begin{bmatrix} 0 & 0 \\end{bmatrix}^{T}$, we have \n", "\n", "$$\n", "\\begin{bmatrix}\n", "-0.525731 & -0.850651 \n", "\\\\\n", "-0.850651 & 0.525731 \n", "\\end{bmatrix}\n", "\\begin{bmatrix}\n", "B_{1} \\\\ B_{2}\n", "\\end{bmatrix}\n", "=\n", "\\begin{bmatrix}\n", "1 \\\\ 1\n", "\\end{bmatrix}\n", "$$\n", "\n", "We can compute $B_{1}$ and $B_{2}$ using NumPy:" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[-1.376382 -0.32492 ]\n" ] } ], "source": [ "B = np.linalg.inv(evectors).dot(np.array([1, 1]))\n", "print(B)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Animating the motion of the masses\n", "\n", "For this case, we can use Matplotlib to crate an animation of the motion of the masses. We evaluate the solution to the ODE for each mass for a range of time steps." ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", "\n", "
\n", " \n", "
\n", " \n", "
\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
\n", " Once \n", " Loop \n", " Reflect \n", "
\n", "
\n", "\n", "\n", "\n" ], "text/plain": [ "" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "%matplotlib inline\n", "import matplotlib.pyplot as plt\n", "import matplotlib.animation as animation\n", "from IPython.display import HTML\n", "\n", "fig = plt.figure()\n", "ax = fig.add_subplot(111, autoscale_on=False, xlim=(0, 6), ylim=(-2, 2))\n", "ax.grid()\n", "line, = ax.plot([], [], 'o-', markersize=10, lw=2)\n", "time_template = 'time = %.1fs'\n", "time_text = ax.text(0.05, 0.9, '', transform=ax.transAxes)\n", "\n", "def x(t):\n", " alpha0 = B[0]*np.cos(np.sqrt(evalues[0])*t)\n", " alpha1 = B[1]*np.cos(np.sqrt(evalues[1])*t)\n", " return evectors[:,0]*alpha0[:,np.newaxis] + evectors[:,1]*alpha1[:,np.newaxis] \n", "\n", "t = np.linspace(0, 8*np.pi, 100, retstep=True)\n", "dt = t[1]\n", "disp = x(t[0]) \n", "\n", "def animate(i):\n", " zero = [0, 0, 0]\n", " x0 = [0, disp[i, 0] + 2.0, disp[i, 1] + 4.0]\n", " x1 = [0, 0, 0]\n", " \n", " line.set_data(x0, zero)\n", " time_text.set_text(time_template%(i*dt))\n", " return line, time_text\n", "\n", "def init():\n", " line.set_data([], [])\n", " time_text.set_text('')\n", " return line, time_text\n", "\n", "anim = animation.FuncAnimation(fig, animate, np.arange(0, len(disp)),\n", " interval=100, blit=False, init_func=init)\n", " \n", "# Call animation function to display the animation\n", "plt.close(anim._fig)\n", "HTML(anim.to_jshtml())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Exercise\n", "\n", "Try changing the spring stiffnesses $k_{1}$ and $k_{2}$, and the masses $m_{1}$ and $m_{2}$ to investigate what effect these have on the vibration response." ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.6.7" } }, "nbformat": 4, "nbformat_minor": 1 }