{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# The Hanging Chain: Part 2\n", "\n", "For mathematical details, see https://scipython.com/blog/the-hanging-chain-part-2/\n", "\n", "Change the definition of the function `f(u)` to change the initial conditions of the chain." ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import numpy as np\n", "from scipy.special import j0, j1, jn_zeros\n", "from scipy.integrate import quad\n", "import matplotlib.pyplot as plt\n", "from matplotlib import animation\n", "from IPython.display import HTML" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Acceleration due to gravity, m.s-2\n", "g = 9.81\n", "# Chain length, m\n", "L = 1\n", "# Vertical axis (m)\n", "N = 201\n", "z = np.linspace(0, L, N)\n", "# Scaled vertical axis\n", "u = 2 * np.sqrt(z/g)" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": false, "scrolled": true }, "outputs": [], "source": [ "# The initial state of the chain, x(z,0): linear from (0,d) to (c,0) then zero\n", "# from (c,0) to (L,0).\n", "d, c = 0.05, L/4\n", "p, q = -d/c, d\n", "def f(u):\n", " z = u**2 * g / 4\n", " if np.isscalar(z):\n", " if z > c:\n", " return 0\n", " return p*z + q\n", " h = p*z + q\n", " h[z>c] = 0\n", " return h\n", "\n", "nmax = 10\n", "# jn_zeros calculates the first n zeros of the zero-th order Bessel function of the first kind\n", "w = jn_zeros(0, nmax)\n", "b = 2 * np.sqrt(L/g)\n", "def func(u, n):\n", " return u * f(u) * j0(w[n-1] / b * u)\n", "A = []\n", "for n in range(1, nmax+1):\n", " An = quad(func, 0, b, args=(n,))[0] * 2 / (b * j1(w[n-1]))**2\n", " A.append(An)" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhAAAAFkCAYAAABxWwLDAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzt3XmcV9Wd5//XKZZiRzYBjUHBCMQNQYGCQom4xC2bJoYk\nY8akY6fNxBl60na659Gd7vQvcexMdNK/STr2dHc0nYRO0t0maEdtVBSqWJRFFHdREpVFZV+KKoo6\n88f91koV8C2ouqeo1/Px+FJV93vu9/u5h4J617n33BNijEiSJBWjJO8CJElS12OAkCRJRTNASJKk\nohkgJElS0QwQkiSpaAYISZJUNAOEJEkqmgFCkiQVzQAhSZKKZoCQJElFKzpAhBBmhRAWhBDeDiHU\nhRA+chT7zA4hrAoh7A8hvBJC+Hz7ypUkSSlozwhEf+AZ4FbgiAtphBBOBx4EHgPOB74H/H0I4fJ2\nvLckSUpAOJbFtEIIdcDHYowLDtPmTuCqGON5TbbNBwbHGK9u95tLkqTcdMY1ENOBR1tsewQo64T3\nliRJHaBnJ7zHKGBLi21bgEEhhNIYY3XLHUIIw4ArgQ3A/g6vUJKkE0cf4HTgkRjj1o56k84IEO1x\nJfDTvIuQJKkL+yzws4568c4IEJuBkS22jQR2tTb6ULAB4Cc/+QkTJ07swNLStXTpUr761V3AP3PL\nLWX8/u//ft4l5W7evHncfffdeZeRO/shYz80si8y9kPmxRdf5HOf+xwUfpZ2lM4IEMuAq1psu6Kw\nvS37ASZOnMjkyZM7qq6kvfPOO2STXO7nlFNO6bb90NTgwYPtB+yHevZDI/siYz8cokMvAWjPfSD6\nhxDODyFMKmwaW/j6tMLzd4QQ7muyyw8Lbe4MIYwPIdwK3ADcdczVn/BOAg7kXYQkSYdozyyMC4E1\nwCqyX5G/C6wG/rLw/CjgtPrGMcYNwDXAZWT3j5gHfDHG2HJmhg4xGAOEJClFRZ/CiDE+yWGCR4zx\n5la2LQamFPtecgRCkpQm18JI2knA2LyLSMbcuXPzLiEJ9kPGfmhkX2Tsh85lgEjUgQMB6MehE1i6\nL/9zyNgPGfuhkX2RsR86lwEiUXv31p9d2pFrHZIktcYAkai9e3sVPtuZax2SJLXGAJEoRyAkSSkz\nQCSqcQTCACFJSo8BIlF79tSPQHgKQ5KUHgNEovbvrw8Qu3OtQ5Kk1hggErV/fw+gGjiYdymSJB3C\nAJGoLEDszbsMSZJaZYBIVBYg9uVdhiRJrTJAJMoRCElSygwQiTJASJJSZoBIVHW1AUKSlC4DRKIc\ngZAkpcwAkShHICRJKTNAJKqmpgTYn3cZkiS1ygCRqIMHA3Ag7zIkSWqVASJRtbUlQG3eZUiS1CoD\nRKLq6gIGCElSqgwQicpOYRggJElpMkAkymsgJEkpM0AkymsgJEkpM0AkKhuBcClvSVKaDBCSJKlo\nBohE9egRgR55lyFJUqsMEIkqKTFASJLSZYBIVDYC0TPvMiRJapUBIlGOQEiSUmaASJTXQEiSUmaA\nSJQjEJKklBkgEuU1EJKklBkgEtW790Ggb95lSJLUKgNEovr0OQj0z7sMSZJaZYBIVGlpHdAv7zIk\nSWqVASJRffrU4giEJClVBohEZSMQBghJUpoMEInKroHwFIYkKU0GiESVlnoRpSQpXQaIRPXtWwsM\nyLsMSZJaZYBI1IABB8gCRK+8S5Ek6RAGiEQNHHig8NnQXOuQJKk1BohEDRpkgJAkpcsAkajGEYhh\nudYhSVJrDBCJ8hSGJCllBohEZRdRggFCkpQiA0SievaMwE48hSFJSpEBImnbcARCkpQiA0TSDBCS\npDQZIJK2FU9hSJJSZIBImiMQkqQ0tStAhBC+EkJ4I4RQFUJYHkK46AjtPxtCeCaEsDeEsDGE8A8h\nBH8yHpEBQpKUpqIDRAjhRuC7wDeAC4C1wCMhhOFttJ8J3Af8X+CDwA3AVODv2llzN+IpDElSmtoz\nAjEPuCfG+OMY40vAl4F9wBfaaD8deCPG+P0Y429jjEuBe8hChA7LEQhJUpqKChAhhF7AFOCx+m0x\nxgg8CpS1sdsy4LQQwlWF1xgJfBL49/YU3L1sAwZQW9sj70IkSWqm2BGI4UAPYEuL7VuAUa3tUBhx\n+Bzw8xBCDbAJ2A78lyLfuxvaCkBVVd+c65AkqbmeHf0GIYQPAt8D/gL4D2A08L/ITmP83uH2nTdv\nHoMHD262be7cucydO7dDak3PNsAAIUlq3fz585k/f36zbTt37uyU9y42QLwHHARGttg+Etjcxj5f\nBypjjHcVvl4XQrgVWBJC+B8xxpajGQ3uvvtuJk+eXGSJJ5IsQOzbZ4CQJB2qtV+qV69ezZQpUzr8\nvYs6hRFjPACsAubUbwshhMLXS9vYrR9Q22JbHRCBUMz7dz9ZgNi/3wAhSUpLe2Zh3AV8KYRwUwhh\nAvBDspBwL0AI4Y4Qwn1N2j8AXB9C+HII4YzCtM7vAStijG2NWgjILhXxFIYkKT1FXwMRY/xF4Z4P\n3yQ7dfEMcGWM8d1Ck1HAaU3a3xdCGAB8hezahx1kszi+foy1dwO1wB6qqvrkXYgkSc206yLKGOMP\ngB+08dzNrWz7PvD99ryXtnkKQ5KUHNfCSN52RyAkSckxQCRvO/v3GyAkSWkxQCRvuxdRSpKSY4BI\n3jYDhCQpOQaI5HkKQ5KUHgNE8nYYICRJyTFAJG8fBw70yrsISZKaMUAkLwsQMeZdhyRJjQwQydtH\njCXU1ORdhyRJjQwQydsLwL59OZchSVITBojkZcnBACFJSokBInkGCElSegwQyTNASJLSY4BIXnb1\nZHV1zmVIktSEASJ5dQBO45QkJcUAkbwsORggJEkpMUAkLxuBqKvLuQxJkpowQCTPEQhJUnoMEMlz\nBEKSlB4DRPJC9mfIuQxJkpowQCSvNPuzNOcyJElqwgCRPAOEJCk9BojkGSAkSekxQCTPACFJSo8B\nInkGCElSegwQyesLQJ8+OZchSVITBojkDQZg0KCcy5AkqQkDRPKGUFq6nx498q5DkqRGBojknUSf\nPvvzLkKSpGYMEMkzQEiS0mOASN4QA4QkKTkGiOSdRN++VXkXIUlSMwaI5HkKQ5KUHgNE8obTt68B\nQpKUFgNE8kYzYMDuvIuQJKkZA0TS+gODGDjQACFJSosBImmjARg4cE/OdUiS1JwBImlZgPAUhiQp\nNQaIpDkCIUlKkwEiaaOBfZSWVuddiCRJzRggkjYa2EQIedchSVJzBoikZQFCkqTUGCCSdgqwMe8i\nJEk6hAEiaY5ASJLSZIBImgFCkpQmA0SiampKgKEYICRJKTJAJGr79tLCZwYISVJ6DBCJ2ratPkBs\nzrUOSZJaY4BI1LZtvQufOQIhSUqPASJR2SmMWmBr3qVIknQIA0SisgCxBYh5lyJJ0iHaFSBCCF8J\nIbwRQqgKISwPIVx0hPa9QwjfCiFsCCHsDyG8HkL4z+2quJvYsaM3Xv8gSUpVz2J3CCHcCHwXuAV4\nCpgHPBJCOCvG+F4bu/0SGAHcDKwnu8GBox+HkY1AGCAkSWkqOkCQBYZ7Yow/BgghfBm4BvgC8Nct\nG4cQPgzMAsbGGHcUNv+ufeV2H9ksDAOEJClNRY0ChBB6AVOAx+q3xRgj8ChQ1sZu1wErgT8OIbwV\nQng5hPCdEEKfdtbcLTReAyFJUnqKHYEYDvTg0J9sW4DxbewzlmwEYj/wscJr/C3ZbRa/WOT7dxu7\nd/cC2jojJElSvtpzCqNYJUAd8JkY4x6AEMIfAr8MIdwaY6xua8d58+YxePDgZtvmzp3L3LlzO7Le\n3NXWQlVVT2Bn3qVIkhI2f/585s+f32zbzp2d87Oj2ADxHnAQGNli+0jaPmG/CXi7PjwUvAgE4H1k\nF1W26u6772by5MlFltj17dpV/9mOwzWTJHVzrf1SvXr1aqZMmdLh713UNRAxxgPAKmBO/bYQQih8\nvbSN3SqBU0II/ZpsG082KvFWUdV2EzsacoMjEJKkNLVnKuVdwJdCCDeFECYAPwT6AfcChBDuCCHc\n16T9z8hup/ijEMLEEMLFZLM1/uFwpy+6s8bRJ0cgJElpKvoaiBjjL0IIw4Fvkp26eAa4Msb4bqHJ\nKOC0Ju33hhAuB/5/4GmyMPFz4M+OsfYTVuMIhAFCkpSmdl1EGWP8AfCDNp67uZVtrwBXtue9uqO9\nexs+y7MMSZLa5N0gE1RT0/BZnmVIktQmA0SCDhyo/8wAIUlKkwEiQY5ASJJSZ4BIUOMIxIHDNZMk\nKTcGiASVNPythDzLkCSpTQaIBPVsmBvTK88yJElqkwEiQb0ackNnLFUiSVLxDBAJcgRCkpQ6A0SC\neveu/6w0zzIkSWqTASJBgwY1fJZnGZIktckAkaCTTmr4LM8yJElqkwEiQYMHN3yWZxmSJLXJAJGg\nxgDhCIQkKU0GiAQNHgw9etQBI/IuRZKkVhkgElRSAkOHVgOn5l2KJEmtMkAkavjw/RggJEmpMkAk\natgwRyAkSekyQCRq2DBHICRJ6TJAJMpTGJKklBkgEjVixH6y+0A4lVOSlB4DRKJOPXVv4bOzcq1D\nkqTWGCASdeqp+wqfGSAkSekxQCSqb9+DwFvA+LxLkSTpEAaIpL2MAUKSlCIDRNJewVMYkqQUGSCS\n9jLwAerq8q5DkqTmDBBJewXox65dLustSUqLASJpLwOwdeuwnOuQJKk5A0TSNgDVvPeeAUKSlBYD\nRNLqgPVs3To870IkSWrGAJG8lx2BkCQlxwCRvJe9BkKSlBwDRPJeZufOk9i378gtJUnqLAaI5GUz\nMV59NecyJElqwgCRvCxAvPxyzmVIktSEASJ52+jbdx8vvZR3HZIkNTJAdAHDhm1l/fq8q5AkqZEB\nogsYPHgnb72VdxWSJDUyQHQBgwbtNkBIkpJigOgCBg3ayZtvQox5VyJJUsYA0QUMGrSLqirYvj3v\nSiRJyhgguoD+/bO7SG3dmnMhkiQVGCC6gN69qwHYtSvnQiRJKjBAdAGlpVmA2L0750IkSSowQHQB\npaU1gCMQkqR0GCC6gPpTGI5ASJJSYYCQJElFM0B0ATEGAHr0yLkQSZIKDBBdQIzZX1OJf1uSpET4\nI6kLqK4uBWDgwJwLkSSpwADRBezf3weAIUNyLkSSpIJ2BYgQwldCCG+EEKpCCMtDCBcd5X4zQwgH\nQgir2/O+3VVVlQFCkpSWogNECOFG4LvAN4ALgLXAIyGE4UfYbzBwH/BoO+rs1nbtGgTAKafkXIgk\nSQXtGYGYB9wTY/xxjPEl4MvAPuALR9jvh8BPgeXteM9ubceOkxg2zGsgJEnpKCpAhBB6AVOAx+q3\nxRgj2ahC2WH2uxk4A/jL9pXZve3YcRJjxuRdhSRJjYodgRgO9AC2tNi+BRjV2g4hhA8A3wY+G2Os\nK7pC8e67I5gwIe8qJElq1LMjXzyEUEJ22uIbMcb19ZuPdv958+YxePDgZtvmzp3L3Llzj1+RyQts\n2TKS88/Puw5JUmrmz5/P/Pnzm23buXNnp7x3sQHiPeAgMLLF9pHA5lbaDwQuBCaFEL5f2FYChBBC\nDXBFjPGJtt7s7rvvZvLkyUWWeKI5g5qaUgOEJOkQrf1SvXr1aqZMmdLh713UKYwY4wFgFTCnflsI\nIRS+XtrKLruAc4BJwPmFxw+Blwqfr2hX1d1KdmlJJ3wvSJJ01NpzCuMu4N4QwirgKbJZGf2AewFC\nCHcAp8QYP1+4wPKFpjuHEN4B9scYXzyWwruPD3HyyVsYPrzloI8kSfkpOkDEGH9RuOfDN8lOXTwD\nXBljfLfQZBRw2vErsbubzemnb+DQs0aSJOWnXXeijDH+IMZ4eoyxb4yxLMa4sslzN8cYLz3Mvn8Z\nY+zuFzYcpdOAcYUAIUlSOlwLI2lZDhsz5rc51yFJUnMGiKR9FFhG//5VeRciSVIzBohE7d9fAlwJ\n3J93KZIkHcIAkag1a4aTTW75Vd6lSJJ0CANEopYuHUk2A/bVvEuRJOkQBogE1dbCihUj8PSFJClV\nBogELVkCu3f3xtMXkqRUGSAS9KtfwbBh+4GVR2wrSVIeDBCJiRF+/WsoK2u5YrokSekwQCRm3Tr4\n7W9h+vR3j9xYkqScGCAS88ADMGAAnHvu1rxLkSSpTQaIxDz4IFxxBfTuHfMuRZKkNhkgEvLuu7B8\nOVx7bd6VSJJ0eAaIhCxenF1EecUVeVciSdLhGSASsnQpjBkDp56adyWSJB2eASIhy5ZBWVneVUiS\ndGQGiETECM89B5Mn512JJElHZoBIxPbtsGcPnH563pVIknRkBohE/Pa32ccxY/KtQ5Kko2GASMQ7\n72QfR47Mtw5Jko6GASIRtbXZx1698q1DkqSjYYBIRP/+2ce9e/OtQ5Kko2GASMTJJ2cfN23Ktw5J\nko6GASIRZ54JPXrA88/nXYkkSUdmgEhE797ZPSCefDLvSiRJOjIDREI+/GF4+GGoqsq7EkmSDs8A\nkZCbboKdO+GXv8y7EkmSDs8AkZAzz4Srr4ZvfQsOHgx5lyNJUpsMEIn51rfglVdgwYL3512KJElt\nMkAkZtIkuO02uPfes4AJeZcjSVKrDBAJuuMOGDVqH/Br4KS8y5Ek6RAGiAT16wff+MYaYBiwgOpq\n728tSUqLASJRp5yyD7gGOJ+f/ewz7NqVd0WSJDUyQCRtBfBhNm8eRXk5/O53edcjSVLGAJG8ZXzx\ni//A7t1w0UXw6KN51yNJkgGiSzj55PdYsQLOPx+uuAL+9E+hpibvqiRJ3ZkBoos4+eTsNtff/jZ8\n5zswZQqsWJF3VZKk7soA0YWUlMDXvw4rV0JpKZSVwa23wnvv5V2ZJKm7MUB0QeefD8uXw113wc9+\nlt0C+667oLo678okSd2FAaIL2L17NzHGZtt69oT/9t/g1Vdh7lz4oz/KgsQ993h9hCSp4xkguoDv\nfve7jBgxgo9+9KN85zvfYdmyZVQXhhtGjIC//Vt4/nkoL4c/+AMYPz7b5rLgkqSOElr+ZpuCEMJk\nYNWqVauYPHly3uXkYuPGjYwdO7YhKLRUWlrK1KlTKS8vp7y8nLKyMoYMGcK6dfBXfwX/8i8wdCh8\n9avZdRLDh3fyAUiScrF69WqmTJkCMCXGuLqj3scAkbD169fzr//6r1RUVFBZWcm2bdsO2/6cc86h\nvLycmTNn8v73z+bnPz+VH/0oUFcHn/40fOUr2b0kJEknLgOEAaKZuro6Xn75ZSoqKhoCxfr16w+7\nz6mnnspFF11FXd1/ZuXKC9m4sZSLLoJbboFPfQoGDeqk4iVJncYAYYA4ok2bNlFZWdkQKNasWcPB\ngwfbaF1C377X06/fH7Ft24WUltZx/fV13HJLL2bNghA6tXRJUgcxQBggirZnzx5WrFjRECiWLVvG\nnj17Wmn5PuAm4AvAOAYN2spll73LbbcN45JLRnRu0ZKk48oAYYA4ZrW1tTz77LMNgWLJkiVs2rSp\nSYsAzAI+C3wSGELv3i9wzjnr+MQn4Prrz2P8+PEEhyckqcswQBggjrsYIxs2bGg47VFRUcHzzz9f\neLY3cBVZmLgG6Aeso2/fh5k69S2uvvpUZs0qZ/LkyZSWluZ1CJKkIzBAGCA6xbZt21i2bFlDoHj6\n6aepri4BrgQ+AVwLDAHeAO6nV68HmTatjlmzZjBz5kxmzJjBkCFD8jwESVITBggDRC6qq6tZtWpV\nQ6CoqFjB9u3nkYWJjwGjgS3Ar4EFwCLOOWdsw/TR8vJyxowZ42kPScqJAcIAkYSm00eXLKnk8cf3\n8vbbU4GPA2cC1cAS4GHgIeAFTj311GaB4rzzzqNHjx45HoUkdR9JB4gQwleArwGjgLXAV2OMT7fR\n9uPAHwCTgFLgeeAvYoz/cZjXN0AkbNOmTVRUVPLggy+zaFEpb755NjAb6Au8SRYmHgYeA3YycOBA\nysrKGgLFtGnT6N+/f34HIEknsGQDRAjhRuA+4BbgKWAe2SX8Z8UYD1lYOoRwN/A2sAjYQTZ38GvA\n1Bjj2jbewwDRhezZs4fFi59m/vy3WLy4P2++eTYxjgdqgWVkYWIhsBo4SI8ePZg8eXJDoJg5cyaj\nRo3K8xAk6YSRcoBYDqyIMf7XwteB7NfOv4kx/vVRvsY64J9jjP9fG88bILqw2tpaHn74Jf7pn96h\nsnIAGzd+kBgHADuBxWRZ8nHgWSD7/hs3blzDuh4zZ85kwoQJXkchSe3QWQGiZzGNQwi9gCnAt+u3\nxRhjCOFRoOwoXyMAA4HDL+ygLqtnz55ce+05XHtt9nVNTWTBgo38/OfvsHTpB9i06XJi7ANsBZ4E\nHmf9+kWsX38f9913HwDDhg1rNkIxZcoUp49KUkKKChDAcKAH2WX4TW0Bxh/la/wR0B/4RZHvrS6q\nd+/ADTecwg03nALA/v2wcOEu5s/fypIlF/D229cRYy9gM/AEsIitW59gwYIFLFiwAGi++qjTRyUp\nf0WdwgghjCa7nqEsxriiyfY7gYtjjIcdhQghfAa4B/hIjHHRYdpNBlZdfPHFDB48uNlzc+fOZe7c\nuUdds9K3bx8sWlTDz3/+Dk8+GXjzzVHE2AN4B6ggm+VRAawBGtf6aLr6qNNHJXVH8+fPZ/78+c22\n7dy5k8WLF0NK10AUTmHsA66PMS5osv1eYHCM8eOH2ffTwN8DN8QYHz7C+3gNRDe2axdUVtaxYME2\nHn+8hvXrh3PwYG9gD7CcxlCxAtjbsJ/TRyWp611E+Tuyiyi/08Y+c8nCw40xxgeP4j0MEGpQUwOr\nVsFvfrOLhx7aw7p1g6mu7k82y2M1jYGiEni3Yb+BAwcyffr0hosznT4qqTtIOUB8CrgX+DKN0zhv\nACbEGN8NIdwBnBJj/Hyh/WcK7W8D7m/yUlUxxl1tvIcBQm2qq4MXX4RHH93PAw9sZ+XKvuzceVLh\n2fVkU0eXFj4+R/1pjx49enDBBRc0m+3h9FFJJ5pkAwRACOFW4HZgJPAM2Y2kVhae+xEwJsZ4aeHr\nRcDFrbzMfTHGL7Tx+gYIFeV3v4MlSw7y4INbWbq0jjffHE6MPclOezxNFibqH1sb9nP6qKQTTdIB\noqMZIHSsqqpg5crIb36znYUL9/DCC4Opqqq/IPdVmo9SrAPqAKePSur6DBAGCB1HMcKGDbBw4R4e\neOA9Vq7syebNo8hmMu8BVpKdkVtReLwNOH1UUtdjgDBAqIPt2wdLl9Zw//0bWby4mtdeG8r+/SMK\nz26kMUw8RRYwdgPZ9NH6UQqnj0pKjQHCAKEcvPVWHfff/za/+c1WnnmmF1u2jCnchrsOeIHmoxTr\ngIOceuqpzQKF00cl5ckAYYBQAg4ehCVL3uWXv9zA4sXZPSmqqsaRnfrYB6yicZRiBfA7p49KylWS\na2FI3U2PHjB79ghmzx7RsO2dd/bws5+9xEMPbWPt2t68884nifFrhWc3s3v3UyxcuIKFC5cAd9Gj\nxx6nj0o64TgCIR2j2tpaFi16gV/+8rdUVBxg/fph1NScD9Tfm+IlshGK+sdaxo07zemjkjqEpzAM\nEOqiYoy8/voG/u3fnuOhh7aydm1vtm07E5gElAI1ZLdPaQwVQ4dupbx8htNHJR0zA4QBQieQbdu2\nsXjxCn71q9dZsqSaDRtGUlc3BZhQaLGD7IZXWaDo3Xst06a9v+HiTKePSjpaBggDhE5g1dXVrFq1\nioULn+bhh99j7dpSqqrOAaYBowut3qTpKMXEiVVcfPEkp49KOiwDhAFC3UhdXR0vv/wyS5ZU8B//\n8TwVFTVs2fJ+YCpwETCQbCrpi9TP+hgxYgOXXDKUiy8uc/qopAYGCAOEurlNmzZRWVnJ4sWVPP74\n27zwwkBivJAsVJwH9AKqgDXAU/Tp8xxTphzk8svHMmuW00el7soAYYCQmtmzZw8rVqygoqKCxYuf\nZtmyqsJpj6mFxwcKLbcCTxHCSsaN28qllw7kiismOX1U6iYMEAYI6bBqa2t59tlnqaiooLKykief\nXMeWLe+jMVBMA04utH4DeIphw15j6tQ6rrvuVGbPLnP6qHQCMkAYIKSixBjZsGEDlZWVVFRUsGRJ\nBS+8sIfmgeJCoC9QDayhT59nOPvs3Vx++QCuvfZ8LrzQ6aNSV2eAMEBIx2zbtm0sW7aMiooKKioq\neOqpNdTUTACmA2WFj+MKrTdSUrKCMWM2MWNGCZ/4xBg+9KHpTh+VuhgDhAFCOu7qp4/WB4rKykq2\nbetJNjpRHyimAv2BA8AzDB36KpMmVXHNNcP4+McncfrpTh+VUmaAMEBIHa5++mjTQLF+/QbgHBoD\nRRlwVmGPLfTp8wwf+MBWZs/uw6c/PY5p085x+qiUEAOEAULKRf300fpAsWbNGg4ePInmoxTTyO5N\nUUtJybO8730bmD79IDfcMJqrr77A6aNSjgwQBggpCU2nj1ZWVrJs2TL27NkHfBCYAcwEyoGxhT1e\nY9iwl7jggj1cd91QPvnJ8xg92umjUmcxQBggpCS1nD66ZMkSNm3aBJxCY5iYSbZ4WA/gPfr1W8PE\niVu5/PK+fPrTZ3HeeU4flTqKAcIAIXUJLaePVlRU8PzzzwMDyE53lBce08kuzqyiZ881nHHG25SX\nw403vp/Zsyc5fVQ6TgwQBgipy2o5ffTpp5+muvog2ahEeZPHSKCOENYxatR6LrqomuuvH851101x\n+qjUTgYIA4R0wmh9+ug24EyaB4rxhT1eYciQZ7nggp189KMn8ZGPTHH1UekoGSAMENIJq/Xpo+vJ\nbr19MTC78Di7sMer9Ov3FOeeu42rr+7PddddwLnnnkvPnj3zOQApYQYIA4TUrbQ+fXQozQPFOYXW\nr9GzZyXjx2/myitLueaa81x9VCowQBggpG6t9emjfWgeKM4ttF5PCE8yZsxvufTSEq6++hxXH1W3\nZYAwQEhqovXpozU0DxTnFVq/DjzBiBEvcvHFdVx55QTKy8tdfVTdggHCACHpMFqfProZmEVjoDi/\n0PoN4En6919JWVk1l1/+AcrLy5kyxdVHdeIxQBggJBXp0NVH11NTM5UsTFxCFihKgA3AInr2rGTS\npO1cdtkrCHX/AAAOi0lEQVRZlJeXM2PGDKePqsszQBggJB2jltNHlyx5nh07zqVxhGJSoeXrwCLg\nCc46ayMf+lA2QlFeXu70UXU5BggDhKTjrOX00cWLn2fDhtPIwsSHaLyG4jXqA8XIkS9xySVnNgQK\np48qdQYIA4SkTtB0+ugTT6zjueeGUFd3MVmgqJ82+grwBLCIfv2eZubMsQ2BwumjSo0BwgAhKQdN\np48uWvQ8K1aUsn//dJrf2Ool6gNFSUkFkyefQnl5OTNnzmTmzJmMHj06p+olA4QBQlISmk4ffeyx\ndSxeHNixYxJZoJhYaPUC9YECnmTcuEENgcLpo+psBggDhKQENZ0++sgjz/LEE5G33jqTLFDUr+Xx\nArAEqACWMHToHmbOnNFw2sPpo+pIBggDhKQuon766EMPPcvChQd47bXR1NXNoPGUx5tkgSILFb17\nv8bUqRc2BAqnj+p4MkAYICR1UfXTRx95ZCUPPbSL554bzP79FwFTgF7ANqCSxlCxirPPPqshUDh9\nVMfCAGGAkHSCqJ8++thjy3jggS08/XRftm8/GygDBgBVwAoaA8UyTjllULNA4fRRHS0DhAFC0gls\n06ZNLF68lF/9agMVFYG33joDmEm2pPlB4BmaBooBA3ZTVlbWcHHmtGnTGDBgQI5HoFQZIAwQkrqR\nPXv2sHz5Cn796xd57LEaXn11FLW104GxhRZvAssLjxWUlDzD5MkTnT6qQxggDBCSurH66aP//u9r\nePjh7axd25e9e88BLgT6AweAtWSnPrJgMXZsZNYsp492dwYIA4QkNaifPrp48VIeeOANli6tY9Om\n04DpNN6PYitNA8WQIa9RXn6O00e7GQOEAUKSDqt++uijj65k4cKdvPTSSRw8eCFZqBhaaPUi9aGi\nV6/VTJ3aj1mzypw+egIzQBggJKko9dNHlyyp4JFH1vPUUyWF0x7TyJYy7wXsBVZSfy3FWWdt50Mf\nGt9wLcXpp5/uaY8uzgBhgJCkY9J09dEnnniKJ57YxcaN7yMboZgOnFZo+TvqRymGDXuNSy4ZyOzZ\n05g5cybnnXee00e7GAOEAUKSjrumq48uWvQK69b1o65uKtkoxUVAPxov0FxOaelapkw5wBVXjKW8\n3OmjXYEBwgAhSR2u6eqjS5YsZ+nSXVRVnUvjKMWEQsvsAs0QnuLMM99jzpyBzJkzxemjCTJAGCAk\nqdM1XX20srKSJ598li1b3k9joJhG8ws0lzN8+OuUlQWuvvr9XHLJTKeP5swAYYCQpNw1XX00G6Wo\n4IUXamgME9PJLtDsCewBVtC37xrOO283V101hMsvn+r00U5mgDBASFKS6qePVlRUUFFRwVNPraOm\n5mxgRuExExgB1AHrKClZzplnbuHSS0u55pqzmTnT6aMdqbMCRElHvbCOj/nz5+ddQjLsi4z9kLEf\nGnV2XwwdOpRrrrmGO+64gyVLlrBr12YqK/+aO+8cznXX/SNDhkwAzgK+CDxFXd0sXnnlz/jhD2/n\nuusuYOjQRxk9+n/yiU/cwb33/pQ33niD4/HLrN8TnatdASKE8JUQwhshhKoQwvIQwkVHaD87hLAq\nhLA/hPBKCOHz7Su3+/EfRCP7ImM/ZOyHRnn3RWlpKTNmzOD2229nwYIFvPfeu7zwwq/5u7+bwU03\nLWHcuOuAYcB1wH3ASDZv/q/cf/+fcPPNH2Ps2NcYPPgOLrvsT/nf//v/sHr1ampra4uuI+9+6G6K\nntwbQrgR+C5wC/AUMA94JIRwVozxvVbanw48CPwA+AxwGfD3IYSNMcaF7S9dkpSikpISJk6cyMSJ\nE/nSl74ENJ8+Wln531m9eh11decD5cBsdu/+Ko89NpDHHtsBPEnv3v/MhRe+x5VXnu7qo4lqz91B\n5gH3xBh/DBBC+DJwDfAF4K9baf8HwOsxxtsLX78cQigvvI4BQpK6gdGjR3PDDTdwww03AM2nj1ZW\n/g1Ll/4n9u6dAMwBLqWm5pssXdqHpUtfAR6kpOR/MmnSbmbNmt5w10ynj+arqAARQugFTAG+Xb8t\nxhhDCI8CZW3sNh14tMW2R4C7i3lvSdKJY8CAAcyZM4c5c+YALaeP3sPixV9i8+YPAtcCN1JX94es\nXr2T1asX8L3v3QPcyLhxZzSsPFpeXp7n4XRLxY5ADAd6AFtabN8CjG9jn1FttB8UQiiNMVa3sk8f\ngBdffLHI8k48O3fuZPXqDruItkuxLzL2Q8Z+aHQi9UV9GLj99sjGjRtZu3Yta9bcxooVO3n77XHA\nFcCdwDusX/8b1q//d3784x8D0LNnTz72sY/x53/+53keQu6a/Ozs05HvU9Q0zhDCaOBtoCzGuKLJ\n9juBi2OMh4xChBBeBv4xxnhnk21XkV0X0a+1ABFC+Azw02IORJIkNfPZGOPPOurFix2BeA84CIxs\nsX0ksLmNfTa30X5XG6MPkJ3i+CywAdhfZI2SJHVnfYDTyX6WdpiiAkSM8UAIYRXZVS4LAEJ2v9I5\nwN+0sdsy4KoW264obG/rfbYCHZaaJEk6wS3t6Ddoz30g7gK+FEK4KYQwAfgh2fJt9wKEEO4IIdzX\npP0PgbEhhDtDCONDCLcCNxReR5IkdUFFT+OMMf4ihDAc+CbZqYhngCtjjO8WmoyicZF5YowbQgjX\nkM26uA14C/hijLHlzAxJktRFJLkWhiRJSptrYUiSpKIZICRJUtFyCRAhhCEhhJ+GEHaGELaHEP4+\nhND/KPb7ZghhYwhhXwhhYQjhzBbPfymEsKjwunUhhEEddxTF64hFyEIInwwhvFh4zbWFe2wk73j3\nRQjhgyGEfym8Zl0I4baOPYLjowP64fdCCItDCNsKj4VHes1UdEBffDyE8HTh/5g9IYQ1IYTPdexR\nHLuOXKwwhPDpwr+Pfzv+lR9/HfA98fnC8R8sfKwLIezr2KM4dh30s2NwCOH7hZ+p+0MIL4UQPlxU\nYTHGTn8ADwGrgQvJFo9/BfjJEfb5Y2Ab2X1NzwF+BawHejdpcxtwe+FxEBiUx/G1Uf+NZPe0uAmY\nANxTOJ7hbbQ/HdhDtr7IeOArwAHg8iZtZhS2/WGhzTeBauCDeR9vDn1xIdnt6T5FdrOz2/I+zpz6\n4Z+ALwPnka2n/I/AdmB03sebQ19cDHy08PwZhf8fmrVJ7dER/dCi7ZvAE8C/5X2sOX1PfL7w72EE\ncHLhMSLvY82hH3oBTwMPkC038X5gFnBuUbXl0BkTgDrggibbrgRqgVGH2W8jMK/J14OAKuBTrbS9\nhPQCxHLge02+DmQzUm5vo/2dwLMtts0HftPk638GFrRoswz4Qd7H29l90eK5N+gaAaJD+6HwfAmw\nE/hc3sebd18U2qwC/jLv4+3sfih8H1QANwM/omsEiI74P/PzwLa8jy2Bfvgy8CrQ41hqy+MURhmw\nPca4psm2R4EITGtthxDCGWTTQx+r3xZj3AWsoO1FvJIRGhcha1p/JDvuYhcha9q+7CjaJKUD+6JL\n6cR+6E/228a2dhfbwTqrL0IIc8hGZZ48lno7Sgf3wzeALTHGHx2fajtWB/fFgBDChhDC70IIvwoh\nfPA4lX3cdWA/XEfhl80QwuYQwnMhhD8JIRSVCfIIEKOAd5puiDEeJPsPbtRh9om0vihXW/uk5HCL\nkB3umNtchOwIbVLuk47qi66ms/rhTrJTOinfd6XD+iKEMCiEsDuEUEM2XPvVGOPjx6fs465D+iGE\nUE428vB7x6/UDtdR3xMvA18APkK2XEIJsDSEcMrxKLoDdFQ/jAU+SXb8V5Gd/v7vwP8opriibyTV\nlhDCHWTXKbQlAhOP1/tJOrwQwtfJrgm5JMZYk3c9OdkNnA8MILvl/t0hhNdjjIvzLatzhBAGAD8G\nvhRj3J53PXmLMS4nOyUAQAhhGfAi8PtkozTdRQlZqLilMKKxJoTwPuBrwF8d7YsctwAB/C+yc2uH\n8zrZ4lonN90YQugBDOXwC3IFsjtfNk1WI4E1re6Rlo5ahKytNm29Zgo6a0G21HVoP4QQvkZ2MfGc\nGOPzx15uh+qwvij85/h64ctnC8PVfwKkGCCOez+EbLmBMcADIYRQeL4EoDAqMz7G+MbxKP4465T/\nJ2KMtSGENcCZrT2fgI7qh01ATeHfR70XgVEhhJ4xxtqjKe64ncKIMW6NMb5yhEct2XmXk0IIFzTZ\nfQ5ZQFjRxmu/QdYpc+q3hWyK5jQ6YcGQYxVjPEB28VbT+usXIWur/mVN2xe0XISstTaXc5iFyvLW\ngX3RpXRkP4QQbicbiryyxbVGSerk74kSIMnTXh3UDy8B5wKTyEZizidbCPHxwudvHqfyj6vO+p4o\nnPM/l+wHanI6sB8qOTQ0jQc2HW14qC8wj6tKfwOsBC4CZpKdl/qnFm1eAj7a5Ovbga1kF3+cSzaN\n81WaT+McSfaP4vfIZnqUF74eksdxtjieTwH7aD4VZyuFKUTAHcB9TdqfTjb8emfhL/ZWoAa4rEmb\nMrJpm/XTOP+CbLpP6tM4O6IvehX+rieRnfO/s/D1uLyPt5P74Y8L3wMfL/x7qH/0z/t4c+iLrwOX\nkU3hnEB2jrcauDnv4+3MfmjlPbrKLIyO+J74M7Jfss4ALiCbnbAXmJD38XZyP7wP2EG2ivYHgGvI\nfkn/elG15dQhJwE/IZteth34v0C/Fm0OAje12PYXZNM595FdVXpmi+e/QRYcDrZ43NSRx1PEcd8K\nbCCbfroMuLDJcz8CHm/R/mKy9FlFFpb+UyuveT1Z2KoCniX7rTP3Y+3sviAbpm3t7/7xjj6WxPrh\njVb64CDw53kfaw598Vdkv5zsJRsKrgBuyPs4O7sfWnn9LhEgOuh74q7Cv5Eqsp8lDwDn5X2ceXxP\n0DiCv6/Q5o8prI91tA8X05IkSUVzLQxJklQ0A4QkSSqaAUKSJBXNACFJkopmgJAkSUUzQEiSpKIZ\nICRJUtEMEJIkqWgGCEmSVDQDhCRJKpoBQpIkFe3/AWf002waE3LgAAAAAElFTkSuQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# To plot the initial condition of the chain, run this cell.\n", "plt.plot(f(u), z, 'k', lw=2)\n", "xfit = np.zeros(N)\n", "for n in range(nmax):\n", " xfit += A[n-1] * j0(w[n-1] / b * u)\n", "plt.plot(xfit,z)\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Set up the Figure and Axes\n", "fig, ax = plt.subplots(figsize=(4,6))\n", "\n", "ax.axis('off')\n", "ax.set_xlim((-2*d, 2*d))\n", "ax.set_ylim((0, L))\n", "\n", "line, = ax.plot([], [], 'k', dashes=[5,2], lw=3)" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# Initialization function for the animation\n", "def init():\n", " line.set_data([], [])\n", " return (line,)" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Delay between frames (ms)\n", "interval = 1000/24\n", "# Total length of animation (s)\n", "anim_duration = 5\n", "# Total number of frames in the animation\n", "nframes = int(anim_duration * 1000 / interval)\n", "\n", "# The animation function, to be called sequentially\n", "def animate(i):\n", " t = i * interval / 1000\n", " x = np.zeros(N)\n", " for n in range(0, nmax):\n", " x += A[n-1] * j0(w[n-1] / b * u) * np.cos(w[n-1] * t)\n", " line.set_data(x, z)\n", " return (line,)" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Set up the animation\n", "anim = animation.FuncAnimation(fig, animate, init_func=init,\n", " frames=nframes, interval=interval, blit=True)" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# To save the animation as a gif, run this cell.\n", "anim.save('chain.gif', writer='imagemagick', fps=1000/interval)" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# To show the animation in a Jupyter Notebook cell, run this cell.\n", "HTML(anim.to_html5_video())" ] } ], "metadata": { "anaconda-cloud": {}, "kernelspec": { "display_name": "Python [conda env:py35]", "language": "python", "name": "conda-env-py35-py" }, "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.5.2" } }, "nbformat": 4, "nbformat_minor": 1 }