{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Duffing Oscillator Solution" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "%matplotlib inline\n", "%load_ext autoreload\n", "%autoreload 2\n", "import scipy as sp\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import mousai as ms\n", "from scipy import pi, sin" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Equation errors (should be zero): [[ 1.11022302e-15 -1.68753900e-14 1.33226763e-15]]\n", "Constant term of FFT of signal should be zero: (-0.10771053458129715+0j)\n" ] } ], "source": [ "# Test that all is working. \n", "# f_tol adjusts accuracy. This is smaller than reasonable, but illustrative of usage. \n", "t, x, e, amps, phases = ms.hb_time(ms.duff_osc, sp.array([[0,1,-1]]), .7, f_tol = 1e-12)\n", "print('Equation errors (should be zero): ',e)\n", "print('Constant term of FFT of signal should be zero: ', ms.fftp.fft(x)[0,0])" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Equation errors (should be zero): [[ 2.04616314e-08 -1.11916120e-08 -4.32641484e-09 -8.27785285e-10\n", " -1.06717036e-08 -2.49418783e-07 -6.82304595e-07 -1.59994936e-07\n", " -3.88299430e-11 9.85229565e-09 1.56362223e-09 -1.71778813e-10\n", " 7.32659703e-08 4.91827027e-07 5.65568172e-07]]\n", "Constant term of FFT of signal should be zero: (7.04098108967e-06+0j)\n" ] } ], "source": [ "# Using more harmonics. \n", "t, x, e, amps, phases = ms.hb_time(ms.duff_osc, x0 = sp.array([[0,1,-1]]), omega = .7, num_harmonics= 7)\n", "print('Equation errors (should be zero): ', e)\n", "print('Constant term of FFT of signal should be zero: ', ms.fftp.fft(x)[0,0])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Sometimes we can improve just by restarting from the prior end point. Sometimes, we just think it's improved. " ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Errors: [[ 7.29971639e-15 2.38697950e-15 1.02140518e-14 1.55431223e-14\n", " 2.88657986e-14 -4.44089210e-16 -1.30007116e-13 -1.18793864e-14\n", " -5.61399885e-15 -1.65978342e-14 -2.33146835e-15 -1.99840144e-15\n", " -2.56461519e-14 7.13873405e-14 1.18960397e-13]]\n", "Constant term of FFT of signal should be zero: (7.08481329931e-06+0j)\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "/Users/jslater/anaconda3/lib/python3.6/site-packages/scipy/optimize/nonlin.py:474: RuntimeWarning: invalid value encountered in double_scalars\n", " and dx_norm/self.x_rtol <= x_norm))\n" ] } ], "source": [ "t, x, e, amps, phases = ms.hb_time(ms.duff_osc, x0 = x, omega = .7, num_harmonics= 7)\n", "print('Errors: ', e)\n", "print('Constant term of FFT of signal should be zero: ', ms.fftp.fft(x)[0,0])" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The average for this problem is known to be zero, we got 4.72320886587e-07\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD7CAYAAABpJS8eAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8VNX9//HXyTZJyEaABAgJEJBN9gyLqDVoUFFcy+Ze\nrMb6035pbWtba1u//doqtsVv+7VaUaBuLIoblSoSNAgIQsKO7AkJAUIgISHrTGbm/P7IRAKEJExm\ncmf5PB8PHpmZOzP3k0Pyzp1zzj1Xaa0RQgjh/4KMLkAIIUTHkMAXQogAIYEvhBABQgJfCCEChAS+\nEEIECAl8IYQIEBL4QggRICTwhRAiQLgl8JVSo1vYNsf5NdMd+xJCCOEa1d4zbZVSGcCrWut+F9l+\nGigDHtFaZ7X0Xl27dtV9+vRxqY7q6mo6derk0mv9kbTHuaQ9ziXtcSFfbZPc3NxTWutubXluSHt3\nprXOUkrltfCUh7XWy9ryXn369CEnJ8elOrKzs0lPT3fptf5I2uNc0h7nkva4kK+2iVKqoK3P7Yg+\n/FSlVIZS6skO2JcQQoiLaHeXDoBSapXWelIrz5kDrDq/W8fZt58JkJiYmLZkyRKXaqiqqiIqKsql\n1/ojaY9zSXucS9rjQr7aJhMnTszVWpvb8tx2d+m0xBnmZc4unVIg9fznaK3nAfMAzGazdvUjla9+\nHPMUaY9zSXucS9rjQoHQJh7p0lFKxTlv5gCNR/T9nPeFEEIYoN2Br5SaCpidXxutBtBabwGmO7cd\nct4XQghhAHfM0lkGLDvvsbQmt+e1dx+iA1QWw7JZMPVfEJ1odDVCCA+QM21FgzUvQOFGWDPH6EqE\nEB7i0UFb4f30/ySg7JazD+TMh5z5WAljWrePMAUHEREWTGRYMN1jw0nuHElyfCRDesbQMzYcpZRx\nxQshLokEfoCx2Ox8fbCUL/aWkFNwmrKaufw6+B1uCM4hQlmpI4wNYRNYFJtJTHgIVpuD8horRadt\nrNl/khqr/bv36hZtYmRyHFf178p1gxPo1TnSwO9MCNEaCfwAoLVm3cFTLMst4os9JVRabESGBTM6\npTOTrh3LmBPrCD+4AR1iItxez8Th/Zg45eZm36es2kpBWQ27jlawrbCc3MLTrPr2BL9fvptB3aO5\ndWRPpo7uRUJMuIwLCOFlJPD9WJXFxns5R3hrYwF5J6vpHBnKTcN6cOPQ7kzo3wVTSHDDE5dUgvlB\nMM+CnIVQdaLZ91NK0SXKRJcoE6NTOnP/FQ2P552sYvWeElbuLuaFz/bx18/3c92gBP47eAHdCzei\n1syBKXM76LsWQlyMBL4fqqu382l+PU+s/ZKyaiujUuJ4ccYIbhrW42zINzXznbO3XQjm1G5RpHaL\n4uHvpZJ3soqUl/sSkmc9+wTnuAAhJni6xIXvSAjhDhL4fkRrzbLcIv7y+T5OnLFy9WVd+emkAYxO\n6dxhNaR2i4IndmL77Dew5xNCHHXU6jC+CZ9AlzteYFiHVSKEOJ8Evp84WFLJUx/uYlN+GaNS4pg1\nSPGjO8cZU0x0d0LCY0BbG8YFbFZOWk384F8HuXN0Lb+6cVBDH78QokPJPHxfVVkMCydjP1PM/60+\nwOS/rWVfcSXP3TmM9380gUHxzXTddKTqEkibhXpoNcr8ILdfFsKj6f34ZPtxMuau4YMtRbhj4T4h\nRNvJEb6vWvMCumADWf98gr+W3cOU4T145tbL6RplMrqyBueNC4QCvwSmm5P5xXvbeeLd7Xy6q5g/\n3jGUhGg52heiI0jg+5pnE8DWcKKUAm6oWcHh8BXoPBMqyvsHRPt27cTSR65gwbp8/vz5Pm7621r+\nftcoJvTranRpQvg96dLxMfq/trM/4UZqdRgAjpAIGDYNNXunwZW1XXCQ4uHvpfLJj68iNiKUe1//\nhpezD+JwSBePEJ4kge9DLDY7P/vsBJuO1WNS9egQE0F2C5hifPLEpgGJ0Xz8+FVMHtaDFz7bx4/e\nzqXGajO6LCH8lgS+j6isq+e++Zv4YMtRxic6UOYHUQ+thrRZFz1RyhdEmUJ46a5R/HbKELL2nOCu\neRs5WWlp/YVCiEsmffg+oKzaygMLNrHn+Bn+NnMk/Ud+dHajH5zBqpTih1f1JSU+kh8v3sKdr6zn\nX7PG0q+b711uTghvJkf4Xu7EmTpmvLqB/ScqmXd/GreNTDK6JI+ZNCSRJZlXUGu1M/WVr9l1tMLo\nkoTwKxL4XqzEGfbHymt548GxXDvI9/rpL9XI5Djef3QCkWEh3P3aRrYfKTe6JCH8hgS+lzpdbeW+\n+ZsoqbTw5g/HMT61i9EldZjeXTqxJHM8Mc4ZPLkFp40uSQi/IIHvhSrr6nlg4SbyS6t5/X4zab07\nbi0cb5EcH8m7j1xBfFQYDyzYJN07QriBBL6XsdjsPPRGDt8eO8PLd49mQv/APSGpZ1wESzLHExsR\nygMLNlFw+BAsnAyVvjsrSQgjSeB7Ea01v35/J9/kl/GXaSPIGOL/ffat6REbwVs/HAtA7ltPoQs2\nyHV3hXCRTMv0Ii99cZAPth7liUkDuH2U/87GuVSpr6aS28x1d2V9fSEujRzhe4nl24/x11X7uWNU\nEj++tr/R5XiX2Ttg6DTswQ2LrFmUCcfQaeBDy0kI4Q0k8L3A7mMV/OK97Yzp05nnvz8MpZTRJXmX\n6O5giibYYcUeFEaow0rOcZtPLichhJEk8A1WUVvPo29vIS4ylFfuTWv+EoTiu/X1gzO/YEvCHZSV\nFPHWhsNGVyWET5E+fAM5HJqfvbuNY+W1LH1kvPesZe+NmqyvP+rRBWS+mUP2v79lYPcYxvaNN7Aw\nIXyHW47wlVKjW9g2VSmVoZR60h378if//OoQWXtKePrmwaT1ltBqq+Agxf/OHElKfCSPL9oii60J\n0UbtDnylVAbw3kW2jQbQWmcB5S39YQg0uQWn+cvKfdwyoicPTOhjdDk+Jzo8lJfvGU1FbT3/tXgr\ndllLX4hWtTvwnWGed5HNM4DGxVDygIz27s8fVFls/HTpNnrGRfCnO4bKIK2LBveI4dnbh7Ihr5QX\nV+03uhwhvJ6nB23jgLIm9wNnQZgW/Pfy3RSdruHFGSOJDg81uhyfNs2czAxzMi99eZB1B04ZXY4Q\nXs3wQVulVCaQCZCYmEh2drZL71NVVeXyazvS5mIb722zcEtqKNWHd5B92DP78ZX2cIdrO2u+6qT4\n8dub+J8rI4gKu/ATUyC1R1tIe1woENrE04FfDjSORsYBpec/QWs9D5gHYDabdXp6uks7ys7OxtXX\ndpSTlRZ+8uIahveKZe4PJxAa7LkPWL7QHu7Uc2AFd7y8nk9PxfLSXaMu6CYLtPZojbTHhQKhTTyS\nOEqpOOfNpUCq83YqkOWJ/fmEymLKX84g0lLK3OkjPBr2gWhYr1h+OmkAK3Yc56NtR40uRwiv5I5Z\nOlMBs/Nro9UAWustzudkAOWN9wPRkY+eoV/NTl5NyaJ/QrTR5filH13TjzF9OvO7j3ZztLzW6HKE\n8DrumKWzTGvdWWu9rMljaU1uz9NaZzm7bgLPswnwTCzJhxYTpDTDji+DZ2IbHhduFRykmDt9JDaH\n5ukPd6K1TNUUoinpV/C02TvYGjeJWh3WcD8kAobJwl+ekhwfyc9vGMiX+06yfPsxo8sRwqtI4HtY\nblkYu085MKl6CAkHuwVMMbLwlwf9YEIfRibH8czy3ZRWyVm4QjSSwPcgu0Pzu493kxRWhX30LHgo\nC9JmQZVcscmTgoMUc74/nCqLjT988q3R5QjhNQyfh+/PFm8qZPexM1Tf/S9Ch/dseHDKXGOLChAD\nu0fzaHp//r76AHeO7mV0OUJ4BTnC95Cyait/XrmPK1K7cPOwHkaXE5Aem9iPvl078czy3dTLWjtC\nSOB7yp9X7qPKYuO/b7tc1soxiCkkmN/fMoT8U9WsPFxvdDlCGE4C3wN2H6tgyeZCfjChDwMSZc69\nkdIHJjBpSCLLD9VzTObmiwAnge9mWmue+89e4iJC+a/rLjO6HAH8bsoQtIY//meP0aUIYSgJfDdb\ns/8k6w6e4sfXXkZshKyE6Q2S4yOZkhrKih3H+fqQrKgpApcEvhvZHZrnP91LSnwk947vbXQ5oonJ\nfUNJiovgT//Zg0MGcEWAksB3o/e3FLG3uJInbxxIWIg0rTcJC1b84oaB7Dp6ho+3y+JqIjBJKrlJ\nrdXO3M/3MyI5TqZheqlbR/RkaFIMf1m5n7p6u9HlCNHhJPDd5K2Nhyk+U8dTkwfJNEwvFRSkeOqm\nwRwtr2Xh+sNGlyNEh5PAd4Mqi41/rsnj6su6Mi5VruLozSb068p1gxJ4+cuDss6OCDgS+G7wxteH\nKau28sSkAUaXItrgV5MHUW218Ur2IaNLEaJDSeC305m6euZ9lce1gxIYldLZ6HJEG1yWGM0do3rx\n1sYCSo4WwMLJUCkL2gn/J4HfTvPX5lNRWy9H9z7mJxmXYXdo8j/4PRRuhDVzjC5JCI+T1TLbobzG\nyoJ1+Vw/JJGhSbFGlyMuQfLLfTgYZoFS5wM58xv+hZjg6RJDaxPCU+QIvx0Wrj9MpcXGT+Xo3vfM\n3kHNwDvlSmQioEjgu6jKYuNfXx8mY3Aig3vEGF2OuFTR3YmMisOk6qnToWhbnVyJTPg9CXwXvbOx\ngIraeh6b2M/oUoSrqkuwjPgBMxzPsr7zbXIlMuH3pA/fBXX1dl5bm8+V/bvIzBxfNvMdIgBz6Lc8\n8HUK2felk2x0TUJ4kBzhu+C9nCOcqrLw2MT+Rpci3CDze6kEK8Ura2RevvBvEviXqN7u4J9r8hiV\nEscVclatX0iMCWf6mF4syynieIVcJEX4Lwn8S/Tv7cc4Wl7LY+n9Zc0cP/LI9/rh0Jp5X+UZXYoQ\nHiOBfwm01ry+Np/LEqK4bnCC0eUIN0qOj+SOUUks3lTIyUpZY0f4Jwn8S7Ahr5Rvj5/hh1f1laN7\nP/Roej+sNgcL1ucbXYoQHtHuwFdKTVVKZSilnrzI9jnOr5nt3ZfR5q/Np0unMG4flWR0KcIDUrtF\nccPl3XlnYwHVFpvR5Qjhdu0KfKXUaACtdRZQ3nj/PJlKqUOAT3eOHiypYvXeEu67ojfhocFGlyM8\n5KGrUzlTZ+O9nCNGlyKE27X3CH8GUO68nQdkNPOch7XW/Zx/FHzWgvX5hIUEybVq/Vxa786k9e7M\n/PX52OwOo8sRwq3aG/hxQFmT+83NU0xtqcvHF5RVW3k/t4g7RyXRNcpkdDnCwx6+OpUjZbWs3C1n\n3gr/4vEzbbXWLwAopSYppTLOP9J39u1nAiQmJpKdne3Sfqqqqlx+bWuWH7JisTkYZjrlsX24myfb\nwxddSnuEaU1CpOKvK7YRWbrXLwfo5efjQoHQJu0N/HIg3nk7jrOLzQLfhXmZ1nqZc1vq+W+gtZ4H\nzAMwm806PT3dpUKys7Nx9bUtsdjs/Hzdl6QP7MY9U8a6/f09xVPt4asutT1+HH6Y3368m6i+IxjT\nJ771F/gY+fm4UCC0SXu7dJZyNsRTgSwApVSc87GcxseAfs77PuXjbcc4VWXhoasu+Fsl/NjUtGTi\nIkN5TU7EEn6kXYGvtd4CoJTKAMob7wOrm2yfrpSaChxqst0naK1ZsC6fQd2jubK/LKMQSCLCgrlv\nfG9W7TlB/qlqo8sRwi3aPQ9faz1Pa53l7JppfCztvO3LGvvyfck3+WXsLa7kwSvlRKtAdN8VvQkN\nCmL+OjnKF/5BzrRtwVsbCoiNCOWWET2NLkUYICE6nDtGJfFeThGnq61GlyNEu0ngX0RxRR0rdxcz\nY0wyEWFyolWgmnVVHyw2B+/KiVjCD0jgN6eymPr5N9JZn+becXKiVSAb1D2GcX3jeWtjAXaHNroc\nIdpFAr8Z9uw5JJ3ZxvNdPiWlS6TR5QiDPTChD0Wna/lib4nRpQjRLnKJw6aeTQCbhcYOnOuq/g3P\nxEKICZ6WX/ZAdf2QRHrEhvPmhsNMGiIXORe+S47wm5q9A4ZOw0LD8gk6JAKGTYPZOw0uTBgpJDiI\ne8alsPbAKQ6WVBldjhAuk8BvKro7ZXYTodqKLSgMZbeAKQai5agu0M0cm0JYcBBvbThsdClCuEwC\n/zzHjhayWGdQc/9KSJsFVbKAloCuUSamDO/BstwiKuvqjS5HCJdI4DdRUVvPtNOPsW3Yb4npMxqm\nzIWZ7xhdlvAS90/oQ7XVzodbjxpdihAukcBv4oMtRdTW23lgQh+jSxFeaGRyHCN6xfLG14fRWqZo\nCt8jge+kteadbwoZkRzH0KRYo8sRXuqBCX04dLKa9QdLW3+yEF5GAt8pp+A0B0uquGdsitGlCC92\n8/AedOkUxhsbDhtdihCXTALfafE3hUSbQpgyoofRpQgvZgoJZpo5mS/2llBcUWd0OUJcEgl8oLzG\nyic7j3P7qCQiw+RcNNGyu8YmY3dolm6W9XWEb5HABz7YchSrzcFd0p0j2qB3l05cfVlXlm4ulPV1\nhE8J+MDXWrNoUyEjk+MY0jPG6HKEj7h7bArHKurI3idLbgjfEfCB3zhYe7cc3YtLkDEkkW7RJhZ9\nU2h0KUK0WcAH/iIZrBUuCA0OYrq5F1/uK+FYea3R5QjRJgEd+OU1VlbIYK1w0cwxKWhgiQzeCh8R\n0IEvg7WiPZLjI/neZd1YurkQm91hdDlCtCpgA18Ga4U73DMuhRNnLHJxFOETAjbwvxusHSdH98J1\n1w5KoHtMOIs2yeCt8H4BG/jfDdYOl8Fa4bqQ4CCmj0lmzf6THCmrMbocIVoUkIHfOFh7x2gZrBXt\nN3NMMgpYslmO8oV3C8jAl8Fa4U494yKYODCBd3OKqJfBW+HFAi7wtdYsdg7WDu4hg7XCPe4el8LJ\nSgtZ38oV0oT3anfgK6WmKqUylFJPurK9o+UWnOaADNYKN0sfmEDP2HBWbNgKCydDpQS/8D7tCnyl\n1GgArXUWUN54v63bjbBokwzWCvcLDlLMGJPCuML56IINsGaO0SUJcYH2jljOAFY5b+cBGcCWS9je\noSpq6lmx4zjTzckyWCvc69kEZtssZ3+jcuY3/AsxwdMyR194h/Z26cQBZU3ud7nE7R3qw61FWGwO\nZo5NNrIM4Y9m74Ch07AoEwA6JAKGTYPZOw0uTIizDD/MVUplApkAiYmJZGdnu/Q+VVVVLb5Wa83r\n62vpGxvEyf1byd7v0m58RmvtEWg6oj0uK62kp7ZSp0MJs9Vx7FQlB3L3AHs8ul9XyM/HhQKhTdob\n+OVAvPN2HHD+lZ1b247Weh4wD8BsNuv09HSXCsnOzqal1+YWnKZo5dc8f+dQ0gNgOmZr7RFoOqQ9\nil/DkfQgmduHcH9YNhmxmiQv/T+Qn48LBUKbtDfwlwJm5+1UIAtAKRWntS6/2HYjLN5USKewYG4Z\n0dOoEoS/m/kOQUBa+AEeyurB2h9ORDoPhTdpVx++1noLgFIqAyhvvA+sbmV7h6qoreeTHce4bVQS\nnUyG92IJPzd9TC+ClJx5K7xPu9PP2SVz/mNpLW3vaB9vO0pdvUOuaiU6RI/YCK4d1HDm7U8yBhAa\nHHDnNwov5fc/iVprFn1TyLCkWIYmxRpdjggQjWfert4jJ2AJ7+H3gb/tSDl7iytl3RzRoa4ZkECP\n2HAWbZKrYQnv4feBv3hTIZFhwdw6UgZrRcdpOPM2mbUHZNlk4T38OvDP1NXz7+3HuW1kT6JksFZ0\nsBmybLLwMn4d+B9vO0ZtvV26c4Qhmg7eyrLJwhv4beA3DtZe3jOGYTJYKwxy19jGwVtZT0cYz28D\nf0dRBXuOn+GusSkopYwuRwSoawZ0o0dsOIvlmrfCC/ht4C/eVEhEaDC3yWCtMFBIcBDTzcl8JYO3\nwgv4ZeBX1tWzfPsxbhnRg+jwUKPLEQGucfB26WaZoimM5ZeB/9G2Y9RYZbBWeIez17w9IoO3wlB+\nF/haa97eUMDQpBhGJscZXY4QQMPgbYkM3gqD+V3gbz58mn0nKrlvfG8ZrBVeI32gDN4K4/ld4L+1\nsYCY8BBuHZFkdClCfEcGb4U38KvAL6ms47Ndx5lmTiYiLNjocoQ4x3Tn4O27OTJ4K4zhV4G/dNMR\n6u2ae8bJYK3wPklxEaQPTGDpZhm8Fcbwm8C3OzSLNhVy9WVdSe0WZXQ5QjTr3vENg7ef75Zlk0XH\n85vA33bSzvGKOu4d39voUoS4qGsGJJAcH8EbGw4bXYoIQP4R+JXFXL/vNwyNqeW6QQlGVyPERQUH\nKe4d15tN+WXsLT5jdDkiwPhF4Fes/COX2/fxXJdPCZHLyQkvN92cjCkkiLc2FBhdiggwvp2OzybA\nM7HE7nqTIKUZdnwZPBPb8LgQXqpzpzBuHdGTD7ce5UxdvdHliADi24E/eweOoVOpI6zhfkgEDJsG\ns3caW5cQrbj/ij7UWO18kFtkdCkigPh24Ed3J8gUgwkbdhUKdguYYiA60ejKhGjRsF6xjEyO482N\nBWitjS5HBAjfDnyA6hKUeRZb0v4MabOgSqa7Cd9w/xW9yTtZzfqDpUaXIgzkcOgO+6Pv+4E/8x2Y\nMpfqqL4wZW7DfSF8wE3DehDfKYw3Nxw2uhRhoGVbirjtH+sprbJ4fF++H/hC+Kjw0GBmjEkma88J\njpbXGl2OMIDWmgXr8rHUO4jvFObx/UngC2GgxmVA3t4oUzQD0YZDpewtruTBq/p0yOq+EvhCGKhX\n50huuLw7i74ppMZqM7oc0cEWrM+nS6cwbhvZMav7tjvwlVJTlVIZSqknL7J9jvNrZnv3JYQ/eujq\nvlTU1vO+TNEMKIdOVrF6bwn3jO9NeGjHrO7brsBXSo0G0FpnAeWN98+TqZQ6BOS1Z19C+KvRKZ0Z\nmRzH/HX5OBwyRTNQvL42j7DgIO6/ouPW/2rvEf4MoNx5Ow/IaOY5D2ut+zn/KAghzqOU4qGr+3K4\ntIasPTKtOBCcrLTw/pajfD+tF12jTB2235B2vj4OKGtyv0szz0lVSmUAo7XWL5y/0dnVkwmQmJhI\ndna2S4VUVVW5/Fp/JO1xLm9vjwiHpku44q+fbCXsZITH9+ft7WGEjmyT9/dbqbc5GB52skP/H9ob\n+K1qDHml1CSlVMb5R/pa63nAPACz2azT09Nd2k92djauvtYfSXucyxfa49HQPJ5dsYf4/iMZ3ivO\no/vyhfboaB3VJtUWG7PXfMH1lycy82azx/fXVKtdOkqpzGb+NXbdlAPxzttxQGkzr53qvFsKpLqr\ncCH8zYwxyUSZQpi/Lt/oUoQHLd18hIraejK/16/D993qEb7zCPxilgKNf6JSgSwApVSc1rocyOHs\nYG0/4FXXSxXCv0WHhzJzTDILvz7ML28cRM84z3ftiI5Vb3cwf10+5t6dSevducP3365BW631FgDn\nEX95431gdZPt051H+YeabBdCNOMHV/ZBa80bXx82uhThAcu3HeNoeS2PXNPxR/fghj785j4BaK3T\nWtouhGher86RTB7Wg0WbCnn82v5Eh4caXZJwE7tD84/sgwzqHm3YlfnkTFshvEzm1alU1tlY9E2h\n0aUIN/p013HyTlbz+LX9CQry/DIKzZHAF8LLjEiO46r+XXltbT519XajyxFu4HBoXvriIKndOjF5\naA/D6pDAF8ILPTaxP6eqLLyXc8ToUoQbrN5bwt7iSh5L70+wQUf3IIEvhFcanxrP6JQ43svOwbFg\nMlTKGbi+SmvNS18cIDk+gltH9jS0Fgl8IbyQUorHr+3PtOrFqMINsGaO0SUJF2XvP8n2ogoevaY/\nocHGRq7Hz7QVQrjg2QSutVnO/obmzG/4F2KCp0sMLU20ndaav36+j+T4CKam9TK6HDnCF8Irzd4B\nQ6dhCw4HaPg6bBrM3mlwYeJSrNxdzK6jZ5h93QDCQoyPW+MrEEJcKLo7mKIJdlixEEqQ3YIjLBqi\nE42uTLSR3aGZu2o/qd06cbvBffeNJPCF8FbVJai0WWy+7j3etl3H8WMyL9+XfLLjGPtPVPHEpAGE\nGNx330j68IXwVjPfAWCCQ/PHLSEsOGMjy+7wmvAQF2ezO3hx1X4G94jhJgPn3Z9PfnKE8HJBQYqf\nZlzG4dIaPth61OhyRBss2XyEw6U1/GzSAMPOqm2OBL4QPmDSkESGJcXy99UHsNocRpcjWlBZV8//\nZu1nbN94rhtszJo5FyOBL4QPUErxxKQBFJ2uZfEm6cv3Zq+uyeNUlZXf3DQYpbzn6B4k8IXwGekD\nuzGubzx/X32Ayrp6o8sRzTheUctra/O4bWRPRiR79qplrpDAF8JHKKV46qbBlFZbmfdVXusvEB3u\nzyv3oYFf3DDQ6FKaJYEvhA8ZkRzHLSN68traPE6cqTO6HNHEzqIKPtx6lAev7EuvzpFGl9MsCXwh\nfMwvrh+I3aF5cdV+o0sRTg6H5rcf76JLpzD+30RjrmbVFhL4QviYlC6R3De+D+/mHGFfcaXR5Qhg\nWW4R246U86vJg4nx4quUSeAL4YN+7Lz84TPLd6O1NrqcgFZRU8/zn+0lrXdn7hyVZHQ5LZLAF8IH\nde4Uxs9vGMiGvFJW7DxudDkB7a+r9lFeY+UPt13uVSdZNUcCXwgfdffYFIb0iOGPK/ZQY7UZXU5A\n2nW0grc3FnDv+N5c3jPW6HJaJYEvhI8KDlL84bbLOV5Rxz++PGh0OYGlshjHgsk8/94a4juZ+Nkk\n75yGeT4JfCF8mLlPPHeOSuK1r/LJP1VtdDmBY80LqMIN3HDqDZ69/XJiI713oLYpCXwhfNyvJg/C\nFBrEUx/slAFcT3s2AZ6JhZz5KDT3hWRx47JBDY/7AAl8IXxcQkw4v548mA15pbybc8Tocvzb7B3o\noVOxYAJAh0T41JXI3BL4SqnRLWybqpTKUEo96Y59CSEuNHNMMuP6xvPsij2UyBm4nhPdnX2nIVRb\nsQeZUHYLmGJ85kpk7Q58pVQG8N5Fto0G0FpnAeUt/WEQQrguKEjx/PeHY7U5+N3Hu40ux28dOFFJ\n4ZECvoyeQtDDWZA2C6pOGF1Wm7U78J1hfrGVnGYA5c7beUBGe/cnhGhe366d+EnGAD7bXcx/ZG6+\n21lsdmYzfS7WAAAK+0lEQVQv2cavQ3/J8Efmo3oMhylzv7symS/wdB9+HFDW5H4XD+9PiID20NV9\nGd4rlqc+3CmLq7nZ3FX7+fb4GZ7//nC6RZuMLsclhl/TVimVCWQCJCYmkp2d7dL7VFVVufxafyTt\nca5Aao+7+zr43fF6Hnz1S35mDieomYtwBFJ7tFVLbbKn1M68zXWk9wohtGQP2SV7OrY4N2k18J2B\nfL48Z1dOa8qBeOftOKD0/CdorecB8wDMZrNOT09vw9teKDs7G1df64+kPc4VaO1h71rAbz7cRUFY\nH2Zd2feC7YHWHm1xsTYpOVPHz/++jr7dOvFy5lVEhhl+nOyyVit3BvIlUUrFaa3LgaWA2flwKtCW\nPxJCiHa6e2wKX+wp4blP9zKhX1cGdo82uiSfZLM7eHzxVqotNhY9PM6nwx7cM0tnKmB2fm20GkBr\nvcX5nAygvPG+EMKzlFLMmTqcmPBQHn0nlyqLrLXjir98vp9N+WU8d+cwBiT6/h9Nd8zSWaa17qy1\nXtbksbQmt+dprbNc+aQghHBd1ygT/3fXKA6fquaX7++Qs3Av0crdxfxzzSHuGZfC7V6+7HFbyZm2\nQvixK/p14ckbB7Fix3EWrj9sdDk+Y/exCn66dBsjesXy2ylDjC7HbSTwhfBzj3wvlUlDEvnTf/aw\n+XBZ6y8IcCWVdTz8Rg4x4aG8dr+Z8NBgo0tyGwl8IfycUoq/TBtBcnwkj7yVS1FhPiO3PgWVvnOG\naEepq7fz8Ju5nK6p5/UHzCTEhBtdkltJ4AsRAGIjQpn/gBm7Q5P71q+IrfgW1swxuiyv4tCany7d\nxvYj5bw4YyRDk7z/giaXyrfnGAkh2iz11VS2awvUOx/Imd/wL8QET5cYWpvRtNb8a7eVr4qKefrm\nwdw4tLvRJXmEHOELEShm74Ch07AFN3RTWJUJ7UNL+3rSCyv38VWRjccn9uehq1ONLsdjJPCFCBTR\n3cEUTYjDSj2hhDisbDpmQ0f5xsU7POWV7EO8kn2I9OQQfnb9AKPL8SgJfCECSXUJpM1ia9oLbE28\nk9MlRcz5bF/AztF/6YsDzPlsL7eO6Mn9Q8JQzaw75E8k8IUIJDPfgSlzqYlOZfSj81mX9iL/XHOI\n/806EFChr7XmxVX7+cvn+7ljVBJzp49odpE5fyODtkIEKKUUf7h1KJZ6B39bfYBqi43f3DzY749y\nHQ7NnM/28upXeUxL68Xz3x9OcJB/f8+NJPCFCGBBQYo53x9OJ1MIr6/Lp7LOxp/uHOZ/AVhZDMtm\nYbljPr/49ATLtx/j3vEp/OHWoQT52/faAgl8IQJcUJDi97cMISYilL+vPkB5rZUXZ4z0+ZUhz7Hm\nBXTBRrLn/ZzlZffw5I0DefSafn7/aeZ8fvQ/KoRwlVKKJyYNoHNkKP/zybdMfWUDrz1gJikuwujS\n2ufZBLBZAFDADTUrOBy+AtaZID3wzj2QQVshxHdmXdmX+T8Yw5GyGm57aT25BT6+9s7sHRxJupla\nHQaAPTgcAvjcAwl8IcQ5Jg5M4MPHJtDJFMyMVzfyjy8PYnf43gyeGquNp1efZE1BHSZVjw42Eeyw\ngikGohONLs8QEvhCiAv0T4hm+eNXccPQ7vx55T7uff0biit856LouQVl3PS3tby9sRBzVxs67UHU\nw6shbRZUBe6icdKHL4RoVmxEKC/dNYprBnTjmeW7mfTiGp68cRB3j03x2lk8lXX1/C3rAPPX55MU\nF8Hih8czqN/NZ58wZa5xxXkBCXwhxEUppZhuTmZMn3ie/mgnv/1oF8tyi/jj7UO9ajVJh0Pz/paG\ns4ZPVVm4e1wKT900mCiTRFxT0qUjhGhV366dePuH4/jbzJEcPV3DlP9bx2OLtnCwpKphjvvCyYas\nr6+1ZtW3J7jtH+v5xbIdJMdH8PFjV/KnO4ZJ2DdDWkQI0SZKKW4bmUT6wARe+yqPBevz+XTncd5M\nXMqV5RtQa+Z0WJeJ1ebgs93FvPzlQfYWV5ISH8nc6SO4fWRSQJ1Idakk8IUQlyQ2IpSf3zCQn226\nGhVkgXLnBuf6+jrYhPptO+e4O8+MZeq/zplRc7Ckkndzing/t4jSaiup3Toxd/oIbh3Rk5Bg6bBo\njQS+EMIl6ic7YOXT6L2foGy11GHiU5uZOdZ7GbBgE9cPSWRc33j6dYu69KPuNS9A4Ubqv3yOzZc/\nTfa+k2R9e4K8U9WEBCkyBicyc2wyV1/WzWsHkL2RBL4QwjXO9fWV3QIh4ZjsVq4a2pe9MaP4dFcx\nT3+0C4C4yFBGJseR2jWKPl0jSe4cSUxECJ1MIUSGhmC126mrd1BtsWF+ezDBDst3uwjdspAJWxYy\nWoeyJ2UFD0zow+Rh3UmI9q9rzXYUCXwhhOuc6+tjnoXKWUi3qhP8+qbB/GryIPJOVZNbcJrcw6fZ\nXlTON3ll1NbbW3y7bszlNyHvcENIDhFYqQ8yUZp8PdG3zuGtLkkd9E35Lwl8IYTrZr5z9naTAVul\nFP26RdGvWxTTzclAw4yak5UWjpyuobLORrXFTo3Vhik0mPCQICLCgukRG0GfDTmEbNsIweGE2q10\n75YAEvZuIYEvhOgQSikSYsJJiGmlO6b25HefGshZGNBnxrqbWwJfKTVaa73lItvmaK1/qZTK1FrP\nc8f+hBB+7CKfGkT7tXsek1IqA3ivhadkKqUOAXnt3ZcQQgjXtfsIX2udpZRqKcwf1lova+9+hBBC\ntE9HnKmQqpTKUEo92QH7EkIIcRHKHVeqV0qt0lpPauU5c4BVWuus8x7PBDIBEhMT05YsWeJSDVVV\nVURFRbn0Wn8k7XEuaY9zSXtcyFfbZOLEiblaa3Nbnttql44zkM+Xd35wt/DaMmeXTimQev5znAO5\n8wDMZrNOT09v7W2blZ2djauv9UfSHueS9jiXtMeFAqFNWg18V2bWKKXitNblQA5nB2v7Aa9e6nsJ\nIYRwD3fM0pkKmJ1fG60GcE7VnO7cduhiUzeFEEJ4nlv68N1FKXUSKHDx5V2BU24sx9dJe5xL2uNc\n0h4X8tU26a217taWJ3pV4LeHUiqnrQMXgUDa41zSHueS9rhQILSJLCAthBABQgJfCCEChD8FvqzT\ncy5pj3NJe5xL2uNCft8mftOHL4QQomX+dIQvRItkeQ8R6Hx+PXznHP9yYLTW+gWj6/EGTc6O7qe1\n/qWhxXgJ56quk4CA/xlRSo3Geda7LGx4Toak+vsS7j59hO/8wcW5zEN54/1A5gy2LOcPbqrzvhBN\n/doZ9KmB/jvj/P4bl4rJ8/f28OnAB2bQ8JcZGpZwkHBrOHJrbIc8mlm/KNA4L9DT6tpPgcB5NLsZ\nQGv9gpz9DsAc59dUf28PXw/8OKCsyf0uRhXiLbTW85p8LB1Nw3pGgS7e6AK8yBigi1JqtIxpfLf8\nS55S6jTnZolf8vXAFxfh/Gi6xd+PWFojR/fNKm38uThvDayAo5SKo6GX4DngNaWUX38i9vVB23LO\nHr3F0bAEs2iQIQO2QEM/dSoNPyfxLV1/OUCUcnYF23IajvgDeeA2E3hOa13uvHLfVPx4YN/Xj/CX\ncraPOhWQIzkaZuk0zlgK9EFbrfWyJjNR4gwtxjss4+zvTBzO/nzx3Yyl8laf6MN8/sQr5xTEPAJg\nSlVbNLmofBkNR7XTpEtDNNV4YSJgjHwK/O78jDwg3t8zxOcDXwghRNv4epeOEEKINpLAF0KIACGB\nL4QQAUICXwghAoQEvhBCBAgJfCGECBAS+EIIESD+P/kvMUrC0SO5AAAAAElFTkSuQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Let's get a smoother response\n", "time, xc = ms.time_history(t,x)\n", "plt.plot(time,xc.T,t,x.T,'*')\n", "plt.grid(True)\n", "print('The average for this problem is known to be zero, we got', sp.average(x))" ] }, { "cell_type": "code", "execution_count": 26, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def duff_osc2(x, v, params):\n", " omega = params['omega']\n", " t = params['cur_time']\n", " return np.array([[-x-.1*x**3-.1*v+1*sin(omega*t)]])" ] }, { "cell_type": "code", "execution_count": 27, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[-0.1854827 0.74408211 1.27869989 1.42305317 1.22766696 0.6400849\n", " -0.31045409 -1.06650367 -1.39410353 -1.36873107 -0.98898229]] [[ 3.18347107e-09 -5.75875847e-10 -6.27391217e-09 -7.88667857e-08\n", " -2.08381912e-07 -1.90325802e-08 1.15335749e-09 3.95564803e-10\n", " 2.19995009e-08 1.48110614e-07 1.50174076e-07]]\n", "Constant term of FFT of signal should be zero: (-0.000670319513207+0j)\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD7CAYAAABpJS8eAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8VNXdx/HPySSZ7AtLQggkkLDIKpAAEQFBcUFwAVlU\ncEGF1rZW69PW1vr0aa3WYt2rfQnWohYEFZdWpC5BAigECAGRNUDIwhpIyE62mfP8kYkkbAmZ5c7y\ne79evDIz9965P4+T79yce+65SmuNEEII7+dndAFCCCFcQwJfCCF8hAS+EEL4CAl8IYTwERL4Qgjh\nIyTwhRDCR0jgCyGEj5DAF0IIHyGBL4QQPsLfEW+ilBqmtc6+wLL5WuvHlFLztNYLL/Y+nTp10j16\n9GhXDVVVVYSGhrZrW28k7dGStEdL0h7n8tQ22bJly0mtdee2rGt34CulJgALgOQLrDJPKTUN+FFr\n79WjRw+ysrLaVUdGRgbjxo1r17beSNqjJWmPlqQ9zuWpbaKUym/runYHvtY6XSmVe5FV5mqtl9u7\nHyGEEPZxRR9+klJqglLq1y7YlxBCiAtweuBrrZ/VWqcDHW3dP0IIIQygHDE9slLqK631ted5fR5Q\norVebjvCLz37xK1tnXkAsbGxKcuWLWtXDZWVlYSFhbVrW28k7dGStEdL0h7n8tQ2GT9+/BatdWpb\n1nXIKJ2zKaWitNalQBbQ1L+fTOPJ3RZsXwALAVJTU3V7T5p46gkXZ5H2aEnaoyVpj3P5QpvY3aVj\nG4GTavvZZBWAbajmDNuyAxcauincSMUxWDQRKo4bXYkQwsEcMUpnObD8rNdSmj2+6Nh74WbWPAsF\nmbBmPkx+wehqhBAO5JQuHeGetNYcL69l55Ey9hVVcry8hqLyWoqrannnyM0EUndm5aw3IetNLH5m\n8h7MpUfHUEx+yrjihRB2k8D3cicqalmTc4LVe4vYmFvMycozoR5m9icmwkzH0EB+HruI2eX/YPjp\n9Zip5TSBfG4Zzp9rZnHi+TUEB5hI7RHN2N6dGdOnE31jw1FKvgCE8CQS+F6oqraBdYfq+fuCDWw6\nWAJATLiZsX06Mzg+kgHxkfTtEk5EUEDLDT9dBdlrwBREsKWOm0b0pfewyew5VsGOw2V8u/8kT6/c\nDSshqVMoU4fFc+vQeLpFhxjwXymEuFQS+F7k0Klq3liby/Ith6iqs9CzUy2PXtuHa/rF0D8uovUj\n8qoiSJkDqXMgaxH+lccZGB/JwPhIpqV0A+Bo2Wky9p7gk62Hee7LHJ77MofxfTszb2wyaUkd5Khf\nCDcmge8Fck9U8urq/fxn2xGUgpsvj6dvwEnm3nrVpQXw7UvOPL7ACdu4yGDuGJHAHSMSKCypZvmW\nQyzOzOeONzIZ3C2Sh67uzYR+MRL8QrghCXwPVlZdz8ur9vHOhjwCTH7cfUUP5o7tSVxkMBkZGU4P\n3e4dQvjFtX14cFwyH2Yf4o21ucx9J4sRPTvw+I39GNI9yqn7F0JcGgl8D6S1ZummQp79Yg/lp+uZ\nOTyBR6/tQ+dwsyH1BAWYmDUykRmp3Vm2uZCX03O49bVvmTo0nicm96dDaKAhdQkhWpLA9zCFJdU8\n9uF21h8oJi2pA/930wD6xUUYXRYAASY/7kpLZMrQeF7POMCCtQfIyDnB7yf355YhXaWbRwiDyR2v\nPEHFMfSiiXyyLpsbXlrL9kNlPDN1EEvnprlN2DcXZvbnl9f3ZcVDY0joEMIj721j7jtZlFTVtb6x\nEMJpJPA9QP3qv6DzN1DxxdMMSYjii1+M5Y4RCW5/xNy3SzgfPjiKJyb1Y23OSW58eR2ZucVGlyWE\nz5LAd2dPxcAfIgnIXoQfmrv801ly6HriX000urI2M/kpHhiTxEc/GUVwoIk738jkb6v24YhZWoUQ\nl0YC341tujmDz/RoTmM76ekfDIOmw8PfG1tYOwyMj+TTh0Zz0+Vdef6rHB5aupWaeovRZQnhUyTw\n3dTyLYe4c1keVnMYQdSDfxBYasEcAeGxRpfXLmFmf16aOYTHbriMz74/yowFGzhxJF9m5xTCRSTw\n3YzWmpfT9/HLD74jLakjN/QwoVLvgwfSG6+CrfTsYFRK8eC4ZBbelcr+okrW/eNX6PwNjbNzCiGc\nSoZluhGtNX/5fA8L1uQyLaUbz0wdRIDp3TMreNF0xdd+NIhdfrVgtb1gm50TfzM8UWRobUJ4KznC\ndxNaa55csYsFa3K5Ky2RZ28bTIDJi//3PLwdBk7H6h8MwGkdyMmet3jk+QkhPIUXJ4rn0Frzf//Z\nyaJv87jvyp48ecsA/Lx97vnwLmAOx89SizaZMat6vsqtJqs4oPVthRDtIoHvBp7/Mod3NuQzb2wS\n/zu5n9uPr3cY2+ycau4qai6/h3j/CuYs2sz2Q6VGVyaEV5LAN9g/vznIq6v3c8eI7vx24mW+E/bQ\nODvn5BegyyBCprxMr4c+ITIkgLve3MTuo+VGVyeE15HAN9DHWw/x5Ipd3DCgC0/dOsi3wv48ukYF\n8+4DaQQHmLjrzU0UllQbXZIQXkUC3yCZucX86oPtXJHUkZduHyL3i7VJ6BjCv+4fQV2Dhfve2kzZ\n6XqjSxLCa0jgG6CwpJoHF28hoWMIr9+VQlCAyeiS3Erv2HBevyuFvOIqHly8hboGa+sbCSFaJYHv\nYpW1DTzwdhZWDW/eM5zIYBmVcj6jkjvxl6mDWX+gmN99/L3MvSOEA8iFVy5ktWoeWbaN/ScqeXvO\nCHp2CjW6JLd2W0o38kuqeWXVPgbGR3LPqB5GlySER5MjfBdasDaX9N3HeWJSP0b37mR0OR7hkWsa\n75H7pxW7yMorMbocITyaBL6LZOWV8NyXe5k0KI575Ui1zfz8FM/PGEJ8dDA/WZJNUUWN0SUJ4bEc\nEvhKqWEXWTZNKTVBKfVrR+zLE52qquOhpVvpFh3MM7fJ8MtLFRkcwOuzUyivqednS7bSYJGTuEK0\nh92Br5SaAHxwgWXDALTW6UDpxb4YvJXWmv/54DuKK+t47c5hRATJSdr26BcXwTNTB7Epr4RXV+83\nuhwhPJLdgW8L89wLLJ4JNF0nnwtMsHd/nmZxZj5f7ynid5P6MTA+0uhyPNqUod2YMjSeV1btY0u+\n9OcLcamc3YcfBTT/zezo5P25lbyTVfx55R7G9unM3Vd4zm0J3dmTtwwgPjqYh5dto7xGLsoS4lIo\nR4xvVkp9pbW+9jyvLwAWaK2zbV0/12qtHztrnXnAPIDY2NiUZcuWtauGyspKwsLC2rWtM1i15s8b\nazhSaeXp0cFEB7n2/Li7tYcj7T9l4c+bahjZxcSPLg9q0zbe3B7tIe1xLk9tk/Hjx2/RWqe2ZV1n\nj8MvBTrYHkcBxWevoLVeCCwESE1N1ePGjWvXjjIyMmjvts7w+poD7C/dw0szh3Dr0HiX79/d2sOR\nxgFV4ft4MT2He67pyw0D41rdxpvboz2kPc7lC23ilMNOpVSU7eF7QJLtcRKQ7oz9uZv9RRW88GUO\nEwd24ZYhXY0uxyv9ZHwy/eMieOKTnZRW1xldjhAewRGjdKYBqbafTVYBaK2zbetMAEqbnnszq1Xz\n+Ec7CDGb+NOtA2UIppMEmPz46/TBlFbX8acVu40uRwiP4IhROsu11tFa6+XNXktp9nih1jrd1nXj\n9ZZvOcSmvBIen9iPTmFmo8vxagO6RvLjq5L5MPsQGXvlPrhCtEautHWgk5W1PL1yNyN6dmB6ajej\ny/EJD13Ti14xYTz+0fdU1TYYXY4Qbk0C34H+/Nluqusa+PMU6cpxFbO/ifm3DeZIWQ2vfL3P6HKE\ncGsS+A6y4UAxH209zI+vSqZXTLjR5fiUlMRoZqR24811B9lfVGF0OUK4LQl8B7BYNU+u2EV8VDA/\nHd/L6HJ80mM3XEZIoInf/3unzJ0vxAVI4DvAe5sL2X20nMdv7Cd3rzJIxzAzv7q+L+sPFLNi+1Gj\nyxHCLUng26nsdD3PfbmXET06cOOgLkaX49PuHJnIgK4RPPXZLjmBK8R5SODb6W+r9nGquo7f39Rf\nTtQazOSnePKWgRwvr2XBmgNGlyOE25HAt0PuiUreWp/H7cO7y0yYbiIlMZrJg+NYuC6XY2VysxQh\nmpPAt8NzX+7F7O/Ho9f2NboU0cxjN1yG1QrPf7nX6FKEcCsS+O30XWEpK78/xgNjkugcLlfUupPu\nHUK498oeLM8+xM4jZUaXI4TbkMBvB11+lMB/TaZPSBVzxya1voFwuZ+O60VkcAB/XrlbhmkKYSOB\n3w5H//MkfWt38HLXLwkzO3uGadEekSEB/Pzq3ny7v5i1+04aXY4QbkEC/1I8FQN/iKTr/nfxU5p+\nhz6AP0Q2vi7czuy0ROKjgnn+y71ylC8EEviX5uHtHOo2idM6sPG5fzAMmg4Pf29sXeK8Av39eHhC\nb7YfKiO7yGJ0OUIYTgL/ElhDY9l63IJZ1aP9zWCpBXMEhMcaXZq4gKlD40nqFMrH++qwWOUoX/g2\nCfxLsHLHUQJqisnveTvqgVWQMgcqjxtdlrgIf5Mfj1zbh0OVmhXbjxhdjhCGkjOObWS1al5O3wfR\n/8vnd40FPwWTXzC6LNEGkwfF8eyn23gpfR+TBsXhb5LjHOGb5JPfRit3HGVfUSU/v6Y3Jj+ZQsGT\n+PkppvYO5ODJKv67YRssmggV8peZ8D0S+G3QdHTfOyaMGwfFGV2OaIehMSb6x0VgyZiPLsiENfON\nLkkIl5MunTZoOrr/2x1D5ejeQ41dN53x1vozL2S92fjP3wxPyP1whW+QI/xWaK35++oDJHUOlaN7\nD7Zx5EL0wGnU0DgNhpYhtcIHSeC3Yk3OCXYdLefHVyXL0b0HqzN3QJkjMFNHjQ6AhhoZUit8jgR+\nK/6ecYC4yCBuHRJvdCnCXlVF6JT7+Gnws3xmnoiWIbXCx0gf/kVsyS9h08ESfj+5P4H+8t3o8W5f\ngh9wfVwhP1seS+iQ4Yw3uiYhXEhS7CL+vvoA0SEB3D6iu9GlCAeaMjSeuMggFqyVu2IJ3yKBfwF7\njpWzak8Rc67sSUig/CHkTQJMftx3ZU8yc0vYfqjU6HKEcBm7A18pNU0pNUEp9esLLJ9v+znP3n25\n0sK1uYQEmrj7ikSjSxFOcPuI7oSb/Vm4NtfoUoRwGbsCXyk1DEBrnQ6UNj0/yzyl1AHAY36zispr\n+PS7I8xI7U5USKDR5QgnCA8K4M6RCaz8/iiFJdVGlyOES9h7hD8TaPqbOBeYcJ515mqtk21fCh7h\nnQ35NFg1c67sYXQpwonuvbIHfkrxz28PGl2KEC5hb+BHASXNnnc8zzpJF+vycTen6yws3pjPdf1j\nSewYanQ5woniIoO5eUhX3ttcSFl1fesbCOHhnH42Umv9LIBS6lql1ISzj/RtffvzAGJjY8nIyGjX\nfiorK9u9bXNfF9RTWl1PSliZQ97PKI5qD29xofYYGmzlozoLf1q6msnJvtN9J5+Pc/lCm9gb+KVA\nB9vjKKC4+UJbmJdorZfblp1zx2+t9UJgIUBqaqoeN25cuwrJyMigvds2sVo1T76whsHdgpl765Uo\n5blX1jqiPbzJxdrjqxObWHO0nKfvGYPZ3+Tawgwin49z+UKb2Nul8x5nQjwJSAdQSkXZXstqeg1I\ntj13Wxk5ReSerOL+0T09OuzFpZk3JokTFbX8e6vcIEV4N7sCX2udDaCUmgCUNj0HVjVbPkMpNQ04\n0Gy5W/rHuoPERQbJJGk+5speHekfF8HCdblys3Ph1ezuw7d1yZz9WsrFlrujnUfKWH+gmN9OvIwA\nuSOST1FK8cCYnjz6/nd8s/8kY3p3NrokIZxCks3mzW8OEhJo4vYRCUaXIgwwaXAcncICeXt9ntGl\nCOE0EvjAiYraHy60igwOMLocYQCzv4k7RySwak8RBcVyIZbwThL4wPtZhdRbNHfJNAo+bVZaIial\neGdDntGlCOEUPh/4Fqvm3Y0FjEruSHLnMKPLEQaKjQhi4qA43ssqpKq2wehyhHA4nw/8jL1FHC49\nzew0OboXcO+oRCpqGvh462GjSxHC4Xw+8Bdn5hMTbuba/nKrOwHDEqIZGB/B2+vzZIim8Do+HfiF\nJdVk5Jzg9uHdZSimABqHaN47qif7iipZf6C49Q2E8CA+nXLvbipAgQzFFC1MHhxHh9BA3pIhmsLL\n+Gzg1zZYeH9zIdf0i6VrVLDR5Qg3EhTQOEQzffdxmStfeBWfDfzPdxyjuKpOTtaK85qdloifUizO\nzDe6FCEcxmcDf0lmAYkdQxjTq5PRpQg31CUyiOv6x/J+ViE19RajyxHCIXwy8Pceq2BTXgl3jkjA\nz09mxRTnN2tkIqeq6/nvjqNGlyKEQ/hk4C/ZmE+gvx/TU7sbXYpwY6OSO9KzUyhLMguMLkUIh/C5\nwK+qbeCj7MNMGtQ4EkOIC/HzU8wamUBW/in2HCs3uhwh7OZzgf/vbUeorG1gdpoMxRStu21YNwL9\n/eTkrfAKPhX4WmsWZ+ZzWZdwhiVEG12O8ADRoYFMHhzHx9mHqZT5dYSH86nA31pYyq6j5cxOS5Rb\nGIo2mzUykao6C//eJvPrCM/mU4G/ODOf0EATtw6NN7oU4UGGJUTRLy6CxZkFMr+O8Gg+E/inqupY\nsf0oU4bFE2a2+86OwocopZidlsDuo+VsLSw1uhwh2s1nAn/5lkPUNVjlylrRLrcMiSc00CQnb4VH\n84nAt1o1Szbmk5oYzWVdIowuR3igMLM/U4bFs2L7UU5V1RldjhDt4hOB/+2Bk+QVV8vRvbDLrJGJ\n1DVY+TD7kNGlCNEuPhH4izPz6RAayMRBXYwuRXiwfnERpCRGs2RjAVarnLwVnsfrA/9YWQ3pu4uY\nntoNs7/J6HKEh5udlsDBk1VycxThkbw+8JduKsCqNbNGSHeOsN/EgXFEhwTIyVvhkbw68OstVpZt\nLmBs784kdAwxuhzhBYICTExP7c5Xu49zvLzG6HKEuCR2B75SappSaoJS6tftWe5Mq3Yf53h5rZys\nFQ51x4gELFbN+5sLjS5FiEtiV+ArpYYBaK3TgdKm521d7myLMwvoGhnE1ZfFuHK3wsv17BTK6F6d\nWLqpAIucvBUexN4j/JlA06WHucCES1zuNAdPVvHN/pPcMSIBk9zkRDjY7LQEjpTVsHpPkdGlCNFm\n9gZ+FFDS7HnHS1zuNEsy8/H3U8wcITc5EY53Tb9YYsLNLNkoJ2+F5zB8Uhml1DxgHkBsbCwZGRnt\nep/Kysoftq2zaJZurGZojIldWzLZ5aBaPUnz9hDOaY+0GCuf7j3BByu/pnOIZ41/kM/HuXyhTewN\n/FKgg+1xFHD24OTWlqO1XggsBEhNTdXjxo1rVyEZGRk0bbt8yyGq6r/jFzelMirZN29S3rw9hHPa\no/eQ06yY/zV5pq5MH3eZQ9/b2eTzcS5faBN7D0veA5Jsj5OAdAClVNTFljvb4sx8kjuHckWSy3qQ\nhA+Kjwrm6stieG9z48R8Qrg7uwJfa50NoJSaAJQ2PQdWtbLcaXYcLmNbYSmzRspNToTzzRqZyMnK\nWr7addzoUoRold19+LYumbNfS7nYcmdasjGfoAA/bkvp5srdCh81tk9n4qOCWbIxn0mD44wuR4iL\n8qwzTa0or6nnk61HuPnyrkQGBxhdjvABJj/FnSMTWH+gmAMnKo0uR4iL8qrA/zj7MKfrLXJlrXCp\nGand8fdTLN1YYHQpQlyU1wS+1prFmfkM7hbJ4G5RrW8ghIN0Djdz/cAufLDlEDX1FqPLEeKCvCbw\nc05Z2VdUyeyRcnQvXG/WyATKTtfz2fajRpcixAV5TeB/XVBPRJA/N13e1ehShA+6IqkjSZ1C5cpb\n4da8IvCLjxbwaMn/cffgYIID5SYnwvWUajx5m11Qyq4j5UaXI8R5eUXgH/nPHxmu9jLX8oHRpQgf\nNi2lG2Z/P97dJEf5wj0ZPpeOXZ6KgYZaBgEoiNz5Dux8B/zN8ITMYihcKyokkMmDu/Jx9mF+M7Ef\nYWbP/vUS3sezj/Af3o4eOA2LKajxuX8wDJoOD39vbF3CZ81KS6CqzsJ/th0xuhQhzuHZgR/eBWWO\nwGStw+IXAJZaMEdAeKzRlQkfNbR7FP3iIliyMR+t5eYowr14duADVBVByhyyh/0VUuZApcxpIoyj\nlGLWyAR2HilnW2Fp6xsIn7fjcBlvfXuQ03XOv4bD8wP/9iUw+QWqwnrC5BcanwthoFuHxhMaaGKJ\nXHkr2uCNdbk8/2UOVhf8Rej5gS+Emwkz+3PL0Hg+/e4IZdX1Rpcj3FhRRQ0rvz/KtNRuhLrgJL8E\nvhBOMGtkArUNVlZu2AqLJkKFdDWKcy3dWEi9RXP3FT1csj8JfCGcYEDXSIYmRBGS+Ty6IBPWzDe6\nJOFm6hqsLNmYz1V9OtOzU6hL9ikDhYVwhqdi+Lih9szzrDcb/8k1IsLmi53HKKqoZf5tPVy2TznC\nF8IZHt5Ow4BpnCaw8blcIyLO8vb6PBI7hnBVn84u26cEvhDOEN4F/6AIzNRTowPQDTVyjYj4wY7D\nZWTln+KutET8/Fx3K1YJfCGcpaqI6sH3MLX+SbJjpso1IuIH72zIIzjAxPSU7i7dr/ThC+Esty8h\nDEis3sL9ub3Y8MA1BBtdkzDcqao6/r3tCFOHdSMyxLW3YpUjfCGc7N5RPSitruff2w4bXYpwA+9n\nFVLbYOWeUa6/WZMEvhBONqJnB/rFRfDW+jyZX8fHWayaf2XmM7JnBy7rEuHy/UvgC+FkSinmjOrB\nnmMVZOaWGF2OMNCXO49x6NRp7h3Vw5D9S+AL4QI3D+lKdEgAb6/PM7oUYaA31uXSvUMw1w3oYsj+\nJfCFcIGgABN3jEjgy13HOHSq2uhyhAG25J8iu6CU+67sicmFQzGbk8AXwkVmpyWilOJfmXILRF/0\n5je5RAT5MyPVtUMxm5PAF8JFukYFc/2AWJZtKqSqtsHocoQLFZZU8/mOY9w5MtEls2JeiN2Br5Sa\nppSaoJT69QWWz7f9nGfvvoTwdPePTqLsdD3vZxUaXYpwoX9+exA/pQw7WdvErsBXSg0D0FqnA6VN\nz88yTyl1AMi1Z19CeIOUxGhSE6N585uDNFisRpcjXKDsdD3vby7kpsu70iUyyNBalD3jgm1H719p\nrdOVUhOAYVrrZ89aZ5rWevlF3mMeMA8gNjY2ZdmyZe2qpbKykrCwsHZt642kPVpyp/bIPt7AK1tr\n+fHlZtLijPnz3p3aw104q01W5tbxfk49fxwVRGKEyeHvP378+C1a69S2rGvvpy0KaD6wuON51km6\n0JcBgNZ6IbAQIDU1VY8bN65dhWRkZNDebb2RtEdL7tQeY62aFYfW8M1JE4/dPhqlXD9iw53aw104\no03qLVZ+s341o5IjuOfmNIe+d3s4/aSt1vpZW5dPR1vwC+HT/PwU88YkseNwOesPFBtdjnCiz7Yf\n5Vh5DQ+M6Wl0KUAbjvAvcLI1t6nfHuhgey0KaPHptW1bYuvSKQaS7CtXCO9w69B4nvsyhwVrc7my\nVyejyxFOoLXm9TUHSO4cyrg+MUaXA7Qh8G1dLhfyHtDUd5QEpAMopaK01qVAFmdO1iYDC9pfqhDe\nIyjAxJwre/DXL/ay+2g5/eJcP6+KcK7Ve4vYc6yC56Zf7tI57y/Gri4drXU2gK2rprTpObCq2fIZ\nSqlpwIFmy4XwebNHJhISaOKNtTKAzdtorXn16/3ERwVzy5CuRpfzA7uHCJzvLwCtdcrFlgshIDIk\ngJnDu/OvDfn88vq+dI2S2fK9xcaDJWQXlPLkLQMIMLnP9a3uU4kQPuj+0T3RwD/WHTS6FOFAr63e\nT6ewQEOnUTgfCXwhDNQtOoRbhnTl3U35nKysNboc4QDbD5Wybt9J7h+dRFCA48fd20MCXwiD/XR8\nL2obrLyxTvryvcFrq/cTHuTP7LQEo0s5hwS+EAZL7hzGTYO78q8N+ZyqqjO6HGGHHYfL+GLnceZc\n2ZPwINfer7YtJPCFcAM/u7oX1XUW/vmt9OV7spfS9xEe5M/9o93jQquzSeAL4Qb6xIYzcWAX3vo2\nj7LT9UaXI9rh+0NlpO8+ztwxSUQGu9/RPUjgC+E2fnZ1LypqG1gkR/ke6cX0HCKDA5hzZQ+jS7kg\nCXwh3MSArpFc1z+WN9cdlL58D7OtsJSv9xQxb2ySW/bdN5HAF8KN/M91famsa+D1tQeMLkVcghe/\nyiE6JIB7DL7BSWsk8IVwI327hHPL5V15e30eReU1Rpcj2mD9gZOsyTnBg+OSCTPw9oVtIYEvhJt5\nZEIfGiyaV1fvN7oU0QqtNX/57x66RgZx9xU9jC6nVRL4QriZHp1CmZ7anaWbCigsqTa6HHERK78/\nxvZDZTx6XV+3u6r2fCTwhXBDP7+mF0opXkzPMboUcQH1Fit//WIPfWPDmTI03uhy2kQCXwg3FBcZ\nzJxRPfh462F2HC4zuhxxHss2FZBXXM1jE/ticpP57lsjgS+Em/rJ+F5EhwTy1Ge70FobXY5oprym\nnpfS9zGiZwfG93WPu1m1hQS+EG4qMjiARyb0JjO3hPTdRUaXI5p5JX0fJdV1/H5yf0NuQt9eEvhC\nuLE7RiSQ3DmUZ1bupt5iNbocAewvquSt9XnMTO3OwPhIo8u5JBL4QrixAJMfv5vUj9yTVSzOzDe6\nHJ+nteZPK3YRHGjil9f3NbqcSyaBL4SbG983htG9OvHiVzlykxSDrd5bxJqcEzx8TW86hZmNLueS\nSeAL4eaUUvzh5gGcrrfwzMo9Rpfjs2rqLTz56S6SO4e6/RQKFyKBL4QH6BUTxgNjkvgw+xCbDpYY\nXY5Pem31fvKKq/njzQPd6sbkl8IzqxbCBz10dS/io4L53092yAlcF9t3vILX1xxg6tB4RvfuZHQ5\n7SaBL4SHCAn05/c39Wfv8QreXp9ndDk+w2rV/Paj7wk1+/O7Sf2MLscuEvhCeJDr+sdy9WUxPP9l\nDgXFMs+OKyzbXEhW/ikev7EfHT3wRG1zDgl8pdSwiyybppSaoJT6tSP2JYQvU0rx1K0DMfkpHvtw\nO1arXIFHfdLxAAALiklEQVTrNBXHqH3jehb+dwNpSR2YntLN6IrsZnfgK6UmAB9cYNkwAK11OlB6\nsS8GIUTbdI0K5neT+rEht5h3NxUYXY7X0mvmE3B4Iz+yfsBfpg72qCtqL8Tu2fq11ulKqdwLLJ4J\nfGV7nAtMALLt3acQvu724d1Zsf0Iz6zczbi+nekWHWJ0Sd7jqRhoqEUBCrjD7yt4tSv4m+EJz57i\nwtl9+FFA8zFkHZ28PyF8glKKv0wdjAZ+8+H30rXjSA9vp7LPFE7rQAC0fzAMmg4Pf29wYfYz/H5c\nSql5wDyA2NhYMjIy2vU+lZWV7d7WG0l7tOSt7TG9l4m3d53k8bfTuaFn22+e7a3tYY+mNrFYNVUH\nKphMPRYVgF9DDYdPVrBvy25gt9Fl2qXVwLcF8tlybf3yrSkFOtgeRwHFZ6+gtV4ILARITU3V48aN\na8PbnisjI4P2buuNpD1a8tb2uEprjv5rCx/uLWL2dSPaPJmXt7aHPZra5Lkv9jKwroy85Jkk3fAz\nyFpEfOVx4r2gvVoNfFsgXxKlVJTWuhR4D0i1vZwEtOVLQgjRRkop5t82mBteXsvPl25lxX29CPn3\nXJj2FoTHGl2ex1mbc4LXMvYzI+UFbpg2uPHFyS8YW5QDOWKUzjQg1fazySoArXW2bZ0JQGnTcyGE\n40SHBvLijCEcLK5i2+LHoSAT1sw3uiyPc6rGyi/e20afmHD+cPMAo8txCkeM0lkOLD/rtZRmjy/5\nLwQhxKUZtaw/B821Z4ZIZL3Z+M8LRpa4QoPFyoLttVTXKV6bNZTgQPe/IXl7yJW2QniDh7djHTiN\nWhqvBLWagrxmZIkrPL1yN3tKrDx160B6xYQbXY7TSOAL4Q3Cu+BnjiCQemoJAEstp/1CpR+/Dd7b\nXMCib/O4LtGf27zgatqLkcAXwltUFaFS51A49VOWWiewfU8OdQ0yq+bFbM4r4YlPdjCmdydm9g00\nuhynk8AXwlvcvgQmv0CvwVcQOuVlZpb9jN98uB2t5aKs8ykorubBxVvoFh3Cq3cMw+Tn+VMntMbw\nC6+EEI5369B4CkqqeeGrHLpFB/PodZ53/1VnOllZy93/3EiDVfPG3alEhrT9ojVPJoEvhJd66Ope\nHDpVzStf76dbdAgzhnc3uiS3UFnbwH1vbeZYeQ1LHkijV0yY0SW5jAS+EF5KKcXTUwZxtKyG3378\nPRHB/twwMM7osgxV22DhwcVb2HmknIV3pZCSGG10SS4lffhCeLEAkx+vz05hSPcoHlq6lVW7jxtd\nkmFqGyz8dEk26/ad5Jmpg7imn++NYJLAF8LLhZr9WTRnOP3iInhwcTZrc04YXZLLNYV9+u4i/nTr\nQGak+mb3lgS+ED4gIiiAd+4bQa+YMOa+k8W2ogajS3KZs8P+rrREo0syjAS+ED4iKiSQxQ+MpG+X\ncP62tZZPth42uiSnK6+p555/bpKwt5HAF8KHdAgN5N25afSJ9uOR97bx1rcHjS7JaY6V1TDj9Q1s\nyT/FSzOH+HzYg4zSEcLnhJn9+UVKEMsPh/OHT3eRX1LN727sh7/Je47/9hwr5/63siitruOf9w5n\nTO/ORpfkFrzn/7AQos0CTYq/zxrG/aN7sujbPOa8tZmy6nqjy3KIFduPMOW19dRbrLz3oysk7JuR\nwBfCR/mb/Pjfyf159rbBZOYWc8tr37DzSJnRZbVbg8XKM//dzc/e3Ur/rhGseGh0m+8A5isk8IXw\ncTOGd2fp3DSq6yxMeW09i7496HHz7xQUVzNzYSYL1uQya2QCS+emERMRZHRZbkcCXwhBao8OfP7I\nWMb07sQfP93F/W9nUVReY3RZrdJas3zLIW58ZR05xyp4ceblPD1lEIH+Em3nI60ihAAaR/D8455U\n/njzAL7Zf5JrXljD4sx8rFb3PNovKK7mvrc288sPvqN/1wj++8gYpgz17vns7SWBL4T4gVKKe0b1\n4ItHxjIoPpInPtnB9AUb+K6w1OjSoOIYLJpIzakjvJy+jwkvrmHTwRKemNSPpXPT6BYdYnSFbk8C\nXwhxjp6dQlnywEien345eSeruOW1b/npkmxyT1SeWckWwFS4Zn4eS8Z8dP4GPvvbI7yYnsN1/WNZ\n9T/jeGBMkk/MZe8IMg5fCHFeSiluS+nGdQNieWPdQf6xLpfPdx5j0qA47h/dk8u/exYKMmHNfJj8\ngtPq0H+KQVlqabqt+G3WL7gt6AvINUOk3KD9UkjgCyEuKjwogEev7cNdaYksWHOAX20ei3lvszH7\nWW82/vM3wxOOCWCtNbuPVvB+ViEZDa/wiOVtbvDfQhC1aP9gVL/JcN3TDtmXL5HAF0K0SedwM09M\n7k9l2lYKP/gV3Y6vIog6TutAvgsbw/G0JxhaXE33DsEodeldLDX1Fr4rLGX13hN8vuMoecXVBJr8\nmDT4Mq609sSckwmmIJSlFswRcoP2dpDAF0JckrBO3enVLQ5d1IDVz4zZUkdhtT+/+uwofHaU+Khg\nBsZHkNw5jOTOYcRGBBEe5E94kD9KKU7XWahpsHCiopaC4mryS6rYeaScHYfLqLdoTH6KUckd+dFV\nyVzXP5aOYWZYVgYpcyB1DmQtgkrfndffHhL4QohLV1WESpmDsgXwtMrjDBk/lg25xWzMLWHPsXJW\n7S6ioQ1DOiODA+gdE8b9o5NITYwmtUc0USGBLVe6fcmZx048X+DtJPCFEJfurABWQG+gd2w4d1/R\nA4B6i5WCkmqKK+uoqKmnoqZxDv6gAD/MASY6hgaS2CHUZ24g7g4cEvhKqWFa6+wLLJuvtX5MKTVP\na73QEfsTQri/AJOfrVvH6EpEE7vH4SulJgAfXGSVeUqpA0CuvfsSQgjRfnYf4Wut05VSFwvzuVrr\n5fbuRwghhH1ccaVtklJqglLq1y7YlxBCiAtQjpgGVSn1ldb62lbWmQ98pbVOP+v1ecA8gNjY2JRl\ny5a1q4bKykrCwsLata03kvZoSdqjJWmPc3lqm4wfP36L1jq1Leu22qVjC+Sz5Z4d3BfZtsTWpVMM\nJJ29ju1E7kKA1NRUPW7cuNbe9rwyMjJo77beSNqjJWmPlqQ9zuULbdJq4LdnZI1SKkprXQpkceZk\nbTKw4FLfSwghhGM4YpTONCDV9rPJKgDbUM0ZtmUHLjR0UwghhPM5YpTOcmD5Wa+lNHssY++FEMIN\nOOSkraMopU4A+e3cvBNw0oHleDppj5akPVqS9jiXp7ZJota6TZe3uVXg20MpldXWM9W+QNqjJWmP\nlqQ9zuULbSJ3vBJCCB8hgS+EED7CmwJfTg63JO3RkrRHS9Ie5/L6NvGaPnwhhBAX501H+EJclMzn\nJHydx98AxXZRVykwTGv9rNH1uINm02Eka60fM7QYN2GbxvtawOc/I0qpYdimOZGZbFtkSJK3Xzfk\n0Uf4tg8utnl9Spue+zJbsKXbPrhJtudCNPdbW9An+frvjO2/v2lusFxvbw+PDnxgJo3fzNA4Z4+E\nW+ORW1M75HKeCet8je2ObK1O9ucLbEezmwG01s/KdCcAzLf9TPL29vD0wI8CSpo972hUIe5Ca72w\n2Z+lw2icwM7XdTC6ADcyHOiolBom5zR+mO8rVyl1ipZZ4pU8PfDFBdj+NM329iOW1sjR/XkVN30u\nzpr00OcopaJo7CV4BnhDKeXVfxF7+knbUs4cvUXROOe+aDRBTtgCjf3USTR+TjrYvgB8+UuwmDNT\nlpfSeMTvyydu5wHPaK1LbbdqnYYXn9j39CP89zjTR50EyJEcjaN0mkYs+fpJW6318mYjUaIMLcY9\nLOfM70wUtv588cOIpdJWV/RgHn/hlW0IYi4+MKSqLWwB/wGN/ZEdgOnSpSGaa7oTHTBc/gr84fqM\nXKCDt2eIxwe+EEKItvH0Lh0hhBBtJIEvhBA+QgJfCCF8hAS+EEL4CAl8IYTwERL4QgjhIyTwhRDC\nR/w/82zqY6D/q3UAAAAASUVORK5CYII=\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "t, x, e, amps, phases = ms.hb_time(duff_osc2, sp.array([[0,1,-1]]), .7, num_harmonics=5)\n", "\n", "print(x,e)\n", "print('Constant term of FFT of signal should be zero: ', ms.fftp.fft(x)[0,0])\n", "time, xc = ms.time_history(t,x)\n", "\n", "plt.plot(time, xc.T, t, x.T, '*')\n", "plt.grid(True)" ] }, { "cell_type": "code", "execution_count": 28, "metadata": { "scrolled": true }, "outputs": [ { "data": { "text/plain": [ "[]" ] }, "execution_count": 28, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD7CAYAAAB68m/qAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHatJREFUeJzt3Xl8FGWeP/DPk6tzp3MBgSRAoggBBJKGgOOM7kx01XVG\nx8HBA8IREkSdg5cuOrs7u+OoP0d2dtz56awa7nAJwhzeOwYddTwgnXBKACFAEgg5Omlyd9Ldz/6R\nioQY0oF0d3VVf96vF6+u6irT38dOfXioeuopIaUEERFpX4DaBRARkXsw0ImIdIKBTkSkEwx0IiKd\nYKATEekEA52ISCcY6EREOsFAJyLSCQY6EZFOBHnzwxISEuS4ceO8+ZFERJpXWlraIKVMdLWfVwN9\n3LhxMJvN3vxIIiLNE0KcGcp+POVCRKQTDHQiIp1goBMR6QQDnYhIJxjoREQ6wUAnItIJBjoRkYc5\nnN55MhwDnYjIg46ca8YtL3yEL89d8PhnMdCJiDzkYLUV96/+Ah1dDoSHeP4+ziEFuhAic5Btc4UQ\nOUKIle4ri4hI20rPNOHB1XsQFRqEHcvmYHxChMc/02WgCyFyALx+mW2ZACClLAZgHSz4iYj8xd5T\njchduwcJUQbsWDYHKXHhXvlcl4GuhHXFZTbPA2BVlisA5LipLiIiTfr8pAUL1+3FqJhQbC+YjdHG\nMK999nBP6hgBNPZZjx/mzyMi0qxPTzQgb2MJUuPCsWXpbCRGGbz6+R6/KCqEKBBCmIUQ5vr6ek9/\nHBGRKj4+Xo8lG0owLj4CW/O9H+bA8APdCiBOWTYCsPTfQUpZKKU0SSlNiYkup/MlItKcvx2rw9Ii\nM9ISI7E1fzYSIr0f5sBVBroQwqgsbgeQpiynASh2R1FERFrxwdFaFBSVYsLISGzLz0ZcRIhqtQxl\nlMtcACbltdduAJBSlin75ACw9q4TEfmD94/UYtmmUkxMisKWvNkwhqsX5sAQLopKKXcC2Nnvvaw+\ny4UeqIuIyKe9d/g8Ht1ahiljYrBxySzEhAWrXRLvFCUiulJvH6zBI1vLcH1yDIryfCPMAS8/U5SI\nSOveOHAOK7bvR2aqEesXz0KkwXdilD10IqIh+vO+s/j5a/uQNTYWG3wszAEGOhHRkOwsrcaKHfuR\nPT4eGxbPRISPhTnAUy5ERC69bq7Cyl0H8a30BKzONSEsJFDtkgbEQCciGsTO0mqs3HUQN17TE+ah\nwb4Z5gBPuRARXdaf9lXjn3ce+Lpn7sthDjDQiYgG9Jf9Z/HYjgOYkxaviTAHGOhERN/QOzRx1vg4\nrF0402fPmffHQCci6uOtg+fw89f2wTQuDusWaSfMAQY6EdHX3jlUg5+9th+msXFYv2imV54D6k4M\ndCIiAO8drsFPtu3DjBQj1vnoOHNXGOhE5Pf+98vzeHTrPkxLjsH6xTN97g7QoWKgE5Ffe/9I7SWz\nJkaF+sZEW1eDgU5Efmt3eS0e3lKKjKRoFOVpO8wBBjoR+akPj9Zh+eYyTEqKRlFeNqI1HuYAA52I\n/NDfjtVh2aZSTBgViU1Lsn1mPvPhYqATkV/5+Hg9CjaV4poRkdicl42YcH2EOcBAJyI/8umJBuQX\nmZGeGIktS7NVfwaouzHQicgv7D3ViKUbzRgbH44tS7MRG6GvMAcY6ETkB/ZXWbFkQwmSjKHYsnQ2\n4nQY5gADnYh07vDZC8hduwdxESHYunQ2EqMMapfkMQx0ItKtY+dbsGDtHkQagrA1PxujYkLVLsmj\nGOhEpEsV9a14cM0eBAcGYGv+bCTHhqtdkscx0IlIdyot7Xhg9R5IKbE1PxvjEiLULskrtDkDDRHR\nZZy1duCBNV+g0+7AtvzZuGZElNoleQ176ESkG7XNnXhw9Re40N6NTUuyMSkpWu2SvIo9dCLShYZW\nGx5cswd1LTZsysvG1OQYtUvyOvbQiUjzrO1dmL9mD6qb2rFu0UxkjY1VuyRVMNCJSNOaO7uRu24v\nKurbULjAhNlp8WqXpBoGOhFpVpvNjsXrS3DkXDNenp+J70xIVLskVfEcOhFpUme3A/lFZuyvsuKl\n+2fge5NGql2S6thDJyLN6XY48ejWMnx20oLf3ns9bp+apHZJPsFlD10IMReAFUCmlHLVINvTpJSF\n7i+RiOgih1PisR0HUFxeh6fvmowfzkhWuySfMWgPXQiRCQBSymIA1t71ftsrlO0V/bcTEbmTlBL/\n9ufDeOPAOay87TosmDNO7ZJ8iqtTLvPQ0/sGgAoAOQPs87zymialLHNXYUREfUkp8Zt3j2Lb3kos\nvzkdD998jdol+RxXgW4E0Nhn/ZLxQEqAVwghmvrtR0TkVn/48ARe/bgCC2aPxcp/vE7tcnzSsC6K\nCiGM6OnBPwdgtRAibYB9CoQQZiGEub6+fjgfR0R+av2np/Dbvx7HPTPG4KkfTIYQQu2SfJKrQLcC\niFOWjQAs/bYXAHhOuViaD2Bu/x8gpSyUUpqklKbERP8eI0pEV+51cxWeevMIbs0YiVVzr0dAAMP8\nclwF+nYAvb3uNADFwNc980tIKXfi4vl2IqJhe/dQDZ7YdRA3XpOAFx+YgaBAjrQezKDDFqWUZUII\nkxAiB4C1z0XP3QCypJSrhBArhRAVAOI4bJGI3OVvx+rw09f2YXqKEYW5WTAEBapdks9zOQ59oJCW\nUmb1Wf7G2HQiouHYe6oRD20uxbUjorB+8SyEh/Cm9qHgv1+IyKccqr6AvA0lGG0MQ1HeLMSEBatd\nkmYw0InIZ3xV24LcdXsQHRaMzXnZSIg0qF2SpjDQicgnVDe1Y8HavQgMCMDmpdkYbQxTuyTNYaAT\nkeosrTbkrt2Lti47NuXNwng/eaizuzHQiUhVrTY7Fq0vwVlrB9Ytmul3zwF1J146JiLVdHY7UFBk\nxpGaZqzOzcLMcXGu/yO6LPbQiUgVDqfEz1/bj89OWvCfc6/HdyfyARXDxUAnIq/rmQb3EN778jx+\neWcG7snknObuwEAnIq/77V+PYdveKjzyD+nIu3G82uXoBgOdiLxqzScV+MOHJ3H/rFQ8fiunwXUn\nBjoRec2u0mo883Y5bp8yCs/cPYXT4LoZA52IvGJ3eS1W7jqIG9Lj8d/3TUcgp8F1OwY6EXlcyelG\nPLylDJNHR6Mw18SZEz2EgU5EHnX0fDOWbCjBmNgwrF80E5EG3v7iKQx0IvKYs9YOLFy3FxEhQdiU\nl414TrblUfyrkog8wtrehYXr9qK9y4GdD92AMZxsy+PYQycit+vsdiBvoxmVlnaszjXhulFRapfk\nF9hDJyK3cjglfrptH8oqm/CHBzIxOy1e7ZL8BnvoROQ2Ukr8+18O469HavEfd2bgjqlJapfkVxjo\nROQ2L31wAlv2VOKhm9Kx6Fu8pd/bGOhE5BY7SqrwX+8fxz0zxuCJ23hLvxoY6EQ0bB8crcUv/nQI\n3742Ac/PvZ639KuEgU5Ew7K/yopHtuxDRlI0Xp6fheBAxopa+H+eiK5aRX0rlmwoQWKUAet4F6jq\nGOhEdFXqWjqxcP1eAMDGJbOQGMW7QNXGQCeiK9Zqs2PJhhI0tHRh3aKZGJ8QoXZJBN5YRERXyO5w\n4pEtZSivacGaXBOmpxjVLokU7KET0ZBJKfHLv3yJj47X45m7p+AfJo5QuyTqg4FOREP2ykcV2La3\nEg/fnI77Z6WqXQ71w0AnoiF588A5PP/eUXx/2mg+C9RHMdCJyKWS0414bMcBzBoXh9/eez0C+Pg4\nn8RAJ6JBVdS3Ir/IjOTYMLy6IIuPj/NhLke5CCHmArACyJRSrhpgeyaANACQUu50e4VEpBpLqw2L\n1pcgUAhsWDwLsREhapdEgxi0h66ENaSUxQCsvev9/EIJ8rTLbCciDersdmBpkRm1zZ1YvdCE1Phw\ntUsiF1ydcpmHnt45AFQAyOm7Uem9lwCAlHKVlLLM7RUSkdc5nRIrtu/H/iorfn/fdGSmxqpdEg2B\nq0A3Amjss97/0SMzAcQLITKFECvdWhkRqea5d8vx7uHz+Nc7JuG2KXxIhVa446KopbdnrvTYLyGE\nKBBCmIUQ5vr6ejd8HBF5UtHnp7H6k1NYdMM45N3Ih1RoiatAtwKIU5aNACz9tlvQcyqmd9+Z/X+A\nlLJQSmmSUpoSExOHUysReVjxkVr86o0vkTNpJH55ZwbnNdcYV4G+HcoIFuW1GACEEL2TN+zss90I\n5Xw6EWnP4bMX8JNt+zBlTAz+//3TEcix5pozaKD3OZWSA8Da56LnbmV7BXpGv8wFEM9hi0TaVNvc\niaUbzYgND8aahSaEh3DePi1y+a1JKQsHeC9rgO0McyIN6uhyIL/IjObObux86AaMiApVuyS6Svxr\nmMiPOZ0Sj72+H4fOXsDqBSZkjI5WuyQaBt76T+THfvf+cbxzqGd4Yk7GSLXLoWFioBP5qT+WVeOl\nD0/gvpkpHJ6oEwx0Ij9UcroRT+46hDlp8fj1XVM4PFEnGOhEfqbS0o5lm0oxJjYML8/PREgQY0Av\n+E0S+ZHmzm7kbSyBwymxdqEJxnDOnqgnDHQiP2F3OPHo1n041dCGl+dnIi0xUu2SyM04bJHITzz9\n1hF8fLwev7lnKm5IT1C7HPIA9tCJ/MDGz05j4+dnUPCdNNzHhzvrFgOdSOc+Ol6Pp97smXDridsm\nql0OeRADnUjHKupb8ejWMkwYGYXf38cJt/SOgU6kUy2d3cgvMiMoQGB1rgkRBl4y0zt+w0Q61PsI\nudOWdmzOy0ZKHJ8H6g/YQyfSoReKj6O4vA7/fmcG5qT3f3Ik6RUDnUhn3j5Ygxc/OIF5phTkzhmr\ndjnkRQx0Ih05cq4Zj79+AJmpRvz67smco8XPMNCJdKKxrQsFm8yIDgvCK/OzYAgKVLsk8jJeFCXS\ngW6HE49sKUNdiw07ls3BiGg+dcgfsYdOpAPPvl2OzysseO6HUzE9xej6PyBdYqATadyOkips+Ow0\n8m4cjx9lJatdDqmIgU6kYWWVTfi3Px/Gjdck4Be387Z+f8dAJ9Ko2uZOPLSpFKNiQvHSAzMQFMjD\n2d/xN4BIgzq7HSjYVIpWmx2rc/mgCurBUS5EGiOlxL/+6TAOVFnxyvwsXDcqSu2SyEewh06kMes/\nPY1dZdX42feuxW1TRqldDvkQBjqRhnx6ogHPvlOOWzNG4mffu1btcsjHMNCJNKLS0o5HtpYhPTEC\nv5s3HQGc25z6YaATaUCbzY78IjOkBFbnmhDJuc1pAAx0Ih/ndEo8/voBfFXXgpcemIGx8RFql0Q+\nioFO5ONe+vAE3j18Hv9yxyR8+9pEtcshH8ZAJ/Jhf/3yPH73/nHcM2MM8m4cr3Y55OMY6EQ+6nht\nC1Zs349pyTH4f/dM5dzm5JLLQBdCzBVC5AghVrrYb9DtRDR0F9q7UVBkRlhIEF5ZkIXQYM5tTq4N\nGuhCiEwAkFIWA7D2rg+wXw6AW9xfHpH/sTuceHRbGc5aO/DqgkwkxYSpXRJphKse+jwAVmW5AkCO\nZ8sholX/ewyffNWAp++agqyxcWqXQxriKtCNABr7rH/j8eFCiEylB09Ew/TnfWdR+HEFcueMxX2z\nUtUuhzTGHRdFB+1CCCEKhBBmIYS5vr7eDR9HpE8Hq614YtdBZI+Pwy/vzFC7HNIgV4FuxcXANgKw\n9N04lN65lLJQSmmSUpoSEzmGlmgg9S02LNtUioRIA/7nwUwEc25zugqu7h/eDsCkLKcBKAYAIYRR\nSmkFkCaESENP6McpAV/msWqJdKjL7sTyzaVoau/CruU3ID7SoHZJpFGDdgN6w1kZxWLtE9a7le07\npZQ7lff4ZFqiKySlxH+8cRjmM034z7nTMHl0jNolkYa5nOFHSlk4wHtZA+zzjf2IaHCb91Ri294q\nPHxzOr4/bbTa5ZDG8UQdkUr2VFjw1Btf4rsTR+CxW69TuxzSAQY6kQqqm9qxfEsZUuPD8d/3TUcg\n5zYnN2CgE3lZm82OpRvN6HY4sSbXhOjQYLVLIp1goBN5Ue/c5sdrW/DSA5lIS4xUuyTSEQY6kRe9\n+MHFuc1vmsD7Msi9GOhEXvLe4Rq8UHwcP8pM5tzm5BEMdCIvKK9pxortBzAj1YhnfziFc5uTRzDQ\niTzM0mrD0o1mxIQF49X5nNucPIePDifyoC67E8u3lKGh1YYdy+ZgRHSo2iWRjjHQiTzoqTe/xN5T\njfj9fdMxLYWzY5Bn8ZQLkYds+PQUtuypxPKb03HX9DFql0N+gIFO5AG7y2vx67eO4NaMkXict/WT\nlzDQidzs8NkL+Mm2fZg8Ooa39ZNXMdCJ3KjmQgfyNpbAGBaMtQtNCA/hZSryHv62EblJq82OvA1m\ntNkc2LmcI1rI+xjoRG5gdzjx0237cKy2BesWzcTEUdFql0R+iKdciNzgmbfL8cHROjz1g8mco4VU\nw0AnGqbVH1dgw2enkf/t8Zg/e6za5ZAfY6ATDcOu0mo8+0457pg6Ck/ePkntcsjPMdCJrtKHR+uw\nctdB3JAejxfmcXgiqY+BTnQVSs80YfmWUkxKikJhrgmGIE64RepjoBNdoeO1LViyoQSjokOxYfEs\nRBo4WIx8AwOd6AqcsbQhd+1ehAQFYFNeNhIiDWqXRPQ1BjrREFU1tuP+wi9gszuwKW8WUuLC1S6J\n6BL8tyLREFQ3teP+1V+grcuBrfnZvHGIfBJ76EQu1FzowAOr9+BCRzc252Vj8ugYtUsiGhADnWgQ\ntc2duL/wCzS1dWFTXjamJjPMyXfxlAvRZVRa2jF/7R5YWm0oysvGdD5xiHwcA51oAEfPNyN37V50\nOZzYkj+bYU6awEAn6qessgmL15cgNDgAO5bNwYSRUWqXRDQkDHSiPj75qh4FRaUYEW3A5rxsDk0k\nTWGgEyl2lVbjyT8eRHpiJIryZmFEFB9QQdrCQCe/53RKvFB8HC9+cAI3pMfj5flZiAkLVrssoivm\nMtCFEHMBWAFkSilXDbC9QFlMl1I+4eb6iDyqo8uBf955AG8drME8Uwqe+eEUBAdyNC9p06C/uUKI\nTACQUhYDsPau99meA6BYSlkIIE1ZJ9KESks77nn5M7x9qAZP3j4Rv/nRVIY5aZqr39556OmdA0AF\ngP6BndbnvQplncjnfXisDne++AnOWTuwbtFMPHRTOoTgfOakba5OuRgBNPZZj++7UemZ98oEsN1N\ndRF5RLfDiRfeP46XPzqJiaOi8er8LKTGcyQL6YNbLooqp2LKpJRlA2wrAFAAAKmpqe74OKKrcsbS\nhp++th8HqqyYZ0rBr34wGWEhfDAF6YerQLcCiFOWjQAsl9kv53IXRJVefCEAmEwmeTVFEg2HlBI7\nzFV4+q1yBAjgfx7MxB1Tk9Qui8jtXAX6dgAmZTkNQDEACCGMUkqrslzQO/pFCJGjXEAl8glnrR14\nctdBfPJVA2anxeG/fjwdY4xhapdF5BGDXhTtPYWijF6x9jmlsrvP+88LIU4KIZo8WinRFbA7nFj3\n91P4xxc+RumZJjx912RsXTqbYU665vIcer8Ln73vZSmvxQBiPVAX0VX7/KQFv3rjSxyrbcF3JiTi\n2bun8BZ+8gu8U5R0o+ZCB559uxxvHaxBcmwYChdk4ZaMkRyOSH6DgU6a19jWhVc+OomNn50GAKzI\nmYBlN6UhNJgjWMi/MNBJsy50dGPNJxVY9/dT6Oh24O4ZY7AiZwJPr5DfYqCT5jS02lD0+Rls+PQU\nmjvt+KepSVhxy7W4ZgTnLSf/xkAnzThR14I1n5zCH/edRZfdiZxJI7Hilmv50GYiBQOdfJrd4cTf\njtVjy54z+PBYPQxBAZiblYy8G8cjPTFS7fKIfAoDnXxSpaUdO8xVeL20CrXNNiREGrAiZwLmz05F\nfKRB7fKIfBIDnXxGU1sX3j18Hm8cOIsvKhoRIICbrxuBX9+Vgu9OHMGpbYlcYKCTqqztXdhdXoc3\nD57D379qgN0pkZYQgRU5E3CvKRmjeWcn0ZAx0MnrKupbsbu8Du+X16L0TBMcTokxxjAs/XYavj8t\nCRlJ0bwZiOgqMNDJ4xrbuvD5SQs+PdmAz0404LSlHQAwcVQUlt+UjpyMkZiWHMMQJxomBjq5XVNb\nF8oqm5QQt6C8phkAEGkIQvb4OCz+1nh8b9IIJMfyBiAid2Kg07A4nBLHzregrLIJZZVN2FdpxamG\nNgBASGAAssbG4vFbJ2BOegKuT47hhU0iD2Kg05C1d9lx9HwLymuaceRcM8prmnH0fAvauxwAgPiI\nEMxIjcW9pmTMSInFjFQj51Mh8iIGOn1DZ7cDpy1tqKhvw8m6Vhyt7QnxUw1tkMozp6JCg5CRFI0f\nm1IwPcWIzNRYpMSF8Tw4kYoY6H6qy+7EOWsHqps6cKZRCe/6Vpysb0V1U8fXwQ0AybFhyEiKxg+m\njUZGUjQmJUUjOZbhTeRrGOg61dntQG1zJ85ZO1Hd1I6qpg5UN7ajuqkDVU3tON/ceUlohwYHYHxC\nJKYlG3HPjGSkJUYgPTESaYkRCA/hrwmRFvBI1RApJdq7HGhs64KlrQt1zZ2obbGh9kInaps7cb65\nE3XNNpxv7sSFju5L/lshgFHRoUiJDcectHgkx4UjOTYMKbHhSIkLw+iYMAQEsMdNpGUMdJVIKdFq\ns+NCRzeaO5TXzm5caO+Gpa0LjW02WFq7lOWePw2tNtjszm/8rMAAgcRIA0ZGGzA2PhyzxsdhZLQB\nI6NDMSqmJ8RHG8MQEsQRJkR6xkC/AlJKdDmcaLc50NZlR3uXA602e591O9psjq9f22z2i6HdeWlw\nN3d0wykv/1mhwQGIjzAgPjIE8ZEhmDAyCvGRIYiL6PmTEBmChEgDRkWHIj7SgED2ron8niYC3dre\nhZZOO7odTjicEt0OCbvTCbtTwu6QsDuUZacT3Q6p7OOEvXfZ2bPc7XDCZneis9sBm90Jm/L69brd\nCZvdgc7unldb98X9O7sdaO9ywD5YCvcTHhKISEMQYsKCER0WjMQoA9ITIxAdFtzzXqjyGhaE6NDg\nr9+PjwzheWsiumKaSI3n3jmK7eYqt/08IQBDUABCgwNhCAqAIUh5DQ5AaFAgQoMDEBMW3G+fAEQY\nghBhCEJ4SGDPckgQwg2BPa9KePeuhwUH8pw0EXmVJgJ9rikZpnGxCA4MQFCgQFCAQFBA73LPa3Cg\nQGBAAIICBIIDAxAY0PNeUGCAsr9AiBLewYGCQ+6ISHc0Eegzx8Vh5rg4tcsgIvJpHPZARKQTDHQi\nIp1goBMR6QQDnYhIJxjoREQ6wUAnItIJBjoRkU4w0ImIdEJIOfS5SYb9YULUAzjT560EAA1eK8A7\n9NYmvbUH0F+b9NYeQH9tGm57xkopE13t5NVA/8aHC2GWUppUK8AD9NYmvbUH0F+b9NYeQH9t8lZ7\neMqFiEgnGOhERDqhdqAXqvz5nqC3NumtPYD+2qS39gD6a5NX2qPqOXQiInIftXvo5GOEEJmDbJsr\nhMgRQqz0Zk3D4aI9zyuvBd6riMhzvBborsJAa2ExhPZoLiyEEDkAXr/MtkwAkFIWA7AOFpS+YrD2\nKAqEECcBVHippGETQhQof56/zHatHUeu2qPJ40j54/XvyCuB7ioMtBYWQ6xXc2GhtOdy9c4DYFWW\nKwDkeKWoYXDRHgDIl1KmK/v5POUvqGIpZSGANGW973atHUeDtkehqeNIacO9yneQ6e2s81YP3VUY\naC0shlKvpsJiCIwAGvusx6tViBulaak3CyANF3/XKpT1vrR2HLlqD6Cx40hKWSylXKaspkkpy/rt\n4tHvyFuB7ioMtBYWQ6lXa2Hhd6SUq5SgiL9M79CnSCkLld4sAGQCMPfbRVPH0RDaA2j0OFLqXTbA\nJo9+R7wo6iFaC4shsALofbCrEYBFxVqGTTlvO1dZtWDg3qFPUv6ZXjZA70+TBmuPVo8jKeUqAMuE\nEEZvfq63At1VGGgtLAatV8th0V+fX8jtuNiONACa+Cdwf33aY8bFNqRj4N6hr8qRUj4xwPtaO456\nDdgeLR5HQoi+580rAPS/mOvR78hbgT5gGGg4LFy1R5NhoRw8pj4HEQDsBoDe3pPSS7JqoXc4hPb8\nWNl2UgvtAXpCTun99X4XWj6OXLVHi8dRDi4N7ArAe9+R124sUoYdVaDnQkGh8l6plDLrctt92RDb\n06hsX6VepaQXfYZhNqInNO6VUhZr9Ti6gvZo5jhSgvvHympW7wVSb31HvFOUiEgneFGUiEgnGOhE\nRDrBQCci0gkGOhGRTjDQiYh0goFORKQTDHQiIp34P0TMHWgUqQVNAAAAAElFTkSuQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "omega = sp.linspace(0.1,3,200)+1/200\n", "amp = sp.zeros_like(omega)\n", "x = sp.array([[0,-1,1,0,0]])\n", "for i, freq in enumerate(omega):\n", " #print(i,freq,x)\n", " try:\n", " t, x, e, amps, phases = ms.hb_time(duff_osc2, x, freq)#, f_tol = 1e-10)#, callback = resid)\n", " amp[i]=amps[0]\n", " except:\n", " amp[i] = sp.nan \n", "plt.plot(omega, amp)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The break is an indicative of a break in the branch and is actually a result of the `solution` being unstable. Not the system, but the solution. By that we mean that while this is considered a solution, it isn't one that will actually continue in a real situation and another solution will necessarily be found. \n", "\n", "A simple solution is to change the starting guess to be away from the solution and see if it finds another one. Indeed that happens. \n", "\n" ] }, { "cell_type": "code", "execution_count": 29, "metadata": { "scrolled": true }, "outputs": [ { "data": { "text/plain": [ "[]" ] }, "execution_count": 29, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXEAAAD7CAYAAACc26SuAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAG4lJREFUeJzt3Xl8VNXdBvDnZF/IvhBjSELYwiaQBSRaWzVa26pURBEE\nwYUAWtFXq9jWvlVri/Jaa7EqxqUWArIEtC6VatBapKgkk4UtJJCQBcjOZJ/JLOf9I6NNEchkMjN3\n7szz/XzyYSYzeH/HyX04Ofece4SUEkREpE5eShdARES2Y4gTEakYQ5yISMUY4kREKsYQJyJSMYY4\nEZGKMcSJiFSMIU5EpGIMcSIiFfNx9AGio6NlcnKyow9DRORWioqKWqSUMYO9z+EhnpycjMLCQkcf\nhojIrQghaqx5H4dTiIhUjCFORKRiDHEiIhVjiBMRqRhDnIhIxRjiREQqxhAnInIAk9k5u6YxxImI\n7Oykthc3v/JvvF96yuHHcvhiHyIiT/J5RTMe3FIMg0nC11s4/HgMcSIiOzCZJdbtrsS6TysxPjYE\nryxKQ0rMCIcflyFORDRMbd19eGBLMfZUtmDujIvx9E1TEOTnnHhliBMRDUNx7Rnct0mDlq4+/P6m\nqVgwcxSEcPwwyjcY4kRENpBSYsO+Gjz94WGMDA3AjpVZmJoQ5vQ6GOJEREPUrTfisZ0H8H7pKVyV\nGovnb52G8CA/RWphiBMRDcGxpk6syNOgqrkLj/xwAlZ+fwy8vJw3fHI2hjgRkZU+KDuF1fllCPD1\nxsa7Z+GysdFKl8QQJyIaTJ/RjDUfHcFf9p5AelIEXlqYhriwAKXLAsAQJyK6oNPtvbhvkwaaWi3u\numw0fvHjVPh6u85id4Y4EdF57D3WglVvF0NnMOHPC2fg+kvilS7pOxjiRERnMZslXvn8OP7w8VGk\nxIzA+kVpGBsbonRZ58QQJyIaoL3HgIe3l6DgSBNumBaPZ+ZORbC/60al61ZGRORkh061Y2WeBqfb\ne/HEDZOwJCvZqasvbcEQJyICsK2wDr9+9yAigvywJWc20pMilC7JKgxxIvJoOoMJT7x3CFv21yFr\nTBTWLZiB6BH+SpdlNYY4EXmsurYerNxUhIMnO3DvD8bg4WsnwFvB1Ze2YIgTkUf67GgTHtxSArOU\neP2ODGRPGql0STZhiBORRzGZJf60uxIvflqJ1LhQrF+UhqSoYKXLshlDnIg8xpnuPjywtQT/qmjG\nvPQEPP3TKQjw9Va6rGGxKsSFEGlSSs2A5/MAaAGkSSnXOqo4IiJ7KavXYmWeBs2dekU2b3CUQW8A\nIITIBrB9wPM0AJBSFgDQfvOciMgVSSmx+atazHtlHwBg+4rZWDgr0S0CHLCiJy6lLBBCVA341nwA\nn1geVwHIBqD5zl8kIlKYzmDC4+8eRH5RPa4YH4M/zZ+OiGBlNm9wFFvGxMMBtA14HmWnWoiI7Ka2\ntQcr8opw+HQHVl09Dg9cPU510wet4ZALm0KIHAA5AJCYmOiIQxARnden5Y14cEsJAOAvSzNxZWqs\nwhU5ji03xdUCiLQ8DgfQevYbpJS5UsoMKWVGTEzMcOojIrKaySzx/CcVuOutQiREBOGD+7/n1gEO\n2NYT3wogw/I4BUCB/cohIrKNO04ftMagIW6ZTpghhJgnpcyXUmqEEBmWWSvagVMPiYiUcKC+HSvy\nitxu+qA1rJmdkg8g/6zv5TqsIiKiIdi6vxa//tshRAf7YfuK2Zg2KlzpkpyKKzaJSJUG3n3w8rHR\nWLdgBiLdbPqgNRjiRKQ69Wd6sDJPgwMn23HflWPw0DXqu/ugvTDEiUhV9lQ2Y9XbxTCaJHIXp+Pa\nyXFKl6QohjgRqcI3mxc/9/FRjI8NwfrF6Rgdrd67D9oLQ5yIXF6HzoCHt5Xik8ONuHFaPJ65eSqC\n/BhfAEOciFzc0YZOrMgrQl1bD35zwyQsVcHmxc7EECcil/Ve6Smszi/DiAAfbF52KWaOjhz8L3kY\nhjgRuRyDyYw1fy/Hm3urkZEUgZdvT0NsaIDSZbkkhjgRuZSmTh1+trkYX1e3YWlWMn7544nw87Hl\nNk+egSFORC6jqOYM7t1UhPZeA/44fxpumpGgdEkujyFORIqTUiLvq1o89f4hxIcH4q07Z2LiRaFK\nl6UKDHEiUpTOYMKv3jmIHZp6XJUaiz/eOh1hQb5Kl6UaDHEiUkxdW//uO4dOdeDB7HFYddU4eHno\n8nlbMcSJSBH/qmjGqi3FMJkl3liSgasnjlS6JFViiBORU0kp8fI//7N8/tXF6Ujm8nmbMcSJyGk6\ndQY8sr0Muw414IZp8XiWy+eHjf/3iMgpjjV1ImdjEWpae/D4Tybi7stHc/m8HTDEicjhdh1swMPb\nShDg6428u2dh9pgopUtyGwxxInIYk1niDx8fxcv/PI5po8KxflEaLgoLVLost8IQJyKH0Pb0YdWW\n/t3nb8schSfnTIa/j/vvPu9sDHEisrvDpzqwPK8Qje16rJk7FQtmJipdkttiiBORXf2t5CRW7yhD\neKAftiy/FGmJEUqX5NYY4kRkF0aTGWs+KscbX1RjZnIkXro9DTEh/kqX5fYY4kQ0bK1dety3WYMv\nq/pvH/urn0yErzdvH+sMDHEiGpayei1WbCxCa3cf/nDLNNycztvHOhNDnIhstr2wDr969yBiRvgj\nf0UWpiaEKV2Sx2GIE9GQ9RnNePrDw9iwrwZZY6Lw4oIZiBrB8W8lMMSJaEiaOnW4b5MG+0+cwT2X\nj8ZjP0qFD8e/FcMQJyKrFdeewco8DbS9ffjTbdMxZ/rFSpfk8RjiRGSVbfvr8Pi7BxEb6o+dKy/D\npHhun+YKbApxIcQ8AFoAKVLKXPuWRESupM9oxlMfHELel7W4fGw0XlwwAxHBfkqXRRZDDnEhRBqA\nKimlRgiRLYRIk1JqHFAbESmsuVOPezcVYf+JM1h+RQoe+eEEjn+7GFuHU54FcA36e+IFdqyHiFxE\nSV3//G9tbx/WLZiBG6fFK10SncOQ/0m19LqrhBBnALTZvyQiUtq2wjrc+uo++HgL7FiZxQB3YbYM\np4Sjfzx8DYDXhBAaKWXVWe/JAZADAImJvHsZkVoYTGY8/cFh/HVfDS4bG4UXF6QhkuPfLs2W4ZQc\nAGuklFohRBWAeQDWDnyD5WJnLgBkZGTIYVdJRA7X0qXHvZs0+Lq6jfO/VWRYUwyllPmWXjcRqdiB\n+nYs31iI1u4+vDB/On46g/O/1WLIIS6lXCuEeNTSC4/kFEMidXunuB6P7TiA6BH+2LEyC1Mu5v1P\n1MSmnriUcu3g7yIiV2Y0mfHMR+V4/YtqzBodiZdvT+P9T1SIKzaJPJC2pw/3v12MPZUtWDI7CY9f\nP4n3/1YphjiRhylv6EDOhiI0tOuw9uZLcGvmKKVLomFgiBN5kF0HT+OhbaUY4e/D/S/dBEOcyAOY\nzRIv7K7Eut2VmD4qHK8uTsfI0AClyyI7YIgTublOnQH/s7UUBUcacUt6An770ykI8PVWuiyyE4Y4\nkRurbunGsg2FqG7pxm9umISlWckQQihdFtkRQ5zITX1e0Yz7N2vg7SWw8a6ZyBobrXRJ5AAMcSI3\nI6XE63uqseajIxg/MgSv3ZGBUZFBSpdFDsIQJ3IjOoMJv9h5AO8Un8SPp8bhuVumIciPp7k746dL\n5CYa2nVYvrEQpfXteOia8bj/qrEc//YADHEiN6CpPYPlG4vQozfi1cXp+OHkOKVLIidhiBOpXH5R\nPX658wDiwgKQd/csTIgLUbokciKGOJFKGU1mrPmoHG98UY2sMVF4aWEaNzD2QAxxIhVq7zHgZ29r\neAMrYogTqc2xpi4s21CI+jM9eGbuVNw2k1sgejKGOJGKfHa0Cas2F8PPxwubl12KzORIpUsihTHE\niVRASonX9lRhzUflmBgXitw70pEQwQU8xBAncnk6gwm/fOcAdmq4gIe+iz8JRC6sqUOHnI1FKKnT\ncgEPnRNDnMhFHahvx7INhWjvNWD9ojRcN+UipUsiF8QQJ3JBH5Sdws+3lyIquH8H+knxoUqXRC6K\nIU7kQsxmiRcKKrDu02PISIrA+sXpiOYO9HQBDHEiF9HTZ8RDW0ux61ADbklPwNM3TYG/D3fgoQtj\niBO5gJPaXtzz10IcbejA4z+ZiLsvH80LmGQVhjiRwopq2rB8YxH0BjPeWJqJKyfEKl0SqQhDnEhB\nO4rq8YudB3BReAC25GRgbCzvQEhDwxAnUoDJLLH2H+V49fMq3oGQhoUhTuRkXXojHtxSjIIjTbh9\nViKeuHEy70BINmOIEzlRXVsPlm0oRGVTF56aMxl3zE5WuiRSOZtCXAiRBiAFAKSU+XatiMhNFZ7o\nv4DZZzLjrTsz8b1xMUqXRG7A1t/hfmEJ7xRLoBPRBewoqsfC175CSIAP3rn3MgY42c2Qe+JCiHkA\n9gOAlHKt3SsiciNms8TafxzF+s+PI2tMFF6+PQ3hQbyASfZjy3BKJvDtkEo2g5zo3Lr1Rjy4tQSf\nHG7EwlmJeJIXMMkBbP2JapVSaoBve+b/RQiRI4QoFEIUNjc3D6tAIjU6pe3FvPX7sPtII35zwyT8\n7qdTGODkELb8VLUCqLI81sLSMx9ISpkrpcyQUmbExHDsjzxLSZ0Wc17ai/q2Hry5NBN3XsYl9OQ4\ntoR4PiwzUwCEwzI+TkTA+6WnMP/VfQjw9cLOe7PwAy6hJwcb8pi4lLJKCKG1DKNEcUycqH8PzD/t\nrsQLBZXITI7A+kXpiOItZMkJbJonLqXMtTzkHHHyeDqDCY/ml+G90lO4OS0Bv5/LW8iS83DFJtEw\nNHfqkbOxEMW1Wqy+LhUrvp/C8W9yKoY4kY2ONnTirrf2o7Vbzz0wSTEMcSIbfHa0CfdvLkaQnze2\nL8/C1IQwpUsiD8UQJxqit/ZW46kPDiM1LhRvLM3ARWGBSpdEHowhTmQlo8mMpz44jA37anDNpJF4\nYf50BPvzFCJl8SeQyAqdOgN+trkYn1c0I+eKFKy+LhXeXryAScpjiBNZoaKxE19Xt2HN3KlYMDNR\n6XKIvsUQJ7JCelIk9qy+EtFcwEMuhnfkIbISA5xcEUOciEjFGOJERCrGECciUjGGOBGRijHEiYhU\njCFORKRiDHEiIhVjiBMRqRhDnIhIxRjiREQqxhAnIlIxhjgRkYoxxImIVIwhTkSkYgxxIiIVY4iT\n25BSKl0CkdMxxMkt1LR2Y85Le3HwZLvSpRA5FbdnI9X7sqoVK/KKAADdeqPC1RA5F0OcVG3L17V4\n/N2DSIoKwhtLMpEcHax0SUROxRAnVTKZJdb8/Qhe/6Ia3xsXjT8vTENYoK/SZRE53bDGxIUQj9qr\nECJrdemNyNlQiNe/qMaS2Un4y9JMBjh5LJt74kKIbADXAFhrv3KILuykthd3v7UflU1d+O2cyVg8\nO1npkogUxeEUUo2SOi2WbSiErs+EN5dm4vvjY5QuiUhxNg2nCCHSpJQF9i6G6Hw+LDuN+a/uQ4Cv\nF3bem8UAJ7KwtSceeaEXhRA5AHIAIDEx0cZDEPUv4Hnps2N47uMKpCdFIHdxOqJG+CtdFpHLGHJP\n3JpeuJQyV0qZIaXMiIlhj4lsozea8PD2Ujz3cQXmTI/HpntmMcCJzmJLTzxFCJGC/t54pCXUNXau\nizxcW3cfVmwswtcn2vBg9jg8cPU4CCGULovI5Qw5xKWU+cC3Qybhdq+IPN7x5i7c9dZ+nG7XYd2C\nGbhxWrzSJRG5LJtnp0gpcwHk2rEWIvz7WAtW5BXB19sLby+7FOlJEUqXROTSOMWQXMY3S+hHRwfj\nzaWZGBUZpHRJRC6PIU6KM5klnt1Vjtx/VeGK8TH488IZCA3gCkwiazDESVHdeiMe2FKCgiONuGN2\nEv73+knw8eYdkomsxRAnxZzS9uLuvxbiaEMHnrxxMpZkJStdEpHqMMRJEQOX0P/lzplcgUlkI4Y4\nOd3fSk7ikfwyjAz1x+Z7ZmHcyBClSyJSLYY4OY3ZLPH8JxX482fHMHN0JNYvSkdksJ/SZRGpGkOc\nnKJbb8TD20qx61ADbsschafmTIGfDy9gEg0XQ5wcrq6tB8s2FKKisRO/vn4S7rosmUvoieyEIU4O\n9e/jLbhvkwYms8Rbd87EFbyASWRXDHFyCCklNn5ZgyffP4zR0cF47Y4MjOYmxkR2xxAnu9MbTfj1\nuwexrbAeV6fG4oXbpiOEKzCJHIIhTnbV2KHD8o1FKKnTYtVVY/Fg9nh4eXH8m8hRGOJkN4Un2rBy\nkwbdeiPWL0rDdVMuUrokIrfHEKdhk1IizzL+fXFEIPLunoUJcVzAQ+QMDHEaFp2hf/x7e1E9rpwQ\ngxfmz0BYEMe/iZyFIU42qz/Tg5V5Ghw42c7xbyKFMMTJJnsqm7Hq7WIYTRK5i9Nx7eQ4pUsi8kgM\ncRoSs1nilc+P47mPj2J8bAjWL07n/G8iBTHEyWranj48tK0Un5Y34cZp8Xjm5qkI8uOPEJGSeAaS\nVcrqtbh3kwaNHTo8eeNk3DE7ifc/IXIBDHG6ICkl8r6qxW/fP4zoEX7Ytnw2ZiRyB3oiV8EQp/Pq\n1Bnw2M4D+LDsNH4wIQbP3zqd9/8mcjEMcTqngyfbcd9mDerP9GL1dalYfkUKpw8SuSCGOP0Xs1ni\nzb3VeHZXOaKC/bEl51JkJkcqXRYRnQdDnL7V0qXHz7eX4p9Hm3HNpJFYe/MliODwCZFLY4gTAOCz\n8iY8kl+GDp0Bv50zGYsu5ewTIjVgiHu4nj4jfvfhEWz6qhapcSHIu2cmUuNClS6LiKzEEPdgmtoz\neHhbKU60diPnihQ8dM14BPh6K10WEQ0BQ9wD6QwmPP9JBV7fU4WLwgKx+Z5LMXtMlNJlEZENbApx\nIUSO5eEYKeVqO9ZDDlZU04ZHtpehqqUbt89KxGM/SuXWaUQqNuQQF0JkAyiQUlYJIbYLIbKllAUO\nqI3sqENnwNpd5dj0VS3iwwKx6Z5ZuGxstNJlEdEw2dITT7F85QKosjwmFyWlxEcHG/DEe4fQ0qXH\nnVmj8dC14zHCnyNpRO5gyGeylDJ3wNM0AFvPfo9luCUHABITE20ujobneHMXnnjvEPZUtmByfChe\nX5KBSxLClS6LiOzI5u6YECINgEZKqTn7NUvQ5wJARkaGtL08skWX3ogXd1fizb3VCPDxxv9ePwl3\nzE6Cj7eX0qURkZ0N53fqbF7UdC1GkxlbC+vwx08q0dKlxy3pCXj0ulTEhPgrXRoROYjNs1OklGst\nj3lhU2FSSuw+0oRndpXjWFMXMpMj8PqSDEwfxaETIndn6+yUZ4UQqwFEArjF7lWRVaSU+OJYC/7w\ncQVK6rRIiQ7Gq4vTce2kkVwyT+QhbLmwWQCAuwIoSEqJvcdasW53Jb4+0Yb4sACsmTsV89IT4Mtx\nbyKPwnlmKmIyS/zjUANe+edxHDjZjpGh/nhqzmTMzxwFfx8ulyfyRAxxFejUGbCjqB5/3VeD6pZu\njI4Oxpq5U3HTjIt5rxMiD8cQd2HHmjqR92Ut8ovq0aU3YkZiOF5amIbrpsTBm7vsEBEY4i6np8+I\nvx9owJava1FYcwa+3gLXXxKPJVnJnG1CRN/BEHcBRpMZe4+34m/FJ/GPQw3o7jMhJToYv/xxKuam\nJSB6BOd5E9G5McQVYjCZ8WVVKz462ICPDzWgpasPIQE+uGFaPG6acTFmjo7kNEEiGhRD3Im0PX34\nvKIZn5U34dPyJnTojAjy88aVE2Jxw7R4XJkaw1kmRDQkDHEH0htNKK7V4t/HWvDFsRaU1GlhlkBE\nkC+yJ47EdVPicMX4GM4wISKbMcTtqFtvREmdFl9Xt6Gwpg2aGi16DSZ4CWDqxWH42ZVj8YPUWExL\nCOfsEiKyC4a4jXQGEyobu3DwVDtKarUordeiorETZgkIAUyMC8X8zFGYPSYKl6ZEISyQu+cQkf0x\nxAdhMkvUn+lBRWMXKho7UdnYifKGTlQ2dcFk7r/LbniQL6YlhOPayXFISwxHWlIEQrnlGRE5AUMc\n/WPXp7U61J3pQW1bD2pbe3CitRtVzd2oae1Bn8n87XvjwwIwIS4E2RNHYlJ8KCbHhyIxMogzSYhI\nEW4d4gaTGWd6+tDS2YeWLj1auvRo7NCjsUOHhnYdTnfocErbi+ZO/X/9PT9vL4yKDERKzAhclRqL\nlJhgjI0NwbiRI9jDJiKX4rIh3ttnQkuXHnqjCTqDGTqDCT19JvT0GdGt7/+zQ2dEl96ITp0BHb1G\ntPca0N5rgLanD23dfejQGc/53w7x98HIsADEhQbgqgmxiA8PRHx4AEZFBiEpKggjQwLgxQuPRKQC\nLhvi75edwqP5ZYO+z8dLICTAB2GBvggL9EVooC8SIgIRFeyHiGA/RAb7IXqEP2JC/BE9wh+xIf4I\n5ibBROQmXDbNZiZH4v/mXYIAX2/LlxeC/LwR6OuDYH9vBPp5IzTAF/4+XhyPJiKP5bIhnhwdjOTo\nYKXLICJyadwGhohIxRjiREQqxhAnIlIxhjgRkYoxxImIVIwhTkSkYgxxIiIVY4gTEamYkFI69gBC\nNAOoGfCtaAAtDj2o87lbm9ytPYD7tcnd2gO4X5uG254kKWXMYG9yeIh/54BCFEopM5x6UAdztza5\nW3sA92uTu7UHcL82Oas9HE4hIlIxhjgRkYopEeK5ChzT0dytTe7WHsD92uRu7QHcr01OaY/Tx8SJ\niMh+OJxCEEKkXeC1eUKIbCHEo86saTgGac+zlj9znFcRkeM4NMQHCwC1BYQV7VFdQAghsgFsP89r\naQAgpSwAoL1QOLqKC7XHIkcIcRxAlZNKGjYhRI7l69nzvK6282iw9qjyPLJ8Of0zcliIDxYAagsI\nK+tVXUBY2nO+eucD0FoeVwHIdkpRwzBIewBgmZRyjOV9Ls/yj1KBlDIXQIrl+cDX1XYeXbA9Fqo6\njyxtuMXyGaQ5O+sc2RMfLADUFhDW1KuqgLBCOIC2Ac+jlCrEjlLU1GsFkIL//KxVWZ4PpLbzaLD2\nACo7j6SUBVLK5ZanKVJKzVlvcehn5MgQHywA1BYQ1tSrtoDwOFLKtZZwiDpPL9ClSClzLb1WAEgD\nUHjWW1R1HlnRHkCl55Gl3uXneMmhnxEvbNqR2gLCCloAkZbH4QBaFaxl2CzjsPMsT1tx7l6gS7L8\nCq45Ry9PlS7UHrWeR1LKtQCWCyHCnXlcR4b4YAGgtoC4YL1qDoizDfgh3Ir/tCMFgCp+vT3bgPYU\n4j9tGINz9wJdVbaUcvU5vq+28+gb52yPGs8jIcTAcfAqAGdfkHXoZ+TIED9nAKg4IAZrjyoDwnLC\nZAw4cQBgNwB800uy9Ia0augFWtGeWy2vHVdDe4D+YLP08r75LNR8Hg3WHjWeR9n475CuApz3GTl0\nsY9lilAV+gf7cy3fK5JSpp/vdVdmZXvaLK+vVa5SchcDpky2oT8obpFSFqj1PBpCe1RzHlnC+lbL\n0/RvLnI66zPiik0iIhXjhU0iIhVjiBMRqRhDnIhIxRjiREQqxhAnIlIxhjgRkYoxxImIVOz/AT0Z\nZg+/9BB+AAAAAElFTkSuQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "omega = sp.linspace(0.1,3,90)+1/200\n", "amp = sp.zeros_like(omega)\n", "x = sp.array([[0,-1,1]])\n", "for i, freq in enumerate(omega):\n", " #print(i,freq,x)\n", " #print(sp.average(x))\n", " x = x-sp.average(x)\n", " try:\n", " t, x, e, amps, phases = ms.hb_time(duff_osc2, x, freq, num_harmonics=4, verbose = False, f_tol = 1e-6)#, callback = resid)\n", " amp[i]=amps[0]\n", " except:\n", " amp[i] = sp.nan \n", "plt.plot(omega, amp)" ] }, { "cell_type": "code", "execution_count": 30, "metadata": { "scrolled": true }, "outputs": [ { "data": { "text/plain": [ "[]" ] }, "execution_count": 30, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD7CAYAAAB68m/qAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XlcXOW9P/DPwzIMy8CwDJNAAmRIyGZMQoiJcalaote9\n3katxrYulbysv1b781dbu93e3l77Sn69v75aW71S27jErabmqjV1IY3WkE0gqyExQELCEjIZGPYZ\nGOb5/cEhQQLMEGbmnDnzeb9evDhnngN8j2Q+PjznnOcRUkoQEVH4i1K7ACIiCgwGOhGRTjDQiYh0\ngoFORKQTDHQiIp1goBMR6QQDnYhIJxjoREQ6wUAnItKJmFD+sIyMDJmXlxfKH0lEFPYqKyvPSCkt\nvo4LaaDn5eWhoqIilD+SiCjsCSHq/TmOQy5ERDrBQCci0gkGOhGRTjDQiYh0goFORKQTDHQiIp1g\noFPY4OpaRONjoFPY+NGmg3jizf1ql0GkWQx0CgvOnj68WdWAKCHULoVIsxjoFBY2VjbA7fFi9bJc\ntUsh0iwGOmmelBKv7DqBwhwz5mUlq10OkWYx0EnzdtQ6UHemG/csZ++caDwMdNK8DbvqkZoQixsW\nTFW7FCJNY6CTpp3ucOGDz1pwe9F0GGOj1S6HSNMY6KRpr316Eh6vxF2X5KhdCpHmMdBJszwDXry6\n+wSumJWBGRmJapdDpHkMdNKsrUfsaG538VZFIj8x0EmzNuyshzU5DsVzM9UuhSgs+Ax0IUSx8rF2\njPa1yueSQBdHkeuEowf/PGrH15bmICaa/Q4if4z7ThFCFAO4XUpZBqBQCFE4ymElQohaAHXBKJAi\n08u76xElBC+GEk3AuItEK0FepuzapJRVoxz2oJRyY8Aro4jl9gzgjYoGFM/NxJQUo9rlEIUNv/6W\nFUI8DmDNGM02ZUjm8TG+tkQIUSGEqLDb7RdaJ0WQ9w6eQmt3Hy+GEk2QX4EupVwHYI0Qwjxam9KT\nT1eGaEa2l0opi6SURRaLZfIVk+5t2FmP3PQEXD4zQ+1SiMKKrzH04ePmdQBKRrSXCCFWKbsOALbA\nl0iR5PCpDnx6vA2rl+UgKopT5RJNhK8eejGANGXbDOXC57CeegXOjbHnK/tEF+yVXSdgiInCqiXT\n1S6FKOz4CvRSDI6RlwDAsIufW5T9KgB3KL302jEumhL5pdvtwZtVjbhxwVSkJRrULoco7Pi6y8WJ\nwVAf+fqSYdvntRNdiLf2NqHL7cE9y3mrItGF4BMbpAlSSmzYWY85U0wozElVuxyisMRAJ03Yc9KJ\nQ80dWL08F4LrhhJdEAY6acLLO08g0RCN2xZnq10KUdhioJPqnD19+Nv+JnxlcTaS4sa9rENE42Cg\nk+o2VjbA7fHyyVCiSWKgk6q8XokXd9SjKDcV87KS1S6HKKwx0ElVH31+Gidae/DNFXlql0IU9hjo\npKoXttcj0xSHf7loitqlEIU9Bjqpps7ehY8/t2P1slzEchELoknju4hU8+KOesRGC9y1jPO2EAUC\nA51U0eX24K+VDbhhwVRkmriIBVEgMNBJFZuqGtDp9uAbl+apXQqRbjDQKeSklHhhRz0WZKegMOe8\nNVOI6AIx0Cnkttc6UHO6C99ckcd5W4gCiIFOIffC9uNISzTgpounql0Kka4w0CmkGtp6UFbdgq8t\nnQ5jbLTa5RDpCgOdQmrDzhMAgNXLOW8LUaAx0ClkXP0DeO3TE7h23hRkm+PVLodId3zOVSqEKFY2\nV0opfzBK+yoATgCFUsp1Aa6PdOTtfU1w9vTjGyvYOycKhnF76EqY3y6lLANQKIQoHNFeCABKu3Nk\nO9EQKSVe2H4cs60mXGpLV7scIl0aN9CllGVSyjXKrk1KWTXikDsx2DsHgDoAxSAaRWV9Gz5r6sA3\nVnCJOaJg8WsMXQjxOIA1ozSZAbQO22fXi0b1/PbjMBlj8JVFXGKOKFj8CnRlbHyNEGLCj/UJIUqE\nEBVCiAq73T7hAin8Nbf34u8HT+HOoulI5BJzREHjawx9+Lh5HYCSEYc4AaQp22YAjpHfQ0pZKqUs\nklIWWSyWydZLYejFHfWQUnIRC6Ig89VDL8YXA7sOAIb11F8HYFO2bQDKAl0ghbfevgG8unvwVsXp\naQlql0Oka74CvRSATQhRAgBSyo3K61uU/Srg7N0wzlEumlKE27SnEc6eftx/+Qy1SyHSvXEHNKWU\nTgyG+sjXlwzbPq+dCBi8VfHP5ccwPysZS/NS1S6HSPf4pCgFzSdHz6DmdBfuv2wGb1UkCgEGOgXN\n+vJjyEiKw00LOasiUSgw0Ckoau1d2HrEjq8vz0VcDGdVJAoFBjoFxfPlx2GIjsLdy3LULoUoYjDQ\nKeDae/qxsbIBtyzKgsUUp3Y5RBGDgU4B93rFCfT2D+C+y/LULoUoojDQKaA8A168sL0ey21pmJ+V\nonY5RBGFgU4B9cGhFjQ6e3HfZXyQiCjUGOgUUH/edgzT0+JRPNeqdilEEYeBTgGzv8GJivo23Lti\nBqKj+CARUagx0Clg1pcfR1JcDO4omqZ2KUQRiYFOAdHc3ot39jXh9qJpMBlj1S6HKCIx0Ckgni8/\nDq+UuJ8XQ4lUw0CnSet09eOVXSdww4KpnPOcSEUMdJq013afRKfbg5Irbb4PJqKgYaDTpPQPePHn\n8mNYbkvDxdMmvOQsEQUQA50m5d39zWhud7F3TqQBDHS6YFJKPPvPOszKTMJVBZlql0MU8XwGuhCi\nRPlYO0b72qHjAl0caVt5jQPVzR148AobovggEZHqxg10ZfHnMmXdUJuyP1KJEKIWQF0wCiTtKv2k\nDhZTHG5dnKV2KUQE3z10G4ChEK9T9kd6UEqZL6UsC2hlpGnVzR345+d23LsijysSEWlEzHiNSs98\nSCGA10c5bKjnXiilXBfI4ki7/vhJHRIM0VjNFYmINMOvi6JCiEIAVVLKqpFtUsp1Su88fbQhGWX8\nvUIIUWG32ydfMamuub0Xb+9twh1F02FOMKhdDhEp/L3LpVhK+YORLyphvUrZdWCUIRkpZamUskhK\nWWSxWCZRKmnFeuUx/wcu52P+RFri110uQ0MpQz1wIcTQEyQVAIbGzvOVfdIxZ08fNuysxy0Ls/iY\nP5HG+HOXy1ohRK0Qom1Y0xYAUIZg7lB66bWjDcmQvqwvP46evgE8dNVMtUshohF8XRQtA5A6yutL\nhm2Xjmwnfepye/D89uNYOc+K2VNMapdDRCPwSVHy2yu76tHe24+Hr2bvnEiLGOjkF1f/AP74yTFc\nPjMDi6ZzEi4iLWKgk1/eqGyAvdONb1+dr3YpRDQGBjr51D/gxbMf12JxjhmX2tLVLoeIxsBAJ5/e\n2deEhrZePHzVTAjBSbiItIqBTuPyeiWe/qgWc6aY8OW5nCKXSMsY6DSuzQebUXO6C9++mr1zIq1j\noNOYBrwSvy07ipmZSbhxwVS1yyEiHxjoNKa/7W/C0dNdeLR4FqK5gAWR5jHQaVQDXonfbTmK2VYT\nbriIvXOicMBAp1G9s68JtfZuPFo8i8vLEYUJBjqdxzPgxW+3HMWcKSZcN3+K2uUQkZ8Y6HSet/Y2\n4diZbnxvZQF750RhhIFOX+AZ8OJ3/ziK+VnJuHaeVe1yiGgCGOj0BRsrG1Dv6MGjxQW875wozDDQ\n6azevgH8puxzFOaYUcynQonCDgOdzlq//RhaOtz44fVz2TsnCkMMdAIAtHX34ZmPalE8NxOXzEhT\nuxwiugAMdAIAPP1RDbrdHnz/ujlql0JEF2jcNUUBQAhRomzmSyl/MEr7KgBOAIVSynUBro9CoKGt\nBy9sr8dXC6dxrVCiMDZuD10IUQygTFkI2qbsD28vBM4uJu0c2qfw8psPjwIC+N7KArVLIaJJ8DXk\nYgMwFOJ1yv5wd2Kwdz7UXgwKK9XNHXhzTwPuW5GHLHO82uUQ0SSMO+Si9MyHFAJ4fcQhZgCtw/bP\nW59MGbIpAYCcnJwLq5KCQkqJX7xzCMnGWDx0FdcKJQp3fl0UVYZSqqSUVRP9AVLKUillkZSyyGKx\nTLhACp73P2vBjjoHHru2AOYEg9rlENEk+XuXS/FoF0QxONwydI+bGYAjIFVR0Ln6B/Cfmw+hwJqE\nuy/hX05EeuAz0IUQJUN3rwxdFBVCmJXm13FuXN0GoCwYRVLg/WnbMZxs7cW/3TwfMdG8e5VID/y5\ny2WtEKJWCNE2rGkLAAwNwSjHOS9kSIZCr6XDhT9srcG186y4bGaG2uUQUYD4uihaBiB1lNeXDNsu\nHdlO2rb2vcPwDEj8+Ma5apdCRAHEv7UjzK46B96sasQDV8xAbnqi2uUQUQAx0COI2zOAH206gGmp\n8fjuNbPULoeIAszno/+kH6Uf16HW3o319y1FvCFa7XKIKMDYQ48Qx89046mtNbhxwVRcPZtznRPp\nEQM9Akgp8ZP/OYi46Cj87OZ5apdDREHCQI8Ab+1twraaM/j+v8yGNdmodjmT1ufx4viZbnS5PWqX\nQqQpDHSdO93hwr+9/RkW55ixelmu2uUExOctnbjq1x+hvOaM2qUQaQoDXceklPjhmwfg9gzgv25f\niOgofSwrl5Y4OO9Ma3efypUQaQsDXcfeqGzAPw6fxuPXzYHNkqR2OQHDQCcaHQNdpxraevCLdw5h\n2Yw03LsiT+1yAsoYG40EQzQDnWgEBroOeb0SP/jrfnilxK9vX4gonQy1DJeWaGCgE43AQNeh57bV\nobzGgR/fOBfT0xLULico0hMNcDDQib6Aga4ze086se69I7huvlXX85ynJhrQ2u1WuwwiTWGg60iH\nqx/febUK1mQj1n11IYTQ31DLEKvJiJYOBjrRcAx0nZBS4ok3D6DJ6cLv7lqElIRYtUsKqixzPOyd\nbrg9A2qXQqQZDHSdeGH7cby7vxn/e2UBluSm+f6CMJedGg8AaHa6VK6ESDsY6Dqws86B/3i3GsVz\nM/HQl/LVLickssyDUxg0OXtVroRIO/wKdCFE4Thta5XPJYEqivzX5OzFwy9XITc9Af/vzkW6vEVx\nNNNTB+/eqW/tUbkSIu3wZ5HoYgBvjHNIiRCiFkBdwKoiv7j6B/DQhkq4PV6Ufr0IyUZ9j5sPl2WO\nhzE2CjWnu9QuhUgzfC5wIaUsE0KMF9YPSik3BrAm8oPXK/HYG/uwr6Edz359CWZm6ufRfn9ERwnk\nW5IY6ETDBGIM3SaEKBZCPB6A70V+Wvv+Yby7vxlPXD8H182fonY5qpiZyUAnGm7SgS6lXCelLAOQ\nrgzPUJBt2FmPZz+uwz3Lc1BypU3tclQzKzMJjc5edHNedCIAkwx0IUSJEGKVsusAcF66KMdUCCEq\n7Hb7ZH4cAfjgs1P42VsHcc2cTPz85vm6fnjIl6Fhpjp7t8qVEGnDBQW6EMKsbFYAKFO285X9L5BS\nlkopi6SURRaL5cKqJADAJ0ft+F+v7MGCaWY8dddixERH9l2nMzNNAIAjLZ0qV0KkDf7c5bIKQNGw\nnjgAbAEAKWUVgDuUtlpln4Lg0+OtKHmxEjZLIl64bykS43xez9a9GRmJSDREY3+DU+1SiDTBn7tc\nNgLYOOK1JcO2S4NQFw2z76QT96//FFNTjHjpgWUwJxjULkkToqMELp5mxt6TDHQigE+Kat6nx1ux\n+rldMCfG4uUHl8FiilO7JE1ZON2M6uYOuPo5pwsRA13DymvO4Bt/2o1MUxz+suZSTE2JV7skzVk0\n3Yz+AYlDzR1ql0KkOga6RpUdasF9z3+KnLQEvM4wH9PinMHr85XH21SuhEh9DHQNenHHcZS8VIE5\nU0x4rWQ5h1nGYU02wpaRiO21Z9QuhUh1vFVCQ7xeiSc3V+O5bcdQPNeK3921CAkG/op8WTEzHZuq\nGtE/4EVshN/KSZGN//o1otvtwbdfrsJz247h3hV5ePbrSxjmfrosPwPdfQPYx7tdKMIxMTSg5nQX\nHtpQiVp7F3560zw8cPkMtUsKK5fmpyNKAB9/bkdRnv4X9yAaC3voKvv7gWbc+vttaO3uw0sPLGOY\nXwBzggFFuWkoqz6tdilEqmKgq8TVP4BfvHMID71chYIpJvztu5fjspkZapcVtlbOs6K6uQMnueAF\nRTAGugoONrbj5qe24c/lg+Plr5fwtsTJKp5nBQB8cKhF5UqI1MNADyHPgBd/2FqD254uR4erHy/e\nfwl+fst8GGL4a5isGRmJmJ+VjLf3NqpdCpFqeFE0RKpOtOHHmw6iurkDN108Fb/8ykWckyXAbluc\njV++W41aexfyLZG1ghMRwB560LX39ONHmw7gq89sR1t3H55ZXYjf313IMA+CmxdmIUoAb+1hL50i\nE3voQdLn8eLlXfX43Zaj6HB58MBlM/DoygIkcdrboLEmG3HZzAxs2tuIR4sLEBUVuYt/UGRiugSY\nlBKbD5zCuvcPo97RgxX56fjJjfMwLytZ7dIiwqol0/DIa3vxSc0ZfKmAC6pQZGGgB8iAV+K9g6fw\nh601ONTcgdlWE9bftxRXFVgiepm4ULv+oqn4j6RqvLj9OAOdIg4DfZL6B7x4e28Tnv6oBrX2btgy\nEvHr2xfitsXZiOaf/CFniInC3ZdMx1Nba3DC0YOc9AS1SyIKGQb6BTrd4cKru0/ild31aOlwY84U\nE566azFuWDCVQa6y1ctz8fRHtfhz+TH8/Jb5apdDFDJ+BboQonCs9UKV9USdAAqllOsCWZzWeL0S\nO+sceGX3Cbx38BQ8XokrZmXgydsW4Jo5mRxa0QhrshG3LsrGq7tP4OGrZ3L6YYoYPgNdCFEM4FkA\n+aO0FQKAlLJMCGEbL/jD2ZFTndi0pxFv7W1Ec7sLJmMMvnFpHu5ZngMb73fWpIevzsemPQ14blsd\nnrh+rtrlEIWEP4tElwkh6sZovhPAh8p2HYBiAGEf6F6vxP7GdpQdasGHh1pwpKUT0VECXyqw4Ikb\n5mLlXCviDdFql0njsFmScNPFWXhpRz0evMKGjCT20kn/JjuGbgbQOmw/fZLfTzWOLjd2HWvFJ0ft\nKKs+DXunG9FRAkvzUvHzm+fh5oVZSGcohJVHimdh84Fm/ObDz/Gfty1QuxyioIvIi6JSSjS1u7Dv\npBO76hzYWdeKIy2dAICkuBh8abYFK+dacdVsC5/oDGP5liSsXpaDl3bW494VeZhlNaldElFQTTbQ\nnQCGVhQwA3CMPEAIUQKgBABycnIm+eMmzu0ZQL2jB7Wnu1Dd3IH9je040NAOR3cfACA+NhpFeam4\ndXEWltvSsSA7hcuY6cgjxQV4c08jntxcjfX3XaJ2OURBdUGBLoQwSymdAF4HUKS8bANQNvJYKWUp\ngFIAKCoqkhdY55j6B7ywd7rR3O5Cc3svTrW70OR0od7RjVp7F0609sCr/NQoARRYTfjy3EwsmGbG\nguwUzM9KZoDrWFqiAd+5Ziae3HwYH3x2CtfOn6J2SURB489dLqsAFAkhVkkpNyovbwGwREpZJYQo\nUu6EcQbrDpe/VJzE9poz6OkbQIerH+29HnT09qO9tx9dbs95x8fHRiMnLQHzs1Jwy8Is5GcmId8y\n+MGLmZHnvstm4M2qRvz0rYNYnp+OZGOs2iURBYU/d7lsBLBxxGtLhm2XBqGuL6izd6PqhBPxsdEw\nGWOQbTZi7lQTUuJjkRIfi0yTEVNTjJhqNmJqcjyS42N4TzidFRsdhXWrLsZX/lCOX20+jF/9Ky+Q\nkj6FxUXRH14/Bz+8fo7aZVAYu3iaGQ9cPgN//OQYrptvxVWzM9UuiSjgOHhMEeOxa2djttWEx/6y\nD6c7XGqXQxRwDHSKGMbYaDx192J093nwvb/shdcb8Gv0RKpioFNEKbCa8POb56O8xoH/+8ERtcsh\nCqiwGEMnCqQ7l07H/sZ2PPNRLWZlJuFfC6epXRJRQLCHThFHCIF/v2U+VuSn44d/PYDK+lbfX0QU\nBhjoFJFio6Pw9OpCZJmN+NYLFfhcmfqBKJwx0ClimRMMeOH+SxAbHYV7ntuFeke32iURTQoDnSJa\nbnoiXv7WMvQPeHH3H3ehoa1H7ZKILhgDnSLeLKsJLz2wDB2uftzx3ztQa+9SuySiC8JAJwJwUXYK\nXitZjr4BL+747x042NiudklEE8ZAJ1LMz0rBX9ZcCmNsNL5WuhNbqlvULoloQhjoRMPYLEnY+NCl\nyMtIwLderMCzH9dCSj5RSuGBgU40wtSUeLyxZgVuuGgqfvX3w3jsjX3o7RtQuywinxjoRKOIN0Tj\nqbsW49HiWdi0pxG3/H4bDp/qULssonEx0InGEBUl8GhxAV66fxnaevpx6+/LsWFnPYdgSLMY6EQ+\nXD4rA39/5Aoss6XjJ/9zEF//026cbOX96qQ9DHQiP1hMcXj+3qX45Vcuwt6TTlz7m39iffkxDHAK\nXtIQBjqRn6KiBO5ZnosPvnclltnS8O/vHMJtT5dzci/SDJ+BLoRYJYQoFkI8Pkb7WuVzSaCLI9Ki\nLHM81t+7FL/92iK0dLjw1Wd24Luv7kGTs1ft0ijCjRvoQohCAJBSlgFwDu2PUCKEqAVQF4T6iDRJ\nCIFbF2XjH49dhe9cMxPvf3YK1/zXR/jV5mo4utxql0cRylcP/U4ATmW7DkDxKMc8KKXMV0KfKKIk\nxsXgsWtnY8tjX8L1F01F6Sd1uGLdVqx77zDauvvULo8ijK9ANwMYPkCYPsoxtvGGZIgiwbTUBPzm\nzkX48HtX4stzrXjm41pcsW4rntxczaEYCplJXxSVUq5TeufpQojzevBCiBIhRIUQosJut0/2xxFp\n2sxME566azHef/RKXD0nE3/adgxXrtuKR17bwwm/KOh8BboTQJqybQbgGN6ohPUqZdcBwDbyG0gp\nS6WURVLKIovFMtl6icJCgXUw2D/+/lX45oo8bKk+jZue2oavPrMdGysbOJUABYWvQH8d50LaBqAM\nAIQQZuW1iqHXAOQr+0SkmJaagJ/eNA/bn7gGP7lxLtp6+vB/3tiHS54sw8/eOohDTZxOgAJH+HqM\nWbkdsQ6ATUpZqrxWKaVcMqy9VWlfN973KioqkhUVzHyKXFJK7D7Wild3n8Dmg6fQ5/GiwJqEWxdl\n45aFWZielqB2iaRBSuYW+TwulPNSMNCJznH29OHtfU14e28TKurbAACLc8y4ZWEWrp0/BdnmeJUr\nJK1goBOFkYa2Hryzrxlv7W3E4VOdAIB5U5Oxcp4VK+dZMT8rGUIIlasktTDQicJUrb0LZYdaUFbd\ngsr6NnglMCXZiGvmZuKKmRm4ND8d5gSD2mVSCDHQiXTA0eXG1iN2fHjoFMprHOhyeyAEsCA7BZfN\nzMDlMzOwJDcVxthotUulIGKgE+lM/4AX+xuc2HbUgfKaM6g60QaPV8IQE4WLs1OwJC8VS3PTsCQ3\nFamJ7MHrCQOdSOe63B7sPubAjloHKurbcLCxHf0Dg+/nfEsiinLTsDjHjAXTUlBgNSE2mpOrhisG\nOlGEcfUPYN9JJyrq21BZ34aK463ocHkAAIaYKMydYsJF2SlYkJ3CkA8zDHSiCOf1StS39mB/gxMH\nG9txoLEdnzV2oNN9LuRnWpJQYE1CwRQTZltNKLCakG2OR1QU76jREgY6EZ1nKOQPNLbjYGM7jpzq\nxOctnWhud509JsEQjVlWEwoyk1BgNWFGRiLyMhKRk5YAQwx79GpgoBOR39p7+1FzuhNHTnXh85bO\nsx9nus5NARwlgOzUeOSlJw6GvPJ5RkYislPjOXwTRP4GekwoiiEibUuJj8WS3DQsyU37wutt3X04\n5ujG8TODH8ccPTh+phubqhrPDt0Ag2E/JdmIaakJyE6NR7Y5Htmp8ZimbGeZ43lrZQgw0IloTKmJ\nBqQmGlCYk/qF16WUcHT3DYb8mW6cbO1Bg7MXDW292H2sFac6XOctoJ2RFIfs1HhkpRhhTTYiMzkO\nVtPgtjU5DpnJRiQbY/hE7CQw0IlowoQQyEiKQ0ZSHIry0s5r9wx40dLpRkNrDxqdvWhs60WjEvhH\nT3dhW80ZdLo8532dMTZqMOBNSuAnG5FpikN6UhzSkwxITzQMbica2OMfBQOdiAIuJjpqcNhlnAnG\nevo8ON3hRkuHCy2dbrS0u85td7hwsLEdZdUtcPV7R/36REM00pPikJZoQEaSAWnDwj49yYDUBAPM\nCQaY42NhToiFyRiLaJ3fvcNAJyJVJBhikJcRg7yMxDGPkVKiy+1Ba3cfznT1wdHlRmt3HxzdfXB0\n9cHRPbjf6HThQGM7HF198HjHvtEj2RgDc4IBKUrIp8THnt02xw++npJw7jWTMRZJcTFIiosJi/8Z\nMNCJSLOEEDAZB4M1N33s4B8ipURHrweObjfaevrQ3tsPZ8/gR3tvv7KvvN7bj8a23rPbI8f8R0ow\nRCMpLgYmYwySjLEwKUGfZBx8zaRsJ8XFDr529tjB4yymOMTFBHeYiIFORLohhBjsYSfETujrhv4S\nGPofQIcS8l0uDzrdnsHPrn50uc/td7k9ON3pOneM24Px7gLf8MAyXD4rY5JnOD4GOhFFvOF/CUxL\n9X38aKSU6OkbQKfLgy53v/LZczbwC6xJgS16FAx0IqIAEEIgMS4GiXExAIyq1MBHu4iIdMJnD10I\nsQqAE0DhaItA+2onIqLQGLeHLoQoBAApZRkA59C+v+1ERBQ6voZc7sRg7xsA6gAUT7CdiIhCxFeg\nmwG0DttPn2A7hBAlQogKIUSF3W6/sCqJiMinoF8UlVKWSimLpJRFFosl2D+OiChi+Qp0J4ChmXfM\nABwTbCciohDxFeivA7Ap2zYAZQAghDCP105ERKHnc8UiIUQJBi942qSUpcprlVLKJWO1j/O97ADq\n/awtA8AZP4/VOj2dC6Cv8+G5aJeezmey55IrpfQ5Zh3SJegmQghR4c+SS+FAT+cC6Ot8eC7apafz\nCdW58ElRIiKdYKATEemElgN93PH4MKOncwH0dT48F+3S0/mE5Fw0O4ZOREQTo+UeOmnAePPzCCFW\nCSGKhRCPh7KmyfBxPmuVzyWhq4gocDQR6L6CIZyCw49zCZvQEEIUA3hjjLawm5htvPNRlAghajF4\nG66mKVN0HWRHAAAB7ElEQVRqlAz9exqlPZzeM77OJWzeM8DgvzPlI+S/G9UDXU8zOvpZa9iEhnIe\nY9UZdhOz+TgfAHhQSpmvHKdZyv+YypTnPmzK/vD2cHrPjHsuirB5zyj13678ty8MdZ6pHujQ14yO\n/tQaFqHhB58Ts4UhW5j0am0492+rDuee1h4STu8ZX+cChNF7RkpZJqVco+zapJRVIw4J6u9GC4E+\n6RkdNcSfWsMlNCKOlHKdEhrpY/QUNUGZ8G7orolCABUjDgmb94wf5wKE4XtGqXXNKE1B/d1oIdAj\nSriEhh90NTGbMoa7Stl1YPSeoqYof65XjdILDDvjnUs4vmeU1dvWDJv3KiS0EOh6mtFx3FrDMTRG\n0tvEbMPOpwLnziEfo/cUtaZYSvmDUV4Pp/fMkFHPJdzeM0KI4ePmdQBGXsgN6u9GC4GupxkdfZ1L\nWIWG8kYqGvaGAoAtADDUk1J6TM5w6CX6cT53KG21Wj8fIUTJ0Bq+Q73WMH3P+DqXsHrPYHBMfHhg\n1wGh+91o4sGiQM7oqDY/z6VVaeei2jRhw26/bMVgeNwupSwLx/fMBM4lLN4zSnDfoewuGbpAGqrf\njSYCnYiIJk8LQy5ERBQADHQiIp1goBMR6QQDnYhIJxjoREQ6wUAnItIJBjoRkU78f4ddJvA8kMsg\nAAAAAElFTkSuQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "omegal = sp.arange(3,.03,-1/200)+1/200\n", "ampl = sp.zeros_like(omegal)\n", "x = sp.array([[0,-1,1]])\n", "for i, freq in enumerate(omegal):\n", " # Here we try to obtain solutions, but if they don't work, \n", " # we ignore them by inserting `np.nan` values.\n", " x = x-sp.average(x)\n", " try:\n", " t, x, e, amps, phases = ms.hb_time(duff_osc2, x, freq, num_harmonics=4, f_tol = 1e-6)#, callback = resid)\n", " ampl[i]=amps[0]\n", " except:\n", " ampl[i] = sp.nan\n", "plt.plot(omegal, ampl)" ] }, { "cell_type": "code", "execution_count": 31, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[]" ] }, "execution_count": 31, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXEAAAD7CAYAAACc26SuAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHwJJREFUeJzt3Xl4VOXB/vHvM5nsCQkkCEHZwqIgiISgCLa1NSi1LlVB\na10Ql1j3vtran639WdvaWlr7tlXfKi6t2reK4k6rYqy4s4QdBBQiO0IgCWRPZuZ5/5gJRgSTTGZy\n5kzuz3XlymzJ3Idk7hyeOed5jLUWERFxJ4/TAUREJHwqcRERF1OJi4i4mEpcRMTFVOIiIi6mEhcR\ncTGVuIiIi6nERURcTCUuIuJi3mg/QW5urh00aFC0n0ZEJK4sWbJkj7W2d1uPi3qJDxo0iNLS0mg/\njYhIXDHGbG7P4zScIiLiYipxEREXU4mLiLiYSlxExMVU4iIiLqYSFxFxMZW4iEg0BPxd8jQqcRGR\nSKvcBI+dDqufi/pTRf1kHxGRbmX1c/DKD4OXPdGvWJW4iEgkNNXCqz+BZU/CUePh/Eeg56CoP61K\nXESks3auhDlXwN4N8LVb4ZTbISGxS55aJS4iEi5rYdEsmHcHpPaCy16C/G90aQSVuIhIOOoq4KXr\nYf2/YfgUOOd/ID2ny2OoxEVEOmrTe/Dc1VC3B6bcAyf+AIxxJIpKXESkvfw+eGcmvPN76DkYriqB\nvDGORlKJi4i0x75t8NxVsOVDGPN9OOP3kJzhdCqVuIhIm9a+Ai/dAAEfnPcwHHeB04kOUImLiBxO\ncwPM+xksfgTyjoepj0HOEKdTfYFKXETkUMrXB4/93rUaTroBTr0TvElOp/oSlbiISGvWwrJ/wKu3\nQWIaXDwHhk12OtVhqcRFRFo07IO5/xWc/2Tw14Pj35l9nU71lVTiIiIA25bAnBnBo1BO/f8w6Yfg\nSXA6VZtU4iLSvQUC8OH98OZdkJkHM16FASc6nardVOIi0n3VlMOLP4ANJTDibDj7L5Da0+lUHaIS\nF5HuqWw+PF8M9VXwnT9C4RWOnTrfGSpxEele/D6Y/1t4917IHQaXvgB9jnU6VdjaVeLGmAJr7dJW\n16cCVUCBtXZmtMKJiERU1dbgqfNbF8DYS+Hbv4OkdKdTdUqba2waY4qAZ1tdLwCw1pYAVS3XRURi\n2tq58ODJsGsNnP8onHO/6wsc2lHiobIua3XThQT3wgndXhSFXCIikdHcAP/+Mcy+OLhc2jVvw+ip\nTqeKmHDGxLOBilbXu34WdBGR9tizAeZcDp+tggnXQ9EvYvLU+c7QG5siEp9WPA1zbwFvMlw0G46e\n4nSiqAinxKuAXqHL2cDegx9gjCkGigEGDBgQdjgRkQ5rrAkOn6z4JwycFDx1PutIp1NFTZtj4ocw\nG8gPXc4HSg5+gLV2lrW20Fpb2Lt3787kExFpv89WwaxTYMVT8I3/B9NfiesCh/YdnTIVKAx9puVQ\nw9BRK1WtDz0UEXGEtbDoYXj4VGishukvwzdvd8XcJ53V5nCKtXYOMOeg22ZFLZGISEfUV8HLN8La\nl2FoEXz3QcjoPiMAemNTRNxrW2lw5sH9O6DoLph4E3jCGSV2L5W4iLhPIAALHoCSX0BmP5jxGvQf\n73QqR6jERcRdavcGZx78ZB6MOAvOvs91Mw9GkkpcRNxj0/vBuU/q9sAZf4DxV7ly5sFIUomLSOwL\n+OHdP8L830DPwXDVm5B3nNOpYoJKXERiW/UueP5q+PRtGD0NzvxvSM50OlXMUImLSOza+FawwBtr\n4Oz7Yewl3X745GAqcRGJPa0Xbuh9dPDMyyNGOJ0qJqnERSS27NsefPNyywfBPe9v/x6S0pxOFbNU\n4iISOz6eFzx8sLkBzp0FYy50OlHMU4mLiPP8zfDmL+GDv0CfUTDt78H1L6VNKnERcVbVVphzBWxb\nFFxx/vTfQGKq06lcQyUuIs5Z/yq88IPgceBTH4NR5zudyHVU4iLS9XxN8OZd8OH90Pe44PBJzhCn\nU7mSSlxEulbl5uDwyfZSGH81nPZrSExxOpVrqcRFpOus+xe8eG1wEYdpf4djz3U6keupxEUk+nxN\nwWljFzwAeWOCBd4rv62vknZQiYtIdFVuDi7csH0JnHANnPar4Ar0EhEqcRGJntbDJxc8ASPPcTpR\n3FGJi0jkfWH45HiY9jcNn0SJSlxEIqtqCzw7I3j0yQnFwaNPNHwSNSpxEYmclpN3bACmPQ7Hftfp\nRHFPJS4inedvDp6888F9Ovqki6nERaRz9m0LDp9sWxRc8/K0u3XyThdSiYtI+D55A54vDu6Ja+4T\nR6jERaTj/D54625474/QZzRc8LjmPnGISlxEOmb/TnjuStj8Poy7HKbco6ljHRRWiRtjpgJVQL61\ndlZkI4lIzNr4VnDptOY6rbwTIzwd/QJjTAFQZq0tAcpC10UkngX8MP8eePJcSO8NxfNV4DEi3OGU\n3wGTCe6Jl0Qwj4jEmppyeP4qKJsPYy6C79wLSelOp5KQDpe4tXapMabMGFMJXB2FTCISKzZ/EJz7\nu74SzvoLFFwGxjidSloJZzglm+B4+G+Bh40xXzqi3xhTbIwpNcaUlpeXRyCmiHQpa+H9P8Pfz4TE\nNLiqBMZNV4HHoHCGU4qB31prq4wxZcBUYGbrB4Te7JwFUFhYaDudUkS6Tn0lvHAtfPwqjPwunH0f\npPRwOpUcRqcOMbTWzjHGFEcqjIg4bPtSeHZ68DDCb88MTmClve+YFs6Y+ExjzG2hvfBeOsRQJA5Y\nC4sfgdd/Chl94IrX4KhCp1NJO4S1J26tndn2o0TEFRpr4JWbYfUcGDoZzpsFab2cTiXtpDM2Rbqz\n3evgmctg7yfwrZ/DybeAp8PHO4iDVOIi3dXKZ+GVm4LHfF/2Egz+utOJJAwqcZHuxtcIr90OpY/C\ngInB2Qd75DmdSsKkEhfpTio3B4dPdi6HiTfBqXdCgmrAzfTTE+kuPn49OPe3tXDh/8KIM51OJBGg\nEheJdwE/vPUbePcP0Hc0XPCElk6LIypxkXhWUx6c+/vTt2HspXDG7zX3d5xRiYvEqy0L4dnLob4C\nznkAxl7idCKJApW4SLyxFhY+CPPugKyj4Mo3IO84p1NJlKjEReJJYzW8fCOseQGOPgO++1dIzXY6\nlUSRSlwkXuxeB89cCns3QNEvYOLNOvuyG1CJi8SDVXPg5ZsgKU1nX3YzKnERN/M1wRs/D46B9z8R\npv0devRzOpV0IZW4iFvt3xE8+mTrQphwHUz+JSQkOp1KuphKXMSNPn0nuPZlU11w7pNR5zudSByi\nEhdxk5a1L9+8C3KGwvS5cMQxTqcSB6nERdyiYT+8dB2sfQVGnhM8gSc50+lU4jCVuIgb7F4Lsy+B\nik/htF/DSTdo7UsBVOIisW/1c/DSjcHDB6e/DINOdjqRxBCVuEis8jfDG3fCggd0+KAclkpcJBZV\n74I5M2Dz+3BCMZx2N3iTnE4lMUglLhJrtiyEZ6dDfRWcOwvGXOh0IolhKnGRWGEtLHoYXr8dsvrD\nVSXQd5TTqSTGqcRFYkFTHcz9IaycDcOnwLkPafZBaReVuIjTKj6F2ZfCrtVwyk/h6z/W7IPSbipx\nESd9UhJcPg0L338Ghp/mdCJxGZW4iBMCAXj3XnjrbuhzLFz4pBYvlrCEVeLGmAIgH8BaOyeiiUTi\nXcM+eOFaWP8vGD0NzvpL8EQekTCEO/B2e6i880OFLiLtUb4eHv4WfPwaTLkHzntYBS6d0uE9cWPM\nVGAxgLV2ZsQTicSrj16GF6+FxFSY/goMmuR0IokD4eyJjwdyjDEFxpjbIh1IJO4E/FByV3D9y97H\nQPHbKnCJmHCHU/Zaa5fCgT3zLzDGFBtjSo0xpeXl5Z0KKOJqdRXwv9PgvT/CuMthxr8h60inU0kc\nCafE9wJloctVBPfMv8BaO8taW2itLezdu3dn8om412erYNYpsOldOOvPwQ9vstOpJM6EU+JzCB2Z\nAmQTGh8XkVZWzYFHJgdnIpzxanAvXCQKOlzi1toyoCo0jJKjQwxFWvH7YN4dwRN4+h0PxfPhqEKn\nU0kcC+s4cWvtrNBFFbhIi7qK4Orzn76t6WOly+iMTZFI2LkSZl8cnAf8nAdg7CVOJ5JuQiUu0lmr\n5sBLN0BqT7jiVThynNOJpBtRiYuEy++Dkjvhw/thwES44HHIOMLpVNLNqMRFwlFXEVw+rWw+jL8a\nTv+Nxr/FESpxkY7aswH+cR5U74Sz74OCy5xOJN2YSlyko9JzIXsAnP8o9P/SuW4iXUolLtJRqdlw\n+VynU4gA4c+dIiIiMUAlLiLiYipxEREXU4mLiLiYSlxExMVU4iIiLqYSFxFxMZW4iIiLqcRFRFxM\nJS4i4mIqcRERF1OJi4i4mEpcRMTFVOIiIi6mEhcRcTGVuMSnnSvh/T87nUIk6lTiEn9WzIZHT4OF\nD0F9ldNpRKJKJS7xw98M/74NXiiGI8dB8fzgKjwicaxTJW6MuS1SQUQ6pWY3PH42LHoIJlwPl70I\nGUc4nUok6sIucWNMETA5gllEwrNtCTz0DdixDM57BKb8BhISnU4l0iW0ULK427J/wNz/gsy+cOU8\nyDvO6UQiXSqsPXFjTIG1tiTSYUTazd8M//oRvHQ9DDgJit9WgUu3FO6eeK+IphDpiJrd8Mx02PIB\nnHQDFN0FCfpPpXRPHf7Nb89euDGmGCgGGDBgQJjRRA5h+1KYfQnUVQTHv4+b5nQiEUeFs/uSb4zJ\nJ7g33itU6ktbP8BaOwuYBVBYWGg7H1MEWP4UvHJz8KiTK1+HvDFOJxJxXIdL3Fo7Bw7sbesgXIk+\nfzPM+zks/CsM+hpMexzSc5xOJRITwh5IbL23LRItFbu3E3hmBrl7FsKE62DyrzT+LdKKXg0Ssz5Z\n/i6ZL84g21ZRfcb9ZJ54qdORRGKOTruXmLTgxb/S/4VzAcu2c19QgYschvbEJaY0NTWx+OEbmFQ+\nm7XJo+l71Wz6HnGk07FEYpZKXGLG7s+2suvR7zOpeSWlfaZx/JUP4E1KdjqWSExTiUtMWFM6n5y5\nVzLM7mN54T0UnnWt05FEXEElLo6y1vLBnD9RuPpuKj3Z7Jr6EsePmuR0LBHXUImLYxrq61j2UDGT\nql7ho9QCjrr6Kfrm9HU6loirqMTFETs2f8L+Jy7mJP96Fh91OeMuvxePV7+OIh2lV410uVXvvMBR\n/7mRHvhYNek+xk++zOlIIq6lEpcu4/f7Wfj4T5mw+SE2Jwwg8aJ/MHqYpo8V6QyVuHSJveU72fro\npUxsWMyS7MmMuPpR0jKynI4l4noqcYm6NQtLyHn1GkbaKpaMvoOC827FeHSysEgkqMQlagL+AB88\ndTcnfPLf7PXksP28Fxl33NecjiUSV1TiEhVVe8vZ+OjlnFz3HiszJjL46ifIy+7tdCyRuKMSl4hb\nt/RtMl+5muMCeyg95lbGXXiHhk9EokQlLhFjAwEWPv1bCtbfS6WnJ5vPfpbCcac6HUskrqnEJSKq\n9uxi02OXM6HuA1akT2DwlU/QJ6eP07FE4p5KXDpt9YJ55Lx2LSNtJR8O/xETLvqZhk9EuohKXMLm\na25mwZN3MGHzLHZ7erP1uy9y0vFfdzqWSLeiEpew7Ni6kb1PzODk5hUsyy5i+BUP0y+rl9OxRLod\nlbh0WOmrTzB04e1k22aWFtxNwdnXgzFOxxLpllTi0m611VWseex6TqicywbvUNK+9zcKhmruExEn\nqcSlXT5a9CaZr95AYWAnHx45ncLLZ5KYlOJ0LJFuTyUuX6mhoZ4lj/+ECTueYLcnl3Wn/5OTJp7h\ndCwRCVGJy2F9smIBnpd+wKTAp5TmfIcRl99PXg+9eSkSS1Ti8iXNTQ2UPnkHhVseo9pksPrrD1L4\nrYucjiUih6ASly9Yt/QdkubeyEmBTZRmTWb49AcYpTMvRWKWSlwAqKmuYtWTt3HCrmeoMNksm/Q/\nFE6+2OlYItKGsErcGFMcujjEWvuTCOYRBywreZq+7/2Mk9jD4txzOObSPzI2O9fpWCLSDh0ucWNM\nEVBirS0zxjxrjCmy1pZEIZtE2fayteyecwtj6z5gs6c/6789h/HjJzsdS0Q6IJw98fzQxyygLHRZ\nXKSmpprlT99F4da/0xMPHw65icLv/UzHfYu4UIdL3Fo7q9XVAmD2wY8JDbcUAwwYMCDscBJZAb+f\nJXNn0X/ZHziZPSzN+hb9L7yXk47U32ERtwr7jU1jTAGw1Fq79OD7QkU/C6CwsNCGH08iZeV7/yJ1\n/p2M933CBu9QPj7tAQpOmOJ0LBHppM4cnVKkNzVj37ol82ma90uOa1zCLnJYNPYeCs8sxpOQ4HQ0\nEYmAsI9OsdbODF3WG5sxaMOK96l+/W7G1r1PFZksGnYLY867lT6pGU5HC4vPH+A/63aTl5XK6KOy\nnI4jEjM6vPxK6OiU3xljNhpjKqOQScJkrWX1wjdYfs9khr5wBkNql7FgwDUk3rKSEy6+k2SXFjiA\nMYZbn13BPxdtdjqKSEwJ543NEqBnFLJImPw+H6v+8xRJix9kVPNqqshkweDrGXnOrUzIznE6XkQk\neAwnDu7FgrIKp6OIxBSdseliVXt3s+bVhxi08UmOt7vYSW9Kj/4Ro86+mQnpPZyOF3ET8nMoWbub\nTXtqGZSb7nQckZig1WxdxgYCrFvyFov+dBHJfzmWSRv+QK23J0tP/BO5P/uIwot+TkocFjjAGaPz\nMAbmLNnmdBSRmKE9cZf4bFsZZW/+jX6bX+CYwFZqbTKrcqbQ+5s/YPjoSU7H6xL9slMpGtGHJz7c\nRPE38umRkuh0JBHHqcRjWPmu7ZS98zTpG15mRMMK+hrLusSRlI78OcOLZnBCnIx3d8TNpw7jzI92\n8cg7Zdxy2tFOxxFxnEo8xuzcvJ5N7z9H+qZ5jGxcwYkmwFbTj8X9r2DAN6/gmCGjnI7oqFFHZnHW\nmH48+HYZZ47px/A+mU5HEnGUStxhDfW1rF9cQu2a1+lb/h75gc3kAVs8R7G0/2UcMeH7DBw5nv4e\nvX3R4hdnjeT9DXu48Z/LeP66iaQn69dYui/99nex+vp6Nq54h+p188nYuYAhDWsYYxppsglsSBnF\nh/3Pp/+E8xgwdDSadebQcjKS+fP3jmf6Y4u46all/PWScSR59UdOuieVeBTZQICdWzfw2doPafx0\nAVl7l5Pf/AmjTDMAZZ6BrO5zNqkjTmPo+NMZmaEzEdvra8N688tzRnHHi6u5+ell/Pl7Y1Xk0i2p\nxCOkqamJHRtXUb5hKU3bV5Je+REDGj+mH/vpBzRZL5uShrIqbyopQyYyYOxk8nPzNI9vJ1wyYSCN\nvgC/mvsRex9dyIOXjKNXepLTsUS6lLE2upMMFhYW2tLS0qg+R1exgQCVe3aya9NH7N++Hl/5BpKr\nNpDbsIl+/h0kGT8AzTaBbd4B7O0xgkDf48kaUsigUSeRnJLm8BbEp5eWb+fHc1bSp0cy911UwPH9\ns52OJNJpxpgl1trCNh+nEg+y1lJbW0PVrq3sL99KbfkWmiu3wr7tpNRtJ7txB338u0g3DQe+xm8N\nOxPy2Js6mIbsoST0GUnOkAL6DzsOrxZY6FLLt1Zx3T+WsKu6ketOGcL13xxKSqJmahT36tYl7vf5\nqNlfRW11JbX7K2morqCppgJfbSWBugpM/V5MfSXehr2kNFWQ7quiZ6CKLFP7pe9VQyrlCX3Yn5xH\nQ0Z/bPYAUvsOI3fACPoOPIaExOQu3TY5vP0Nzdz18kc8t3Qb/Xul8rMzRnD6sX0xxjgdTaTDXF/i\ny7ZU8vzS7fgCAXx+iy9gafYH8AcsU8r/xuCGNXgDjXgDjSTaJpIDDaRQT6ptJNU0feX39lvDfpPJ\nfk8Wtd6eNCb3wpeai03vgycrj6TsfmT3HUjukfmkZfYKd9PFIe9+Us6v565l/a5qCgf25PpvDeWU\n4b1V5uIq7S3xmH1jc2tlPXNX7sCb4CHRY0hIMCR6PHgTDIGG/aQE6vB5kmn0plOXkILfm0rAm0bA\nm0ogKQNPciYJqZl407JIzuhFckYv0nr0Ir3nEaT3yKFnQoKmYoxTXxvWm3/dlMPTi7fywFsbmPG3\nxYzM68E138hnyqi+JHs1zCLxI2b3xEUiockX4MXl23nw7Y2UldeSnZbIeWOP4sLx/Tm6r872lNjl\n+uEUkUgKBCzvb9zD04u3Mm/NZzT7LSPyevDtUX2ZMqovw47I0HCLxBSVuMhhVNQ28eKy7fx71U6W\nbKnEWsjvnU7RiD5MGprLCYN6kZqkIRdxlkpcpB1272/g9TWf8dqaz1j0aQXNfktSgodxA3syaWgO\nhYN6cdxRWaQlxezbRxKnVOIiHVTX5GPRpxW8v2EP723Yy9qd+4Hg0nBH98mkYGA2Y/v3ZEz/LAbl\npONN0Gn+Ej0qcZFOqqhtYvnWSpZtqWLZliqWb62iptEHQJLXw/A+GRzdpwcj8jI5um/wo3dGssbW\nJSJU4iIR5g9YNpbXsGrbPtbvqmbtzv2s/6ya3dWNBx7TI8XL4Nx0BuakMyg3ncG5aQzMSWdwTjrZ\naYkqeGk31x8nLhJrEjyG4X0yv7QQRUVtE+s+28+6ndV8uqeWTXtrWbqlkrkrdxBotY/UI8XLkT3T\n6JeVQl52Cv2yU+mXlUpeVvBynx4pmolROkwlLtJJvdKTmDgkl4lDcr9we6PPz9aKejaFin3z3jp2\nVNWzY18DS7ZUUlXX/IXHGwO5GcmhjyR6ZySTm5kc+pzU6r5keqUnkeDRXr2oxEWiJtmbwNAjMhh6\nRMYh769r8rGjqoGd++qD5V7VwK79DeypaaS8upGy8lrKaxpp8gW+9LUeAz3TkshKSyQ7NZHstCSy\nUhPJSk0ku/VtrS5npyaSmeLVG7JxRiUu4pC0JO9XljwEZ9esbvSxp7qRPTVNBwp+T00jFbVNVNU3\ns6+umd3VDXy8q5p99c1UN/i+8nlTEj1kJHuDHyle0pO8ZKZ4SW+5rfV9yV4yk0P3pXhJS0ogNTH4\nkRK6nKg/Co5SiYvEMGMMPVIS6ZGSSH7v9n2Nzx9gf4OPqrrPS76qvomqumDB1zb6qG70UdPq8s59\nDdQ0hq43+Gg8xN7/4Xg95gul/vllT/B6UgIpiZ/f13I9JTGBZK+HJK8n+DnBQ3Kih6SEBJJa395y\n34HbgvdrOCkorBI3xkwFqoACa+3MyEYSkc7wJnjolZ7UqVWOmv2BA4Ve2xQs/OpGHw1NfuqbQx9N\nfhoOXA5Q3xy63uSnrtlPQ5OfPTVNX35ss59IHBSX4DGtit9zoPiDtyWQnOAh0WvwejwkJgQ/exMM\niQkevB4TnFwvodX9X7gcfExiQuhrQl/bMiGf96DbD/X9ExM8ZKUF/wBHU4dL3BhTAGCtLTHG5Btj\nCqy1SyMfTUSckpjgCY6jp0V+uTtrLY2+AA3Nfpp8ARp9AZr8ARqbg5+bfKEPv//AbY2h21o+t9x/\n8G2NX/g+fhqaA/j8Ppr99sC01s0tn1vf5g/gC1j8gcgecv3DomH8sGh4RL/nwcLZE78QeCN0uQwo\nAlTiItIuxpgDwymxJhAIrl3gCwSCJR8q92Z/y7oGLbd//sfA5w/QHAh9PugPw4i8HlHPHE6JZwMV\nra7nRCiLiIijPB5DkseQhHverI1KUmNMsTGm1BhTWl5eHo2nEBERwivxKqBlzbJsYO/BD7DWzrLW\nFlprC3v3budb6iIi0mHhlPhsID90OR8oiVwcERHpiA6XeMuRKMaYIqBKR6aIiDgnrOPErbWzIh1E\nREQ6zj1vwYqIyJeoxEVEXEwlLiLiYlFf2ccYUw5sbufDc4E9UYzT1eJpe7QtsSuetkfb8rmB1to2\nj9GOeol3hDGmtD3LEblFPG2PtiV2xdP2aFs6TsMpIiIuphIXEXGxWCvxeDv+PJ62R9sSu+Jpe7Qt\nHRRTY+IiItIxsbYnLjGiZfGPw9w31RhTZIy5rSszhauNbfld6HNx1yUSiRzHSrytInBTUbRjW1xV\nFKF5cZ49zH0HVnYCqr6qIGPBV21LSLExZiPBBU5iWmiK5+KW36dD3O+a1wy0a3tc87oJ/bsXOfGz\ncaTE2yoCNxVFO7O6pijgwLYcLuuFBKcjhs9XdopZbWwLwNXW2iGhx8Ws0B+jktC8Rfmh663vd81r\nBtrenhBXvG5C2aeF/u0LurrPnNoTb6sI3FQU7cnqiqJop3hb2SnfJXuv+Xz+u1XG59NBt3DTawba\n3h5wyevGWltirb0mdDX/EDO7RvVn41SJt1UEbiqK9mR1S1F0O9bamaGSyDnM3mBMCC200nK0QwFQ\netBD3PSaac/2gMteN6Gc1xzirqj+bPTGZhdwS1G0U5srO7lFaDx2aujqXg69NxhTQv8VXxov8/h/\n1fa47XVjrZ0JXGOMye7K53WqxNsqAjcVxVdmdWNRHEqrX0zXr+zUaltK+Tz/EA69Nxhriqy1PznE\n7W56zbR2yO1x0+vGGNN6HLwMOPiN2Kj+bJwq8UMWgUuLoq1tcV1RhF48ha1eRABvgvtWdmrHtlwQ\num+jC7alOLS31/Lv79bXDNDm9rjpdVPEF0u6DLruZ+PYyT6hw4bKCL4RMCt02xJr7bjD3R+r2rkt\nFaH7ZzqXVNyq1aGSFQQLY5q1tsTFr5n2bk/Mv25CZX1B6Oq4ljc5u+pnozM2RURcTG9sioi4mEpc\nRMTFVOIiIi6mEhcRcTGVuIiIi6nERURcTCUuIuJi/wdHMvt65GO3UwAAAABJRU5ErkJggg==\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plt.plot(omegal,ampl)\n", "plt.plot(omega,amp)\n", "#plt.axis([0,3, 0, 10.5])" ] }, { "cell_type": "code", "execution_count": 32, "metadata": { "collapsed": true }, "outputs": [], "source": [ "from scipy.optimize import newton_krylov" ] }, { "cell_type": "code", "execution_count": 33, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def duff_amp_resid(a):\n", " return (mu**2+(sigma-3/8*alpha/omega_0*a**2)**2)*a**2-(k**2)/4/omega_0**2" ] }, { "cell_type": "code", "execution_count": 34, "metadata": { "collapsed": true }, "outputs": [], "source": [ "mu = 0.05 # damping\n", "k = 1 # excitation amplitude\n", "sigma = -0.9 #detuning\n", "omega_0 = 1 # driving frequency\n", "alpha = 0.1 # cubic coefficient" ] }, { "cell_type": "code", "execution_count": 35, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array(-0.5478691201747141)" ] }, "execution_count": 35, "metadata": {}, "output_type": "execute_result" } ], "source": [ "newton_krylov(duff_amp_resid,-.1)" ] }, { "cell_type": "code", "execution_count": 36, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[]" ] }, "execution_count": 36, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWwAAAD7CAYAAABOi672AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHOxJREFUeJzt3XlcVWXiBvDnZVdEEFkUUOGyuKMirrmVOGVmy6RmTaZZ\notVkY/3SqWmaZnK0rJnWadLcSh03NEdTW9RyT1lUVFyAC4iILLLIvtz7/v7gauYgF/Tee+659/l+\nPn7geq/wfC7nPLy+57znCCkliIjI+jkoHYCIiJqHhU1EpBIsbCIilWBhExGpBAubiEglWNhERCrB\nwiYiUgkWNhGRSrCwiYhUwsmUX8zHx0cGBweb8ksSEdm8xMTEQimlr7HXmbSwg4ODkZCQYMovSURk\n84QQWc15HadEiIhUgoVNRKQSLGwiIpVgYRMRqQQLm4hIJVjYREQqwcImIjIBS9y9i4VNRHSb9HqJ\nXSl5eGzxYWxMvGj272fShTNERPaguk6HzUk5WHpAC21BBQI83eDiaP7xLwubiKiZrpTXYNXPWVh1\nOAtXKmrRK7AtPprcF/f37ghnFjYRkfLSC8qxdH8GNiddRE29HqO7+eHZ4RoM1nhDCGGxHCxsIqJG\nSClxJKMIS/drsetMPlycHPBoVCCeGRaCMD8PRTKxsImIbqDTS3x3+jIW703HiYul8HZ3wezR4Xhq\nSBf4tHFVNBsLm4gIvxxI/GK/FhmFFQhu3xrzH+6FCf2D4ObsqHQ8ACxsIrJzpVV1WP1zFlYczERh\neQ0igzzx2e+icG/PDnB0sNz8dHOwsInILl0urcbygxn4z5ELKK+px4gIX8waocGQ0PYWPZDYEixs\nIrIrafllWLxXiy3Hc6DTSzwQGYCZIzXoGeCpdDSjWNhEZBcSs4rx+d50/JCSBzdnBzwxsDOeHa5B\nJ+/WSkdrNhY2EdksvV7ix3P5+HxvOuIzi+HV2hmzR4dj6pAuaK/wGR+3g4VNRDanXqfH9pO5+OzH\ndJzLK0OgVyu8+UAPPDagE9xd1Vt76k1ORHSTmvqGU/M+35uOrCuViPBvgw8e64MHIgMssnTc3FjY\nRKR6VbU6rD16AUv2aXH5ajUigzyxeEp/jOnuDwcrOzXvTrCwiUi1rlbXYdXhLCw/kIErFbUYGOKN\nRRMiMTzcx2pPzbsTLGwiUp2iilqsOJiBlYcyUVZdj5ERvvj9PWEYEOytdDSzYmETkWrkXa3GF/u0\nWHPkAqrqdLivZwe8cHcYegdZ/znUpsDCJiKrl11Uic/3pmNjwkXopMSDfQLw/KhQhPsrc9U8pbCw\nichqXbhSiX/9mIZNSRfhIAQe7R+E50aGonN79Sx2MSWjhS2EmACgBIBGSrnE/JGIyN5lXanAp3vS\nsPlYDhwdBJ4c3AWzRoaig6eb0tEU1WRhCyGiAGillElCiBghRJSUMslC2YjIzmQWVuDTH9Pw9bEc\nODkITBncBc+NCoV/W/su6muaMyXyLoAxaBhh7zJzHiKyQxmFFfhkTyr+e/wSnBwEpg4JxqyRGvix\nqH+lycI2jKy1QohiADMslImI7IS2oByf7knDluM5cHFywNNDgxE7UgM/DxZ1Y4xNiXihYf56IYAv\nhBBJUkrtTa+JBRALAJ07dzZXTiKyIWn55fh0Tyq2nrgEFycHPDMsBLEjQuHrob4LMlmSsSmRWAAL\npZQlQggtgAkAFt34AsOByCUAEB0dLc2SkohsQnpBOT7e3VDUbk6OmDFcgxkjNIrfK1Etmn1an5Qy\nzjCaJiJqkQtXKvHR7lR8fewi3JwdETtCg9jhGlVe4lRJxuawFwkh5hpG1948rY+IWiK3tAqf7EnD\nhvhsODoIPDMsBLNGhrKob5PREbaUcpGx1xAR3aigrAaf/ZSGNUcuQEqJJwZ1xgt3h/H0vDvElY5E\nZDLFFbVYvE+LLw9lolanx4SoILw4OgxB7exzZaKpsbCJ6I5dra7Dsv0ZWHYgAxW19XioTwBeiolA\niI+70tFsCgubiG5bZW09Vh7KxOK9WpRW1WFsrw6YMyYCEXZ2USZLYWETUYtV1+mw5sgF/PunNBSW\n1+Kebn54eUwEegXax2VOlcLCJqJm0+klNiVdxIc/nMel0mrcFdYei8d0Rf8u7ZSOZhdY2ERklJQS\nP6Tk4b3vziE1vxx9gjzx3sQ+uCvMR+lodoWFTURNis8swjs7zyIxqxgaH3d89rsojO3VwSbvmWjt\nWNhE1Kizl6/ivW/PYffZfPh5uGLBI70xMToIzo4OSkezWyxsIvqV7KJKfPDDeXx9PAdtXJ0w976u\neHpoCFq5OCodze6xsIkIQMOdyD/dk4bVP2cBAogdrsFzo0Lh1dpF6WhkwMImsnNVtTosO6DF53u1\nqKytx8T+nfCHMeHo6NlK6Wh0ExY2kZ3S6yW+PpaD978/h9zSavymhz/m3tcVYX5c9GKtWNhEduhQ\nWiHmbz+DlNyr6BPkiY8m98PAEG+lY5ERLGwiO5KaV4aFO89iz9l8BHq1wkeT+2J8ZAAcHHiKnhqw\nsInsQEFZDT7YdR7rjl6Au6sTXhvbDVOHBsPNmWd+qAkLm8iGXTug+O+f0lFTr8dTQ4Ixe3Q4vN15\n5ocasbCJbNDNBxTv7emPefd1g8a3jdLR6A6wsIlsTHxmEf667TRO5fCAoq1hYRPZiJySKizccQbf\nJOeio6cbDyjaIBY2kcpV1tbj871aLN6bDiGAl0aHY9bIUC4lt0EsbCKVklJi64lLeGfnWeSWVmN8\nnwD8cWw3BHpxhaKtYmETqVDyxRL8bVsKErKK0SuwLeep7QQLm0hF8suq8d635xCXdBHt3V2w6NFI\nPNo/CI6cp7YLLGwiFajT6bHiYAY+3p2GmnodYodr8Pt7wuDh5qx0NLIgFjaRlTuUVog3t55GWn45\nRnfzwxsP9ECIj7vSsUgBLGwiK3W5tBrzt6fgm+RcdPZujWVTozG6u7/SsUhBLGwiK1Nb3zD98dHu\nVOj0EnNiIjBzpIbX/SAWNpE1OZhWiDf/ewrpBRWI6e6Pv4zvgU7erZWORVaChU1kBXJLqzB/+xls\nN0x/LJ8WjXu6cfqDfo2FTaSgOp0eyw5k4GPD9MfLYyIQO4LTH9Q4FjaRQhIyi/D61ydxPq8cY3r4\n480HOP1BTWNhE1lYaWUd3vn2LNYevYBAr1ZY+lQ0Ynpw+oOMY2ETWci1a3+8/U0KiivrEDtCg5dG\nh8PdlbshNQ+3FCILyCyswJ//ewr7UwvRp5MXvpzeCz0DPJWORSrDwiYyo9p6PZbsS8fHe9Lg6uiA\ntx/qiScGdeG1P+i2sLCJzORoRsNBxbT8cozr3RFvju8B/7ZuSsciFWNhE5lYaVUdFu44g3Xx2Qj0\naoUV0wbg7m5+SsciG8DCJjKh705fxp+3nMKVilrMHKHBSzHhaO3C3YxMw+iWJISIAqABACllnNkT\nEalQQVkN3tp6GttP5qJHx7ZYPm0AegXyoCKZVnN+9b8mpZwohJgrhIiSUiaZPRWRSkgpsTkpB3/7\nJgVVtTq8em9XxI7QwNnRQeloZIOaLGwhxAQA8QAgpVxkkUREKpFTUoXXN5/E3vMF6N+lHd59NBJh\nfm2UjkU2zNgIewBwfVokprHSFkLEAogFgM6dO5s8IJG10eslVh/Jwrs7z0IC+OuDPTFlcBc48FQ9\nMrPmTIlckVImCSFihBATbp7HllIuAbAEAKKjo6U5QhJZi8zCCrwadwLxmcUYHu6DBY/05vU/yGKM\nFfYVAFrD5yVoGHHzwCPZHb1eYtXPWXhn51k4Owq8NyESE/oHQQiOqslyjBV2HIAJhs+9YJjPJrIn\n2UWVmBuXjMPaKxjV1Rfv/DYSHTy5AIYsr8nCllJqhRAlhoOP7XngkeyJlBJrj2bj79tTIITAu4/2\nxqToThxVk2KMzmEb5qgBToWQHcktrcK8TSex73wBhoa2x6IJkQhqx7lqUhaXYBHdQEqJTUk5+Ou2\n06jXSfztoZ54chDPACHrwMImMigoq8Frm09i15k8DAhuh/cm9EGwj7vSsYiuY2ETAdh9Jg9z45JR\nVlOPN8Z1x9N3hfASqGR1WNhk16pqdZi/PQVrjlxA945tsXZyX0T4eygdi6hRLGyyW8kXS/CH9ceR\nUViB2BEavPKbCLg68W7lZL1Y2GR3dHqJz/em44MfzsPXwxVrnhmEoWE+SsciMoqFTXYlu6gSL284\njvjMYoyL7IgFD/eGZ2tnpWMRNQsLm+yClBJbjufgzS2nIQH8c1IfPNIvkItgSFVY2GTzymvq8cbX\nJ7Hl+CVEd2mHDx7ryws2kSqxsMmmncopxYtrjyHrSgXmxETghbtD4cSbC5BKsbDJJkkp8eWhTCzY\ncRbt3J2xdsZgDNK0VzoW0R1hYZPNKamsxdy4ZHyfkod7uvnh/Yl94O3uonQsojvGwiabkphVhNlr\njyO/rBpvjOuOZ4aF8MAi2QwWNtkEvV7i833p+Mf35xHo1Qpxs4aiTycvpWMRmRQLm1SvsLwGc9Yf\nx/7UQoyL7IiFv+2Ntm48t5psDwubVC0hswjPr0lCaVUdFjzSG48P5A0GyHaxsEmVpJRYfjATC3ec\nQVC7Vvhy+kB079hW6VhEZsXCJtUpr6nHvLhkbD+Zi9/08Mf7k/pwCoTsAgubVCU1rwwzVycis7AC\nr43thtgRGk6BkN1gYZNq/Pd4Dl7bfBKtXZyw5tnBGBLKhTBkX1jYZPVq6/VYsOMMVh7KxIDgdvj0\niSj4t3VTOhaRxbGwyaoVlNXghTVJOJpZhGeHhWDe2G5w5rVAyE6xsMlqnbxYithVCSiurMXHj/fD\ng30ClI5EpCgWNlmlLcdyMG9TMnzauCJu1lD0CvRUOhKR4ljYZFXqdXq8s/Mslh7IwKAQb3z2uyi0\nb+OqdCwiq8DCJqtRUlmLF9cew/7UQkwd0gVvPNCD89VEN2Bhk1U4d7kMM75KwOXSaix6NBKTBnRS\nOhKR1WFhk+L2nM3Di/85BndXJ6ybORhRndspHYnIKrGwSTFSSqw4mIn521PQI6Atlj41AB08eX41\n0a2wsEkR9To93tp2Gqt/voB7e/rjg8f6orULN0eipnAPIYu7Wl2HF9YkYX9qIWaNDMXce7vCwYHX\nAyEyhoVNFpVdVInpK+ORUVjBg4tELcTCJotJzCpC7FeJqNdLrHpmEC/eRNRCLGyyiG0nLuGVjScQ\n4OmG5dMGQOPbRulIRKrDwiazW7pfi/nbz2BgsDcWT+mPdu4uSkciUiUWNpmNXi/x9x1nsOxABu7v\n3QH/nNQXbs6OSsciUq1mr/sVQsw1ZxCyLTX1OsxedwzLDmRg2tBgfPJ4FMua6A41q7CFEDEAxpg5\nC9mIq9V1mLr8KL5JzsVrY7vhL+N7wJGn7RHdMU6JkEldLq3GtBVHkV5Qjg8f64uH+wUqHYnIZhgd\nYQshoqSUuywRhtQtNa8Mv/3sIC4WV2HFtIEsayITa84I29vsKUj1ErOK8fSKo3B1dsT6mYPRM4A3\nHCAytSZH2M0ZXQshYoUQCUKIhIKCAtOmI1XYn1qAJ5cegbe7CzY/N5RlTWQmxkbYGiGEBg2jbG9D\ngSfd+AIp5RIASwAgOjpamicmWatvT+Vi9trjCPVrg6+mD4SvB+8OQ2QuTY6wpZRxUso4w0MvC+Qh\nFdmYkI3n1yShd5An1sUOZlkTmVmzTuuTUi6RUobePLom+7X8QAZejUvGXWE+WPXMQHi2clY6EpHN\n42l91CJSSny4KxUf7U7F2F4d8OHkvnB14oIYIktgYVOzSSmxcOdZLNmnxaToICx4pDeceJNcIoth\nYVOzSCnx120pWHkoE1OHdMFbD/aEEFy9SGRJLGwySq+XeHPrKaz++QKeHRaCP43rzrImUgALm5qk\n10u8tvkk1idkY9bIUMy7ryvLmkghLGy6JZ1e4tW4E9iclIPZ94RhzpgIljWRgljY1Kh6nR4vbziB\nrScu4eUxEZg9OlzpSER2j4VN/6Nep8ecDSew7cQlzL2vK54fFaZ0JCICC5tuotNLzI1LxrYTlzDv\nvm54blSo0pGIyIAn0dJ1er3EHzclY/OxHLwyJoJlTWRlWNgEoOE86z9tOYWNiRcxe3Q4XuScNZHV\nYWETpJT4y9bTWHv0Ap4fFYo5MSxrImvEwrZzUkq8/c0ZfHU4CzOGh+DVe3meNZG1YmHbufe/P4fl\nBxvubP76/VzBSGTNWNh2bPHedPzrx3Q8PrAz/jK+B8uayMqxsO3U2qMXsHDnWTwQ2RHzH+7FsiZS\nARa2Hfom+RJe//okRnX1xT8n9YWjA8uaSA1Y2Hbmp3P5mLP+OAZ08ca/f9cfLk7cBIjUgnurHYnP\nLMKs1YmI8PfA0mnRaOXCO8UQqQkL206cvlSK6SvjEeDZCl9OH4i2brwHI5HasLDtgLagHE8tOwoP\nVyesenYQfNrw7uZEasTCtnG5pVWYsuwoAGD1s4MQ6NVK4UREdLt4tT4bVlpVh2nL41FaVYd1sYOh\n8W2jdCQiugMcYduomnodZq5KgLawHIun9EevQE+lIxHRHeII2wbp9RIvbziBn7VF+GhyX9wV5qN0\nJCIyAY6wbdDfd5zB9uRcvH5/NzzUN1DpOERkIixsG7N0vxbLDmTg6buCMWO4Ruk4RGRCLGwbsvXE\nJczffgbjenfEn8fxYk5EtoaFbSMOpRXilQ3HMTDEG/+Y1AcOvD4Ikc1hYduAM7lXMXNVIkJ83PHF\nlGi4OXPJOZEtYmGrXE5JFaatOAp3VyesfHogPFtzyTmRreJpfSpWUlmLqcuPorJWh42zhiCAqxiJ\nbBpH2CpVXafDs18m4MKVSiyZEo1uHdoqHYmIzIwjbBXS6SX+sO44Ei8U45PH+2FIaHulIxGRBXCE\nrTJSSvxt22l8e/oy3hjXAw9EBigdiYgshIWtMov3afHl4SzMGB6CZ4aFKB2HiCyIha0iW47l4J2d\nZzG+TwBeG9td6ThEZGEsbJU4mFaIV+NOYLDGG+9PjOTCGCI7xMJWgZRLDQtjND5tsHhKNFyduDCG\nyB4ZPUtECBFr+DRUSjnPzHnoJheLKzFtxVF4uDlh5fQB8GzFhTFE9qrJEbYQIgbALinlEgAaw2Oy\nkJLKWkxbEY+qOh2+nD4QHT25MIbInhmbEtEAuFbSWsNjsoDqOh1mfNWwMOaLp6IR4e+hdCQiUliT\nUyKGkfU1UQDWmzcOAUC9To/Za48hPrMYnz7RD4M1XBhDRM086CiEiAKQJKVMauS5WCFEghAioaCg\nwOQB7Y1eL/HHzSfxfUoe3hrPhTFE9IvmniUSc6sDjlLKJVLKaClltK+vrwmj2R8pJf6+4wziEi9i\nTkwEpt3FhTFE9AujhS2EiJVSLjJ8zoOOZvSvH9Ou395r9ugwpeMQkZVpzlki7woh0oUQxRbKZJdW\nHc7E+9+fx2+jAnl7LyJqlLGDjrsAtLNQFru1IT4bb249jZju/lj0KFcxElHjuNJRYRviszFvczKG\nh/vi0yf6wcmRPxIiahzbQUEbEn4p6yVT+vNejETUJBa2QjYkZGPepmQMC/NhWRNRs/COMwq4Ng0y\nLMwHXzzFu5wTUfOwsC1s6X4t5m8/g+HhLGsiahkWtoVIKfHed+fw2U/puL93B3zwWF9eJpWIWoSF\nbQF1Oj3+vOUU1sVn4/GBnTH/4V5w5Kl7RNRCLGwzK62qwwtrknAgrRC/vzsMr/wmgotiiOi2sLDN\nKLuoEk+vjEfWlQosmhCJSdGdlI5ERCrGwjaTvecL8NK6Y5AS+Gr6IAwJ5SVSiejOsLBNTK+X+GRP\nGj7cfR5d/T3w7yf7I8THXelYRGQDWNgmVFBWg//beAJ7zxfgkX6BWPBIb7Ry4ZkgRGQaLGwT+SEl\nD3/clIzymnq8/XAvPDmoMw8uEpFJsbDvUGlVHRZsP4P1Cdno0bEt1k3ui3Def5GIzICFfZuklNhx\n8jLe2nYaV8prMGtkKF4eEwEXJ16ehYjMg4V9G7QF5Xj7mxT8eK4AvQLbYvnUAegd5Kl0LCKycSzs\nFiiprMVHu1Ox6nAW3Jwd8ca47pg2NJjXsCYii2BhN0NZdR1WHMzE0v1alNfUY/LAzpgTEwFfD1el\noxGRHWFhN6G0sg6rj2Thi/1alFTWIaa7P/7v3gh069BW6WhEZIdY2I3QFpRjxcFMxCVeRFWdDvd0\n88MfYsIRGeSldDQismMsbAOdXmJ/agG+OpyFPWfz4eLogAf7BmD6XSHoEcARNREpz+4LOy2/DHGJ\nOfj62EXkXa1Be3cXvDQ6HE8O7sI5aiKyKnZZ2NlFlfju9GVsS87FiewSODoIjIrwxVvjg3BPdz/e\nWICIrJJdFLaUEufyyvD96Tx8e+oyUnKvAgB6dGyLP93fHQ/1C4Cfh5vCKYmImmazhX2lvAYH0gqx\n73wh9qcWIL+sBgDQv0s7/On+7ri3Zwd0bt9a4ZRERM1nE4UtpcTF4iokZhUjPrMICZnFOJdXBgDw\nau2MYWE+GBHui1FdfeHXliNpIlInVRb21eo6pFy6ilM5pTieXYKEzGJcvloNAPBwdUJUl3YY36cj\nhof7olegJ++fSEQ2waoLW6+XuFRahdT8cpzNLcOpS6U4nVOKzCuV118T4OmGgSHeGBDcDv27eKNr\nBw8WNBHZJKso7HqdHheKKpGWX47U/HKkGz6m5Zejqk53/XWdvFuhV4AnJkZ3Qs+AtugZ4MlT74jI\nblhFYW9OysHcTcnXH3f0dEOYXxtMHtgJYX5tEO7ngQj/NvBq7aJgSiIiZVlFYQ8JbY/3JkQi3N8D\nob7u8HBzVjoSEZHVsYrC7uTdGp28eYodEVFTeCFnIiKVYGETEakEC5uISCVY2EREKsHCJiJSCRY2\nEZFKsLCJiFSChU1EpBJCSmm6LyZEAYCs2/znPgAKTRbGdJir5aw1G3O1DHO1zJ3k6iKl9DX2IpMW\n9p0QQiRIKaOVznEz5mo5a83GXC3DXC1jiVycEiEiUgkWNhGRSlhTYS9ROsAtMFfLWWs25moZ5moZ\ns+eymjlsIiJqmtWMsIUQUU08N0EIESOEmGvJTKRe1ro9Gcn1ruFjrOUSkZpYRWELIWIAbLzFc1EA\nIKXcBaCkqQ3eDLma3LEttYM1I4ciBWQt708j39dat6db5jKIFUKkA9BaKBKAhp+P4c+7t3heqe3L\nWC7Fti/DH4u/X1ZR2Iad51Yb6WMASgyfawHEWCJTM3dss+9gxnIoVUDW8v40xhq3J8BoLgCYIaUM\nNbzOIgy/RHZJKZcA0Bge3/i8UttXk7kMLL59GXJMNLwfUZbeH62isI3wAlB0w+P2Fvq+zdmxLbGD\nGcuhVAFZy/vTUkptT82hUWAkq8EvPzut4fGNlNq+jOUCFNi+pJS7pJQzDQ81Usqkm15i1vdLDYWt\nlObs2JbYwYzlUKqArOX9sRlSykWG8ml/ixGlOb7nEsMoFgCiACTc9BJFtq9m5AIU3L4M33NmI0+Z\n9f2yyD0dbzHHpG3mb8YSAN6Gz70AXDFZsDskpVwEAEKIMUKIGCsbSSrOSt8fq9yeDPtIkZQyDg2Z\nGhtRmvP7RwFIamTEqKimcim5fUkpFwkhNhpWN5YY/xemYZHCvuE3ZbMJIbwMb8R6ANeWe2oAmOyH\nYuQXSZM7tgV3MGMFo1QBWcv70yyW2J5uxw25EvDLXGwogMUWjhIjpZzXyN8r/Quu0VxKbV83zFEn\noeHnFQtg0Q0vMev7ZRVTIkKICQCiDR+v2Q1cf2OuTfaXmHIEcO2/XTf9ubYDr8cvG8H1HVsI4WX4\nuwT8srOHovH/spmCsRyNPm8B1vL+/A+lticT5JpkeC7dwrlibxitxhg+Kr19Gcul1PYVg18Xsvam\nXGZ9v7hwpgmG3+JaNBxcWGL4u0QpZf8bni8yPL/o1l/JIjl+9bwlWMv7Q7fvhlMNi9BQRBOllLuU\n3r5akMui25ehmCcZHva/dgDSUu8XC5uISCWsYkqEiIiMY2ETEakEC5uISCVY2EREKsHCJiJSCRY2\nEZFKsLCJiFTi/wGUUWNQIREFwwAAAABJRU5ErkJggg==\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "sigmas = sp.linspace(-1,3,200)\n", "amplitudes = sp.zeros_like(sigmas)\n", "x = newton_krylov(duff_amp_resid,1)\n", "for i, sigma in enumerate(sigmas):\n", " try:\n", " amplitudes[i] = newton_krylov(duff_amp_resid,x)\n", " x = amplitudes[i]\n", " except:\n", " amplitudes[i] = newton_krylov(duff_amp_resid,0)\n", " x = amplitudes[i]\n", "\n", "plt.plot(sigmas,amplitudes)" ] }, { "cell_type": "code", "execution_count": 37, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/Users/jslater/anaconda/lib/python3.6/site-packages/scipy/optimize/nonlin.py:474: RuntimeWarning: invalid value encountered in double_scalars\n", " and dx_norm/self.x_rtol <= x_norm))\n" ] }, { "data": { "text/plain": [ "[]" ] }, "execution_count": 37, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD7CAYAAAB68m/qAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHxdJREFUeJzt3XlwHOd95vHvO7jvwTEACAo8QIqkSIqiANA6Ql02bMmW\nY0sJJXkdHxuvBTtaO5vdxJRVZe8mXq+9orLlxE580LcSx5HEkg/ZjhSBsm4pIkiREknxBEmBIAni\n4OC+8e4f0wBBCMCAxACN7nk+Vajpmbcx8+sazDMv3n6721hrERER7wu4XYCIiMSGAl1ExCcU6CIi\nPqFAFxHxCQW6iIhPKNBFRHxCgS4i4hMKdBERn1Cgi4j4ROJcvlhBQYFdsmTJXL6kiIjn7dy5s9la\nG4q23rQC3RhTbq3dNVkbUAZgrd021fMsWbKE2tra6bykiIg4jDEnprNe1CEXY0wV8NgUqzzgBHmZ\nE+4iIuKCqD10a22NMaZuojZjzCZgh7PelhjXJiIiF2GmO0U3APnGmHJjzOZYFCQiIpcmFrNcWkbG\n150eu4iIuGCmgd4CjAzHhIn02C9gjKk2xtQaY2qbmppm+HIiIjKZSwp0Y0zQWdyGM8MFCOKMp49l\nrd1qra201laGQlFn3YiIyCWaziyXTUDluOGU7QDW2jog7LTlR5u2KCIis2c6s1y2EemJj32sYszy\nVmdx1sL8H545zI7j5/jpp941Wy8hIuJ5njj0/1z3ADuOt6Lrn4qITM4TgV6YlUJ3/xCdfYNulyIi\nMm95ItCLslMBONvR53IlIiLzlycCvTArBYCz7Qp0EZHJeCPQs51A7+h1uRIRkfnLE4EeynKGXNRD\nFxGZlCcCPTs1kdSkgHroIiJT8ESgG2MozEqlUT10EZFJeSLQIbJjVD10EZHJeSbQi7JTNYYuIjIF\nzwR6KCtF89BFRKbgmUAvzE6hs2+QLh0tKiIyIc8EelGWjhYVEZmKZwJ99OCidu0YFRGZiHcCXT10\nEZEpeSbQi5weeqN66CIiE/JMoOekJZGcGKBJPXQRkQl5JtCNMYQyNXVRRGQyngl0iAy7aMhFRGRi\n0wp0Y0z5NNbZPPNyplaYlaoeuojIJKIGujGmCnhsGuu8N1ZFTaYwO0XTFkVEJhE10K21NUDdHNQS\nVVF2Ku29g/QODLldiojIvDPjMXRjTLkT+rMupEvRiYhMKhY7RfNi8BzTMnJt0UadRldE5B1mFOjT\n6Z0bY6qNMbXGmNqmpqaZvByL8tIBONHSPaPnERHxo0sKdGNM0FksM8ZsMsZUA3kTzYax1m611lZa\naytDodBMaqU0L52EgOFYc+eMnkdExI+mM8tlE1Dp3I7YDmCt3Wat3eY8FnzHL8dYUkKARXnpHGvu\nmu2XEhHxnMRoKziBvW3cYxXj7m8Ftsa2tIktLcigrkmBLiIynqeOFIVIoB9v6WJ42LpdiojIvOK5\nQC8LZdA7MMwZHWAkInIBzwX60oIMAI2ji4iM47lALyvIBKCuSTNdRETG8lygF2WnkJ6cQJ166CIi\nF/BcoBtjWFqQoSEXEZFxPBfogAJdRGQCngz0soIM6lu76R8cdrsUEZF5w5OBvjSUwbCFt1t1ThcR\nkRGeDPSRmS4adhEROc+Tgb7EmYuuqYsiIud5MtBz0pIoyEzmqAJdRGSUJwMd4IoF2bx1usPtMkRE\n5g3PBvrqkmwOnunQTBcREYdnA31tSQ79Q8McPqteuogIeDjQ15RkA7DvVLvLlYiIzA+eDfQl+Rlk\nJCewr6HN7VJEROYFzwZ6IGBYXZLNXvXQRUQADwc6wJqSHN463c6Qrl4kIjK9QDfGlE/RVu38PBi7\nsqZnTUk23f1DOmJURIRpBLoxpgp4bIq2Guci0WXO/TmzdmEOAPtOaRxdRCRqoFtra4C6SZrLgJEQ\nr3Puz5nlhZkkJwY000VEBEicyS87PfMR5cAjMyvn4iQlBFhVnMVezXQREYnNTlFnjH2XtXZXLJ7v\nYqwpyWbfqXas1Y5REYlvsZrlUmWtvX+iBmeHaa0xprapqSlGL3fe+tIgbT0DOlGXiMS9Swp0Y0xw\nzHK1tXaLs/yOnaLW2q3W2kprbWUoFLr0SiexYUkeADuOn4v5c4uIeMl0ZrlsAiqd2xHbnbYq4EFj\nzFFjjCuJurQgg4LMZHYca3Xj5UVE5o2oO0WttduAbeMeq3Bua4Dc2SlteowxVC7OY8cJBbqIxDdP\nHyk6onJJLvWtPZxp63W7FBER1/gi0N+1dGQcXb10EYlfvgj01QuySU9OoFaBLiJxzBeBnpgQ4OpF\nQV7TTBcRiWO+CHSITF88cKad9t4Bt0sREXGFrwLdWtipXrqIxCnfBHr5olySEwO8eKTZ7VJERFzh\nm0BPS07gmqV5PH8o9qcXEBHxAt8EOsCNl4c4fLaTU+Eet0sREZlz/gr0FZFzxaiXLiLxyFeBvqIo\nk+LsVJ4/rEAXkfjjq0A3xnDD5QW8eLiZwaFht8sREZlTvgp0iAy7tPcOsuekrmIkIvHFd4G+cXkB\nAaNxdBGJP74L9NyMZK4qDfLMgbNulyIiMqd8F+gA71tdzJsNbTRo+qKIxBFfBvqta4oA+Pd9Z1yu\nRERk7vgy0MtCmawsyuLJvQp0EYkfvgx0gFvXFrPjeCvNnX1ulyIiMiemFejGmPIp2jYZY6qMMZtj\nV9bM3bqmiGELNfsb3S5FRGRORA10Y0wV8NgkbeUwerHo8FTBP9dWL8imNC+NJzWOLiJxImqgO2Fd\nN0nzPUDYWa4DqmJU14wZY7htTTEvHWmmrUcXvRAR/5vpGHoQGHshz/wZPl9MfXBdCQNDln9787Tb\npYiIzDrf7hQFWHdZDmWhDB5/vcHtUkREZt1MAz0M5DnLQaBl/ArGmGpjTK0xprapaW4PxzfGcOf6\nhbx2rJWT57rn9LVFRObaJQW6MSboLD4ClDnLZUDN+HWttVuttZXW2spQKHRpVc7AHVcvBOBXu0/N\n+WuLiMyl6cxy2QRUOrcjtgNYa3c561QB4ZH780lpXjrvWpLH47tOYq11uxwRkVkznVku26y1udba\nbWMeqxizvNVaW2Ot3TpbRc7UHVcv5GhTF3sb2t0uRURk1vh6p+iI269cQEpigEdr690uRURk1sRF\noOekJ3H7lQv4xesNdPUNul2OiMisiItAB/joNYvo7BvkiT3aOSoi/hQ3gV6xOJeVRVn8y2tvu12K\niMisiJtAN8bw0WsW8cbJNt7U9UZFxIfiJtAhMtslNSnAz/7jhNuliIjEXFwFek5aEnesX8gvdzfQ\n2tXvdjkiIjEVV4EO8KmNS+kdGOZnr6qXLiL+EneBvqIoi5tWhPjpKyfoGxxyuxwRkZiJu0AH+PQN\nS2nu7NP5XUTEV+Iy0DcuL2BVcRY/fOGYzu8iIr4Rl4FujOHTN5RxsLGDZw6cdbscEZGYiMtAB/jw\n+hIuy03jm88cUS9dRHwhbgM9KSHAfTcvZ099mBcON7tdjojIjMVtoAP8ccVCSnJS+fvth9VLFxHP\ni+tAT0lM4LM3L2PniXO8fPQdV88TEfGUuA50gLsrSynOTuWhpw6qly4inhb3gZ6alMBfVF3O7vow\nT+1rdLscEZFLFveBDrCp4jKWhTJ46KkDDA4Nu12OiMglUaADiQkBvnDrSo42dbFt50m3yxERuSQK\ndMeta4pZXxrkGzWHdJk6EfGkqIFujNlkjKkyxmyO0l4d+/LmjjGGL39wNY3tfXz72SNulyMictGm\nDHRjTDmAtbYGCI/cH9de57TXjW/3morFudx59UK+/8Ix3m7pdrscEZGLEq2Hfg8QdpbrgKoJ1nnQ\nuS2z1u6KVWFu+eL7V5EYMHz1t/vdLkVE5KJEC/Qg0Drmfv7YRifA64wx58atN8oYU22MqTXG1DY1\nNc2o2LlQlJ3Kf71lOf++v5FnD+rEXSLiHTPaKWqMCRLpwX8d+L4xpmz8OtbardbaSmttZSgUmsnL\nzZlP37CUZaEMvvyrvfT06yIYIuIN0QI9DOQ5y0Fg/PHx1cDXrbVbgHuBTbEtzx0piQl87c4rqW/t\n4e+2H3K7HBGRaYkW6I8AI73uMqAGRnvmF7DWbuP8eLvnXVOWzz2VpfzghWPsO9XmdjkiIlFNGegj\nOzmNMVVAeMxOz+1O+xag2pm6WG2t3Tqr1c6xBz6witz0JDZve4MBHUEqIvNc1DF0Zwy8ZmxYW2sr\nxixvsdZu81uYAwTTk/nqHWvZd6qdf/y95qaLyPymI0WjuG3tAu68eiH/8MwR9jZo6EVE5i8F+jT8\n9R+uIT8zmf/x6G56BzTrRUTmJwX6NOSkJ/HQpqs41NjJ//ntW26XIyIyIQX6NN24IkT1jWX806sn\neGrfGbfLERF5BwX6Rfir963kyoU5bN72Bg3hHrfLERG5gAL9IiQnBvjmf7qaoWHLfT/bRd+gxtNF\nZP5QoF+kpQUZPLRpHXvqw3z1NxpPF5H5Q4F+Cd5/5YLR8fTHd+kKRyIyPyjQL9HmW1dyzdI8vvj4\nm+yu980ZD0TEwxTolygxIcB3PlZBYVYK9z5cy+k27SQVEXcp0GcgLyOZH35yA919g9z7cK2uRSoi\nrlKgz9DK4iy+9dGreet0B5/95530D+okXiLiDgV6DLx7VRFfv/NKXjjczBe27WF42LpdkojEoUS3\nC/CLuzeU0tTZx0NPHSQ/I4Uvf/AKjDFulyUicUSBHkP33byM5s4+fvTSMQqzU/jsTcvcLklE4ogC\nPYaMMXz59tW0dPbzf//tAHnpydy9odTtskQkTijQYywQMPztXVdxrruf+x9/AwzcXalQF5HZp52i\nsyA5McD3P1HJxuUFbN72Bj9/7W23SxKROKBAnyWpSQl8/xOV3LwyxAOPv8k/v3rC7ZJExOeiBrpz\nAegqY8zmSdrLnXU2xb48b0tNSuB7H6/gPasK+dIv9/LTl4+7XZKI+NiUgW6MKQew1tYA4ZH74zxg\nrd0GlE3SHtdSEhP4zscqeO/qIv7Xr/fx7WePYK3mqYtI7EXrod8DjJx5qg6oGtvo9Mp3AFhrt1hr\nd8W8Qh9ITgzw7T8p50NXlbDlyYP8z1/tY0gHH4lIjEWb5RIEWsfczx/XvgFGe/JV1totMazNV5IS\nAvzdPetZEEzle8/Vcaa9l29+5GrSkhPcLk1EfCIWO0VbRnrmE42jG2OqjTG1xpjapqamGLycdwUC\nhgfefwV/86E11LzVyEd/8CqtXf1ulyUiPhEt0MNAnrMcBFrGtbcQGYoZWXfD+Cew1m611lZaaytD\nodBMavWNT16/hO/8SQX7T7XzR99+iSNnO90uSUR8IFqgPwKUOctlQA2AMSboPLZtTHsQZzxdortt\nbTH/cu81dPQOcsc/vsTT+xvdLklEPG7KQB8zlFIFhMfs9NzutNcRmf2yCch3ZrvINFUszuOJz2+k\nLJTBvQ/X8o2nD+lMjSJyycxcTqGrrKy0tbW1c/Z6XtE7MMSXfrmXbTtP8p5VhXzjI+vJTk1yuywR\nmSeMMTuttZXR1tORovNAalICD21ax1c+vIbnDjXxoW+9yJsn29wuS0Q8RoE+Txhj+MR1S/h59bX0\nDQ7zR995iR+8UKchGBGZNgX6PLNhSR6/+/MbuGVlIV/97Vv86U920NTR53ZZIuIBCvR5KDcjme99\nvIL/fcdaXq1r4f1//wK/P3DW7bJEZJ5ToM9Txhg+fu1ifv25jeRlJPGnP9nBFx7bQ1vPgNulicg8\npUCf51YWZ/Hrz23kvpuX8fjrDbzvG8+pty4iE1Kge0BqUgKbb1vFL+67npy0SG/9Lx/dQ7hbpw0Q\nkfMU6B6y7rIgT3x+I5+7ZTm/3N3Au//fczxaW6+ZMCICKNA9JyUxgb+6dSW//fONlBVksHnbG9z1\nvVd463S726WJiMsU6B61qjibRz9zHQ9tWsex5i4++K0X+coT+2nr1k5TkXilQPewQMBwV2Upz/zl\nTdyzoZQfv3yMm/729/zkpWMMDA27XZ6IzDEFug8E05P52p1X8pvPb2T1gmz++on93PqN53l6f6Mu\ndycSRxToPrKmJIefffoafvjJSoyBex+u5SNbX6X2eGv0XxYRz1Og+4wxhvdcUcSTf3EjX/nwGo42\ndbHpu6/wn3/8mk74JeJzOn2uz/X0D/HTV47z3eeOEu4e4LY1xfy3qsu5YkG226WJyDRN9/S5CvQ4\n0d47wI9ePMYPXjhGZ98g715VyJ/dvIwNS/Ki/7KIuEqBLhNq6x7g4VeO8+OXj9Pa1c+GJbncd/Ny\nbl4ZwhjjdnkiMgEFukypu3+QR3bU8/3n6zjV1ssVC7L5s5uX8YG1xSQmaNeKyHwSs0B3rhcaBsqt\ntVumWG/zVO2gQJ+P+geH+dXuBr773FGONnWxMJjGx65dzEc2lJKbkex2eSJCjC5BZ4wpB7DW1hC5\nGHT5JOtVAe+9lELFXcmJAe6qLOXp/34T3/1YBaV5aTz45AGu/fp2Nm/bw94GzYwR8YrEKO33AE87\ny3VAFbBrVisSVwQChtvWFnPb2mIOnGnn4VdO8ItdDTxae5LKxbl88vol3La2mCQNx4jMW9E+nUFg\n7FEp+eNXMMaUOz148YlVxdl87c4refWB9/Cl26/gbEcfn//561z39e187XdvceRsh9slisgEovXQ\np0Pz3nwqJz2JT99Qxqf+YCnPHWriX3e8zY9ePMbW5+u4elGQeypLuX3dArJSk9wuVUSIHuhhzgd2\nEGgZ2zid3rkxphqoBli0aNElliluCgQMt6wq5JZVhTR19PHL1xt4pLaeLz7+Jn/zxH4+cOUC/rhi\nIdcszSchoKmPIm6ZcpaLsxO00lq71RizGaix1u4yxgSttWFnBgxEQv8zwL3W2knH2DXLxT+steyu\nD/NobT1P7DlNZ98gRdkpfHBdCR9eX8KVC3M0r10kRmI5bbGayA7RMmvt1jFPXjFunfuBuxTo8aen\nf4iatxr59Z5TPHvwLANDliX56XzoqhI+tL6E5YVZbpco4mk6sEhc0dY9wJP7TvPrPad45WgLwxZW\nFWdx29pibl1TzKriLPXcRS6SAl1cd7a9l9+8cZon955hx4lWrIXF+encuiYS7leXBglozF0kKgW6\nzCtNHX08vb+Rp/ad4eWjzQwMWQqzUqhaXcQtKwu5flk+GSmxmHQl4j8KdJm32noGePbgWZ7ce4bn\nDzXR1T9EckKAa8ryuGlFiFtWFVJWkKGhGRGHAl08oX9wmNrjrfz+4Fl+f7CJI2c7AViUl87NK0Pc\nsrKQa8vySUtOcLlSEfco0MWT6lu7efZQE88eOMtLR5vpHRgmJTFA5ZJcrl9WwHXL8lm3MEdnhJS4\nokAXz+sdGOK1Y608e7CJl482c+BM5JQDmSmJvGtpHtcvy+f6ZQWsKs7SzlXxtekGuvZCybyVmpTA\njStC3LgiBEBLZx+v1rXy8tFmXjnawjMHzgKQm57EdcvyubYsn4rFuawqztYRqxKXFOjiGfmZKdy+\nbgG3r1sAwOm2Hl452sJLR1p45Wgzv3vzDABZKYlcvTiXysW5VC7JZX1pkPRk/amL/2nIRXzBWsvJ\ncz3sPHGOHcdb2XniHAcbO7AWEgKGtSXZVCzOGw34BTmpmkUjnqExdIl7bT0D7Hr7HLXHW6k9fo7d\n9WH6BocBCGWlcNVlQdaX5nBVaZB1C4PkpOuskTI/aQxd4l5OWhK3rCzklpWFQGSK5P7T7eypD7On\nPszuk2Fq3mocXb+sIIOrSoNcdVkOV16Ww6ribB3sJJ6iv1aJG8mJAdaXBllfGhx9rK1ngDdPtrHn\nZJjd9WFePNLML15vAMAYWJqfweqSbFaXZLOmJIc1JdkUZKa4tQkiU1KgS1zLSUti4+UFbLy8AIiM\nxZ9p72VvQzv7T7Wz71Qbr78d5jdvnB79ncKsFNaMC/nS3HRNnRTXKdBFxjDGsCAnjQU5abx3ddHo\n423dA+w73cb+U5Gg33+6necPNzM0HNkHlZaUwOVFmawoymKFc7uyOIvibO18lbmjnaIil6h3YIhD\njR3sP9XOwcYODjd2crCxg6aOvtF1slISWVF8PuRXFGVxeVEmocwUBb1Mm3aKisyy1KQE1l0WZN1l\nwQseP9fVz6HGDg6d7eTQmQ4ONnbw5N4z/Py1+tF1slITKQtlUlaQEfkJZbK0IIOlBRk6b41cMgW6\nSIzlZiRzTVk+15Tljz5mraWps4/DjZ0cauygrqmLuuZO/qOuZXQn7IiFwTTKQpFwHwn7slAGJTlp\nGqeXKSnQReaAMYbCrFQKs1L5g+UFF7R19w9yrLkrEvJO0B9r7uLxXQ109g2OrpeSGGBRXnrkJz9y\nu9i5vSw3ndQk9ezjnQJdxGXpyYnObJmcCx631tLU0cfRkZBv6uJEazf1rd28fLSFnoGhC9Yvzk49\nH/RjQn9RXjp5Gckas48DCnSRecoYQ2F2KoXZqVy3LP+CNmstzZ39vN3axdut3bzd0sOJ1i7qW7t5\n/lATZ8fsmAVIT06gJJjGwmAaJcE0LstNoySYysJgOgtz0yjKStEpiX0gaqAbYzYBYaDcWrtlgvZq\nZ3GZtfb+GNcnIhMwxhDKSiGUlULF4rx3tPf0D1F/rpu3W7o50dpNw7keToV7aAj38GZDG61d/Res\nnxAwFGenOiGfxsLctNEvgJEvAR01O/9N+Q4ZY8oBrLU1xpgyY0y5tXbXmPYqoMZaW2eMecwYU2Wt\nrZnlmkUkirTkhNFpkhPp6R+iIXw+5E+Fe2g4F1muPXGO37xxmsHhC6c0Z6UkUpSTSnF2KsXO7ej9\n7FSKclIoyEjRjlsXRfvKvQd42lmuA6qAXWPay5yfrU57WawLFJHYS0tOYHlhJssLMydsHxq2nO3o\nHQ35U+FeGtt7OdPWy5n2Xl460szZjr7RA6tGJAYMhVkpo0FfNCb8C53/KEJZKeSkJWlMfxZEC/Qg\n0Drm/gUDedbarWPulgOPxKguEXFRQuD8EbOTHc0yNGxp7uwbDfmxgd/Y3suhxg5eONx8wUydEckJ\nAQoyk0cDPpSVQigzZdz9VEJZKZqXfxFiMijmDM3sGjscM6atGqgGWLRoUSxeTkTmgYSAocjphV81\nxXodvQM0tvdytqOPppGfzvPLDeFedte30dLVx0QHrmemJF4Q+AWZyeRnppCXkUx+RnLkNjOZvIwU\ngmlJcT3kEy3Qw8DIHpcg0DLJelWT7RB1evFbIXLo/6UUKSLelZWaRFZqEssLJx7PHzE4NExrd//5\n0B8X/E0dfbx1pp2mjj46et/Z6wcIGMhNj4T8+aCPhP1o+GckkzfyeHqyr2b3RAv0R2D0P64yoAbA\nGBO01oad5eqR2S/aKSoilyoxITB68FU0fYNDnOsaoKWrj9auflq7+mnpdG67+ml1Hj9wpoPWrn7C\n3QOTPldOWhL5GcnkZiQTTEsimJ5MMD0psjz6WBK56cnkOMuZKYnzch/AlIFurd1ljKl0ZrOExwyp\nbAcqnMcfNMbcT6Qnf9fslisiAimJCRTnJFCcEz38IdL7P9c94AR+JOzPjYZ/5Dbc3c/ptl4OnOkg\n3N1PV//QpM+XGDAE05PISYsEfWQ5mdz0SODnpEe+CM63JVGUnUpy4uz+NxB1DH3cjs+Rxyqc2xog\ndxbqEhGJmcSEwOjOVph66GdE3+AQbT0DtHUPcK57gHB3P+Ee57Z74ILlU+Fe9p9qJ9wzQPckXwT/\n9F/exQ2Xh2K4Ve+kIwVERCaQkphAYVbCtIaAxhr5Igh3j/xEQn/lJMcExJICXUQkhi71iyAW/LN7\nV0QkzinQRUR8QoEuIuITCnQREZ9QoIuI+IQCXUTEJxToIiI+oUAXEfEJYyc6X+VsvZgxTcCJOXip\nAqB5Dl5nPoiXbY2X7QRtq1/NZFsXW2ujnjdgTgN9rhhjaq21k52X31fiZVvjZTtB2+pXc7GtGnIR\nEfEJBbqIiE/4NdDfccpfH4uXbY2X7QRtq1/N+rb6cgxdRCQe+bWH7kvOxbgna9tkjKkyxmyey5pm\nS5RtfdC5rZ67ikTmP88HerQg88uH37nc32OTtJXD6BWkwlOFoRdMta2OamPMUaBujkqaNcaYaufn\nwUnaffFFPY3t9MXnFCJ/v87PnL+nng70aQaZLz78zjZOtg33AGFnuQ6ompOiZkmUbQW411q7zOsX\nJHe+uGqcyzyWOffHtvviizradjp88Tl1tu0u5z0rH/+ezfZ76ulAZ3pB5osPfxRBoHXM/Xy3Cpkj\nZX7otQJlnP+brXPuj+WXL+po2wk++Zxaa2ustZ9x7pZZa3eNW2VW31OvB/p0gswvH35xWGu3OB/8\n/El6e55grd065iLs5UDtuFV88UU9je0En31One34zARNs/qeej3Qo/LLhz+KMJDnLAeBFhdrmVXO\nOOwm524LE/f2PMX5t3vXBL05X5lqO/32ObXWbgE+Y4wJzuXrej3QpwwyP374xxrzx/II57etDPD0\nv60TGbOttZzfvmVM3Nvzmipr7f0TPO63L+oJt9NPn1NjzNhx8zpg/E7eWX1PvR7oEwaZHz/8zh98\n5Zg/fIDtACM9HqdnE/Z6T28a23q303bUB9ta7fTmRt4/X35RR9lO33xOiYyJjw3sOpi799TzBxY5\n05zqiOyA2Oo8ttNaWzGmvdVp3+JepSIXGjM9s5VICNxlra2Z4O/3gr9vr7mI7fT859QJ7ruduxUj\nO0jn6j31fKCLiEiE14dcRETEoUAXEfEJBbqIiE8o0EVEfEKBLiLiEwp0ERGfUKCLiPjE/wcqKU9C\n1kJvFAAAAABJRU5ErkJggg==\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "sigmas = sp.linspace(-1,3,200)\n", "sigmasr = sigmas[::-1]\n", "amplitudesr = sp.zeros_like(sigmas)\n", "x = newton_krylov(duff_amp_resid,3)\n", "for i, sigma in enumerate(sigmasr):\n", " try:\n", " amplitudesr[i] = newton_krylov(duff_amp_resid,x)\n", " x = amplitudesr[i]\n", " except:\n", " amplitudesr[i] = sp.nan#newton_krylov(duff_amp_resid,0)\n", " x = amplitudesr[i]\n", " \n", "\n", "plt.plot(sigmasr,amplitudesr)" ] }, { "cell_type": "code", "execution_count": 38, "metadata": { "scrolled": true }, "outputs": [ { "data": { "text/plain": [ "[]" ] }, "execution_count": 38, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWwAAAD7CAYAAABOi672AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAH0RJREFUeJzt3Xl4lNWhx/HvyWQHQgwEkCVCUIG6IUZwxwV36opU3KpW\nUW+1t7Yu9ba2V2/VK1pt3XqlVau4iwstVWuhKCJuIYIgWEAWWQQCMQSyz8y5f5wJRIVMAjPzzjvz\n+zzPPDOZmSS/Z5j55eW85z2vsdYiIiLJL8PrACIi0j4qbBERn1Bhi4j4hApbRMQnVNgiIj6hwhYR\n8QkVtoiIT6iwRUR8QoUtIuITmbH8Yd27d7f9+/eP5Y8UEUl5c+bM2WitLY72vJgWdv/+/SkvL4/l\njxQRSXnGmJXteZ6GREREfEKFLSLiEypsERGfUGGLiPiECltExCdU2CIiPqHCFhGJhQScvUuFLSKy\nq0JBWPQ3eOJ0+OTpuP+6mB44IyKSFuq/hoqn4KM/w+YvoWs/CGTH/deqsEVE2mvjUvjwjzD3WWiu\ng72OglPuhH1PhUD861SFLSLSFmth+Tvw/iOw5B9uS/qAsXDY1dDrgIRGUWGLiOxIcwPMfxE++CNs\nWAidiuHYW6Dscujcw5NIKmwRkda2rIeP/wzlj0PdRui5P5z5CBwwBjJzPI2mwhYRAfhqntuanj8Z\nwkEYdCocdg30PxqM8TodoMIWkXQWDsG/33BFvXIWZHVyQx4jroJuA71O9x0qbBFJP0218Mkz8MHD\n8PUK6FoCJ90BB18EeYVep9spFbaIpI+tG+DDR90YdUM19BsBJ94Og05PyLS83ZX8CUVEdlflYnj/\nQZj3AoSaYMhoOPw6KBnhdbIOUWGLSGqyFlbOhtkPwuI3IDPXDXkc/uOkHJ9uDxW2iKSWUBA+/xu8\n9wCsrYD8bm7+9KFXQKfuXqfbLSpsEUkNLTsS338IqldC0UAYfT8cNA6y8rxOFxMqbBHxtx3tSDz5\nTjePOiPgdbqYUmGLiD9t+gJmPwBzn/P1jsSOUGGLiL98NQ9m3Q8Lp0BGFhx8IRx+rW93JHaECltE\nkp+1sPI9V9RLp0F2FzjiJ3DYf0CXnl6nSxgVtogkr3AYFr/pinr1R27FvBN+DWU/SuojEuMlamEb\nY8YA1UCptXZi/COJSNoLNcOCl2HW76FyERSWwGn3unnUKTLjY1e0WdjGmGHAMmtthTFmlDFmmLW2\nIkHZRCTdNNW5cyPOftCdeqt4CJw9EfY/BwJZXqfzXHuGRO4GTsRtYU+Lcx4RSUf11W5a3gd/dGtQ\n9x0Op02AfU6GDJ0rvEWbhR3Zsl5mjPkauDJBmUQkXdRVwQePuHnUjTWw9yg46nrY68ikWYM6mUQb\nEinEjV/fBfzJGFNhrV32reeMB8YDlJSUxCuniKSSrZVuMaaPH4OmrTB4NBxzI/Qe6nWypBZtSGQ8\ncJe1ttoYswwYA0xo/YTIjsiJAGVlZTYuKUUkNdR85Q52KX8Cgg2w39lwzA3Qcz+vk/lCu6f1WWsn\nR7amRUQ6pnoVvPd7qJjkTr914Fg4+ufQfR+vk/lKtDHsCcaYmyJb10Wa1iciHVK13M2hnvssYGHo\nBW6MuqjU62S+FHUL21o7IdpzRES+YeMSePc++PQFtwDTIT+EI38Khf28TuZrOtJRRGKn8t/wzgT4\n7BUI5LiT2R7xEyjY0+tkKUGFLSK7b+NSeOdumP8SZOXDEde5BZk69/A6WUpRYYvIrqtaBu/cA58+\n707BdeRP3Ba1z8/skqxU2CLScV+vhJn3uJ2JgSy3at6R/6kt6jhTYYtI+1Wvgnd/B59MAhOA4Ve6\nWR9denmdLC2osEUkupq1btZHxZNubepDLoWjfgZd+3idLK2osEVk57asc/Ooy58AG4KDL3YHvGh6\nnidU2CLyXbWbYNZ9bgW9ULM74OWYG2GPvbxOltZU2CKyXeMWeP9hmP0QNNfCgT+AkTfpyMQkocIW\nEWhugPLH3A7Fuk0w5Ptw/K1QPMjrZNKKClsknYWCMO9ZePtuqFkNpcfBCbdCn0O8TiY7oMIWSUfh\nMCyaAv+6AzYtcQV91iNQOtLrZNIGFbZIOrEWvpgO02+Hr+ZB8WD4wTMw+HSd4cUHVNgi6WLVRzDt\nNlg5y52F/Kz/c+tSZwS8TibtpMIWSXUbl8C0/4bPp0KnHnDqPe7Al8xsr5NJB6mwRVLV1g1uBb3y\nJyArD477JRz+Y8ju5HUy2UUqbJFU01Tr5lK/9wdoroeyy2DkL6BzsdfJZDepsEVSRSgIc5+BGXfC\n1nVuLvUJv9F5E1OIClvE76yFJW/BP38DlYug73AY+ySUHOZ1MokxFbaIn62pgH/+Gla8C0UDYewk\nt2WtKXopSYUt4kfVX7opegsmQ343OO1eN/MjkOV1MokjFbaInzRudavozX4ITAYcfYM700tugdfJ\nJAFU2CJ+EA7DvOdg+m2wdT0cMBZG/Qa69vU6mSSQClsk2a2cDW/+wh1K3qfMHUre71CvU4kHVNgi\nyerrFW7mx8LXoKAPnPNnOGCMdiimMRW2SLJp3OLOn/j+w26dj2P/C464DrLzvU4mHlNhiySLcNgd\n+DL9dqjdAAeeDyf8Wie6lW1U2CLJYPUceP0GWFvhDnwZ9zz01UkE5JtU2CJeqt3oVtL7ZBJ07gnn\n/AkOOE/j1LJDKmwRL4SCUP44zPitW6zpiOvgmJs0n1rapMIWSbSVs+H1G2H9Aig9Fk6doJPdSruo\nsEUSpeYrt+7H/Behaz8Y+xQMOUPDH9JuKmyReAs2wYd/hHcmQKjZDX0cdb2m6UmHqbBF4mnFLJj6\nM9j4b9j3VDjlTigq9TqV+JQKWyQeajfCW7fCvGehcC+44EXY92SvU4nPqbBFYikchrlPu7Hqxi1w\n1M/gmBs1/CExEbWwjTHDgFIAa+3kuCcS8asNi2Dq9fDl+1ByBIy+D3oM8TqVpJD2bGHfYq09zxhz\nkzFmmLW2Iu6pRPykqQ5mToDZD0JOAZz5MAy9ULM/JObaLGxjzBjgYwBr7YSEJBLxk8Vvwes/d2eA\nGXoRnHg7dOrmdSpJUdG2sA+FbcMio1TaIhE1X8GbN8PCKdB9EFz6OvQ/0utUkuIy2vGcTS3DIJEt\n7m8wxow3xpQbY8orKytjHlAkqYTDUP4EPDwcFv/DraZ39SyVtSREtMLeBCyL3K4mssXdmrV2orW2\nzFpbVlxcHOt8Islj0xfw1Bkw9aew50FwzWw4+ueQme11MkkT0YZEJgMtW9WFRMazRdJKKAgfPAwz\n7oRADnz/ARh2iXYqSsK1WdjW2mXGmOrIUEg3jWFL2lk3H6ZcC1/NhcGj4bR7oWBPr1NJmoo6rc9a\nOzFyU3OwJX00N7ipeu/9AfKK4Lwn4XtnaqtaPKUjHUW+7csP3Fb1piVuPvVJv4X8Iq9TiaiwRbZp\nrod//dad/LawH1z0Cux9gtepRLZRYYsArC6H166BjYuh7EfuAJiczl6nEvkGFbakt2AjvH2XG6su\n6AMXvwYDj/M6lcgOqbAlfa39BF69BioXuWl6J92hcypKUlNhS/oJNsG798LMe6FzD7hwMuxzotep\nRKJSYUt6WbcAXrvaza8+aBycchfk7eF1KpF2UWFLegiH4f2HYPrtrqDPfw4Gn+Z1KpEOUWFL6tu8\nGl69Gla8645W/P4DWgJVfEmFLaltwcvuLDChIJzxEBx8kY5WFN9SYUtqaqiB12+ET5+HvofCORN1\ntnLxPRW2pJ6V78Or42HzGjj2Fjj6BgjorS7+p3expI5QM7z9vzDrPigsgcvfhH7DvU4lEjMqbEkN\nVctg8o9gbYUbpz7lfyGni9epRGJKhS3+t+Bl+Ot/QkYAxj7llkEVSUEqbPGvpjp48xdQ8ST0GwHn\nPuZW2RNJUSps8acNi+Cly6Dyc3dexWNvgUCW16lE4kqFLf5iLVQ8BW/c7MaoL34FBh7vdSqRhFBh\ni3801Lgzli94GUqPhbMnQpeeXqcSSRgVtvjDugXw4sXw9Uo4/lY46meQkeF1KpGEUmFL8pv7LEz9\nGeR2hUv/Dnsd7nUiEU+osCV5NTfAmzfDnL9A/6NhzONu/WqRNKXCluT09Qp48RL4ap4b/jjulzq8\nXNKePgGSfBb/A14Z72aEjHseBp3qdSKRpKDCluQRDsGMO93pu3od6I5aLBrgdSqRpKHCluRQ/7Vb\nC+SL6XDwxXDaPZCV53UqkaSiwhbvbVgEz18A1avg+3+AQy71OpFIUlJhi7c+/7sbr87Kh0unQslh\nXicSSVoqbPFGOAwzJ8Dbd0HvYXD+M1DQ2+tUIklNhS2J17jFnRT386lw0AUw+n7IyvU6lUjSU2FL\nYm36wo1Xb1ziTjIw4mqdFFeknVTYkjjLZ8ILF7uCvvhVKB3pdSIRX1FhS2JUTHIr7XXbGy54Afbo\n73UiEd9RYUt8hcMw/TZ47/du3erz/uIWcRKRDlNhS/w01cGr42HR36Dscjj1Hq0HIrIb2r2gsDHm\npngGkRSzZR385TRYNBVOvgtOv09lLbKb2lXYxphRwIlxziKpYt18+NPxULkYxj0Hh/+HZoKIxIA2\neSS2Fr8Fky9z49SXvwl7Huh1IpGUEXUL2xgzzFo7LRFhxOcqnoLnzoduA+GK6SprkRhrzxZ2UdxT\niL9ZCzPvgRl3wMAT3LKoOZ29TiWSctos7PZsXRtjxgPjAUpKSmIYTXwhHILXb4Dyx+GgcXDGgxDI\n8jqVSEqKtoVdaowpxW1lF0UKvKL1E6y1E4GJAGVlZTY+MSUpNdfDy1e4NUGOuh5O+I12LorEUZuF\nba2dDNu2ogsTkkj8oa4KnhsHqz5086tHjPc6kUjKa9cskdZb0SJUr4Knz4Wvl7sjF/c7y+tEImlB\n0/qkYzYsgklnu6MYL3oFBhztdSKRtKHClvZbM8dtWQdy4PI3oOd+XicSSSvtPjRd0tzyd+HJMyCn\nwB0Qo7IWSTgVtkS3+B/wzBgo6OPKumiA14lE0pIKW9q24GV3hpjiwXDZGzrvooiHVNiyc3OehMk/\ngr7D4Yd/hU7dvE4kktZU2LJjH/wf/O0nsPcJcNHLOumASBJQYct3zX4I3rwZhnwfzn8OsvO9TiQi\nqLDl22b9Ht76JXzvLBjzBGRme51IRCI0D1u2m3kv/Ot/YP8xcPajOkOMSJLRJ1Kct++Gt++EA8+H\nsx6BjIDXiUTkW1TY6c5amHEnzJwAQy90y6OqrEWSkgo73c24w518YNglMPoPkKHdGiLJSp/OdPbO\nPSprER/RJzRdvfcAzPitG7NWWYv4gj6l6ejDR+Gft8J+58CZD6usRXxCn9R0U/4EvHETDB4N50zU\n1D0RH1Fhp5O5z8LU62Gfk91BMTpZroivqLDTxfzJMOXHUDoSxj6lIxhFfEiFnQ4W/hVeGQ8lh7u1\nQbJyvU4kIrtAhZ3qlkyDyZdDn0Pgghe0kJOIj6mwU9mXH8ALF0GPwXDhS5DTxetEIrIbVNipat0C\neHasO0PMRa9CXqHXiURkN6mwU9GmL2DS2ZDVCS55DToXe51IRGJAk3BTTc1XMOksCAfh0qlQWOJ1\nIhGJERV2KqmrclvWdVXuHIzFg7xOJCIxpMJOFY1b3Zh11Rdw4WQ3K0REUooKOxUEG91skDVzYOwk\nd3CMiKQcFbbfhUPwypWwbIZbyGnIaK8TiUicaJaIn1kLU38KC6fASXfAwRd5nUhE4kiF7WfTb4OK\np+DoG+CIa71OIyJxpsL2q/cegFn3wyGXwfG/8jqNiCSACtuPPnk6cgKCs+H034ExXicSkQRQYfvN\noqnw1+ug9Dg4e6LOcC6SRlTYfrL8XbfyXu+D4QdPa01rkTSjwvaLtXPhuXGwR393YExOZ68TiUiC\nqbD9YONSePpct+Lexa9CfpHXiUTEA1EPnDHGjI/cHGitvTnOeeTbNq9xizkBXPwadO3jbR4R8Uyb\nW9jGmFHANGvtRKA08rUkSl0VPH0O1FfDRZOh+95eJxIRD0UbEikFWkp6WeRrSYTGrfDMeVC1HMY9\n53Y0ikhaa3NIJLJl3WIY8EJ84wgAzfXw/AWwtsIt5jTgaK8TiUgSaNdOR2PMMKDCWluxg8fGG2PK\njTHllZWVMQ+YdoJN8OIlsHwmnPmIFnMSkW3aO0tk1M52OFprJ1pry6y1ZcXFOhXVbgk1w+TLYMlb\nMPp+GDrO60QikkSiFrYxZry1dkLktnY6xks4BK9eBZ9PhVMnQNllXicSkSTTnlkidxtjvjDGfJ2g\nTOknHHaHmy94GUbdBiOu8jqRiCShaDsdpwF7JChLegoFYcqP4dPn4dhb4Kifep1IRJKUzjjjpVAz\nvHwFLHwNjvsVjLzR60QiksRU2F4JNsJLl8G//w4n/RaOuM7rRCKS5FTYXmiudyfNXToNTrsXhl/p\ndSIR8QEt/pRo9dXw9BhYOh3OeMh3Zb14/RZ+/uI8tjYGvY4iknZU2Im0eTU8fgqs+hDO/TMMu9jr\nRB22pSHIyxWrmTJ3jddRRNKOCjtR1n8Gfz4RatbARS/DAWO8TrRLhpUUMrhXF5798EustV7HEUkr\nKuxE+Pcb8NhJgIXL3oDSkV4n2mXGGC4cUcJna2v4dPVmr+OIpBUVdjxZC+/+zp0pptvecMV06LW/\n16l225kH9yEvK8BT76/0OopIWlFhx0tDjVsXZPrtsP+5cPmbKXPygYLcLMYNL+GVT1YzX1vZIgmj\nwo6Hr+bBxJGwcAqM+m+3gzErz+tUMfXTE/ehW6ccfvXafEJhjWWLJIIKO5ashY/+5HYuNjfApX+H\no64HY7xOFnMFuVncOnoI81Zv5pEZS72OI5IWdOBMrGxeDVOuhWUzYO8T4exHoVM3r1PF1RkH9WbG\n5xu4b9pi9u/bleMG9fA6kkhK0xb27gqHoeIpeORwWPURnH4fXPhSypc1uBkjd51zIIN7FXDtMxXM\nXVXtdSSRlKbC3h3rFsATp7qlUXsdANe8B4f+KCWHQHYmLzvAE5ceSrfOOVzy2IcqbZE4UmHviroq\nePMWePQY2LjYHWL+w6lQNMDrZJ7o1TWXZ68cQdf8LMZN/IBpC9d7HUkkJamwO6K5Ad77AzwwFD74\nozu0/Lo57jojvV/Kvnvk88o1R7JPz85cOamc+/65WLNHRGJMOx3bo7kBPpkEs+53h5bvc5Kbrtdz\nP6+TJZXiLjm8MP5wbp2ygAemL2H20o3ce95B9O/eyetoIinBxHI9iLKyMlteXh6zn+e5hs1uh+Ls\nB2Hreuh3GBz/SxhwjNfJkpq1ltfmruHXUz6jKRjmqmNKufrYgeRna/tAZEeMMXOstWXRnqdP0I5s\nWOTmU897HpprXUGf+xj0PyqtdijuKmMMZx/cl8NLu3PXG4t44F9LebF8NTedMogzh/YhkKHXUGRX\naAu7RXO9W6RpzhOwfCYEcuCA82D4FdD7YK/T+Vr5iipu+9tC5q/ZTP9u+Vw9ciDnDOtLdmZ6j/uL\ntGjvFnZ6F3Y4DCtnwacvwMK/QmMNFPR1U/OG/TAt5lInSjhseWvheh6esZT5azbTsyCHccNLOP/Q\nEnp1zfU6noinVNg7E2qGle/B56/D51PdTsTszvC9M+HAsdD/aMgIeJ0yZVlrmblkI4/PWs7MJZVk\nGMNxg3pw5tDenDCkh8a5JS1pDLu1rZWw/B035LHkn9C4GTJzofQ4OPF2GHQaZOd7nTItGGMYuW8x\nI/ct5stNdTz38ZdMnrOaaYvWk5cV4PghPRh9wJ4cO6gHedn6wynSWmpuYTfUuMPEl82AZe/A+vnu\n/vxusO8prqAHHgfZmm6WDEJhy8crqpj66VremL+OTbVNZGdmMGJAESP3LeaYfYvZp0dnjHb4SopK\nnyERa6FqmTtP4qqP3GXDQsBCIBtKDoMBI93WdO+hGu5IcsFQmA+XV/GvzzfwzuJKlm7YCkCvglxG\nlBZR1r+I4f2L2KdHZzI020RSRGoWdjjkynndp7Buvrus/QTqNrnHcwqg76HQb3jkcpiGOnxuTXU9\n7y6u5N0lG/l4RRUbtjQCUJCbycEle3BAn67s36eA/ft0pU9hnrbCxZf8XdjhEGxeBZWL3VodGxe7\nreb1n0FznXtORiYUD4E9D4J+h0K/EdB9UNofIp7KrLWsqqrn4xVVfLyiirmrqlmyYeu2Q+AL87PY\nv3dXvte7gL17dN52KcjN8ji5SNv8VdjrF8LC1yLlvAQ2LYVgw/bH84qgeDDseSD0OtCtjFc8CDJz\nYpZd/KmhOcTn67awYM1mPlu7mflrNrN43VaaQuFtz+nRJYeBxdsLvH/3TpQU5dOnME9zwSUp+GuW\nyMbFMPMe2KM/dN/X7RDsvq+7dNtH86Flp3KzAgztV8jQfoXb7guFLauq6li6YStLK7e66w1bee2T\nNWxpDG57njFubLxfUT799sinX1Ee/fbIZ8+uufTsmkuvglw65STHR0QEkmULO+jGJbXFLPFkrWXD\nlkZWbqpjVVUdX1bVserrOlZX1fNlVR3rtzTw7Y9D55xMehbk0KtrLj0L3KVX5LpHQQ7dOmVT1Cmb\nzjmZGj+XXeavLWwVtSSAMWZb6Q4fUPSdxxuDIdZWN7BucwPraxpYV7P99vqaBj5cVsX6mgaCO1g2\nNjszg6J8V97dOkeuO+Vsu91y6ZqXRde8LApys8jNylDJS4ckR2GLJIGczAADundiQBvLwYbDlk21\nTayvaaBySyObapuoqo1cb22iqraJjbVNrNhUS9XWJmqbQjv9WVkBQ0FuFgV5kUtuZuQ6Uup5mdsf\nz82kS24m+dmZdM7JJD87QKecTHIyVfrpRIUt0gEZGYbiLjkUd2nf/wobmkNU1TZtu9Q0NLO5vpma\n+iA1Dc3U1Ee+bghSU9/Mmur6bfc1h6IPVwYyDJ0i5Z2fHYiUeSadcjLplOPu75Qd2F70OQHyswPk\nZQXIyQqQmxkgLztAblYGuZkBcrNaHsvQH4MkpMIWiaPcrAC9C/PoXZjXoe+z1tIYDEfK3RV4bVOI\n2sYgtY1B6ppCbG0MUtcUpLYx9J371lbXUxt5rK7JPdZRxhAp8QxyswKtLhnfKfqcVkWfHcggO7PV\ndUdvt7ovM8Poj0YrKmyRJGSM2VaQPQt2fzXDUNhS3xyirjFIbVOIhmZ3qW8O0dgcdl8HQ9Q3bb/d\n0BSiIRhu9dzttxuaQ1RuCX7j+xqbQzQGw9+YUrm7jIGsQAY5rYs9UupZgQyyAobMgCv2rEAGmQFD\nZsb2+7MyjLtv2233nKyMyHXkezNbflZGy89o63u2Py+r1e/ump8V9zn/KmyRNBDIMHTOccMi8Wat\npSkUpinoLs0h626HIoXecgm1PB7efn+r79v22Lfu+8bPDluCoTDBkKWuKUgwbGkORe4LW5ojjwXD\n4W33t3xPrE85eu1xe3PDyYNi+0O/Jeq/njFmDFANDLPWTohrGhHxPWMMOZkBcjKTe92ecNjSHI4U\nemj77eZI2Qe3Xbd+3vbCb478IWj5niF7FsQ9c5uFbYwZBmCtnWaMKTXGDLPWVsQ9lYhInGVkGHIy\nAvjp2Khox+X+ALd1DbAMGBXfOCIisjPRCrsQqGr19XeOETfGjDfGlBtjyisrK2MaTkREttvtlW+s\ntROttWXW2rLi4uJYZBIRkR2IVtjVQMsxvIXApvjGERGRnYlW2C8ApZHbpcC0+MYREZGdabOwW2aE\nGGNGAdWaISIi4p2oE1qstRMTEURERNqm022IiPhETE9gYIypBFbu4rd3BzbGLEzsKFfHJWs25eoY\n5eqY3cm1l7U26jS7mBb27jDGlLfnjAuJplwdl6zZlKtjlKtjEpFLQyIiIj6hwhYR8YlkKuxknY2i\nXB2XrNmUq2OUq2PinitpxrBFRKRtSbOF3bKU604eG2OMGWWMuSmRmcS/kvX9FCXX3ZHr8YlLJH6S\nFIUdOZLypZ08tm1NbqC6rTd8HHK1+cFO1AesHTk8KaBkeX128HuT9f2001wR440xX+CWMk6YyIqb\n41v+vXbwuFfvr2i5PHt/RS4Jf72SorAjH56dvUk9WZO7nR/suH/AouXwqoCS5fXZkWR8P0HUXABX\nWmsHRp6XEJE/ItMiRzSXRr5u/bhX7682c0Uk/P0VyXFe5PUYlujPY1IUdhRR1+SOk/Z8sBPxAYuW\nw6sCSpbXp6O8ej+1R6kHW7KlbP+3W8b2xd5aePX+ipYLPHh/WWunWWuvinxZuoP1leL6evmhsL3S\nng92Ij5g0XJ4VUDJ8vqkDGvthEj5dNvJFmU8fufEVusFDQPKv/UUT95f7cgFHr6/Ir/zqh08FNfX\nKyFnM9vJGNOydv5lTNo1uVtOSmyMOdEYMyrJtiQ9l6SvT1K+nyKfkSpr7WRcph1tUcbz9w8DKpJt\nRc62cnn5/rLWTjDGvBQ5urE6+nfERkIKe1dW/DPGFEZeiBeAlsM9Y7omd5Q/JG1+sBP4AYtWMF4V\nULK8Pu2SiPfTrmiVq5ztY7EDgUcTHGWUtfbmHdzv9R+4Heby6v3Vaoy6AvfvNR6Y0OopcX29kmJI\nxBgzBiiLXLeYDvFdk7vlv13furR8gHd48gZjTGHkvnK2f9gHsuP/ssVCtBxenWQiWV6f7/Dq/RSD\nXGMjj32R4FzjW22tjopce/3+ipbLq/fXKL5ZyMu+lSuur5cOnGlD5K/4MtzOhYmR++ZYaw9p9XhV\n5PEJO/9JCcnxjccTIVleH9l1raYaVuGK6Dxr7TSv318dyJXQ91ekmMdGvjykZQdkol4vFbaIiE8k\nxZCIiIhEp8IWEfEJFbaIiE+osEVEfEKFLSLiEypsERGfUGGLiPjE/wOrOGY8C5aMmwAAAABJRU5E\nrkJggg==\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plt.plot(sigmasr,amplitudesr)\n", "plt.plot(sigmas,amplitudes)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Using lambda functions\n", "As an aside, we can use a lambda function to solve a simple equation without much hassle. For example, $\\ddot{x} + 0.1\\dot{x}+ x + 0.1 x^3 = \\sin(0.7t)$" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[ 1.47796291]\n" ] } ], "source": [ "def duff_osc2(x, v, params):\n", " omega = params['omega']\n", " t = params['cur_time']\n", " return np.array([[-x-.1*x**3-.1*v+1*sin(omega*t)]])\n", "_,_,_,a,_ = ms.hb_time(duff_osc2, sp.array([[0,1,-1]]), 0.7, num_harmonics=1)\n", "print(a)" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([ 1.47796291])" ] }, "execution_count": 20, "metadata": {}, "output_type": "execute_result" } ], "source": [ "_,_,_,a,_ = ms.hb_time(lambda x,v, params:np.array([[-x-.1*x**3-.1*v+1*sin(0.7*params['cur_time'])]]), sp.array([[0,1,-1]]), .7, num_harmonics=1)\n", "a" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Two things to note: \n", "1. Remember that the lambda function has to return an `n by 1` array. \n", "2. Time must be referenced as params['cur_time']" ] } ], "metadata": { "hide_input": false, "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.4" }, "latex_envs": { "LaTeX_envs_menu_present": true, "autoclose": false, "autocomplete": true, "bibliofile": "biblio.bib", "cite_by": "apalike", "current_citInitial": 1, "eqLabelWithNumbers": true, "eqNumInitial": 1, "hotkeys": { "equation": "Ctrl-E", "itemize": "Ctrl-I" }, "labels_anchors": false, "latex_user_defs": false, "report_style_numbering": false, "user_envs_cfg": false }, "toc": { "nav_menu": {}, "number_sections": true, "sideBar": true, "skip_h1_title": false, "title_cell": "Table of Contents", "title_sidebar": "Contents", "toc_cell": false, "toc_position": {}, "toc_section_display": "block", "toc_window_display": false }, "varInspector": { "cols": { "lenName": 16, "lenType": 16, "lenVar": 40 }, "kernels_config": { "python": { "delete_cmd_postfix": "", "delete_cmd_prefix": "del ", "library": "var_list.py", "varRefreshCmd": "print(var_dic_list())" }, "r": { "delete_cmd_postfix": ") ", "delete_cmd_prefix": "rm(", "library": "var_list.r", "varRefreshCmd": "cat(var_dic_list()) " } }, "types_to_exclude": [ "module", "function", "builtin_function_or_method", "instance", "_Feature" ], "window_display": false } }, "nbformat": 4, "nbformat_minor": 2 }