{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Duffing Oscillator Solution- State Space Form" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "ExecuteTime": { "end_time": "2017-10-02T12:54:45.688613Z", "start_time": "2017-10-02T12:54:44.737696Z" }, "init_cell": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The autoreload extension is already loaded. To reload it, use:\n", " %reload_ext autoreload\n" ] } ], "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.19348975e-15 -2.04281037e-14 -8.88178420e-16]]\n", "Constant term of FFT of signal should be zero: (-0.107710534581+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": 4, "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": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The average for this problem is known to be zero, we got 4.69398739293e-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": 6, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def duff_osc_ss(x, params):\n", " omega = params['omega']\n", " t = params['cur_time']\n", " return np.array([[x[1]],[-x[0]-.1*x[0]**3-.1*x[1]+1*sin(omega*t)]])" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[-0.00966614 0.5234018 0.85235438 0.92199611 0.73011656 0.29118919\n", " -0.27313706 -0.71945762 -0.91989617 -0.85870059 -0.53820131]\n", " [ 0.09574541 0.08002063 0.03915022 -0.01219453 -0.06042603 -0.09126401\n", " -0.09182321 -0.06181887 -0.01390861 0.03753843 0.07898057]] [[ -3.04478665e-14 7.78543896e-15 -1.23720478e-14 -2.27248775e-15\n", " -1.64035452e-14 3.11278781e-14 -9.07607323e-15 -1.67782455e-14\n", " 3.28383154e-15 3.51108032e-14 1.57235336e-14]\n", " [ 1.41884117e-15 -4.33490049e-14 -6.47579507e-11 -2.28513881e-10\n", " -2.48465137e-12 9.80704233e-14 -3.05710318e-14 4.03753073e-12\n", " 2.38120123e-10 5.53404943e-11 2.89274681e-13]]\n", "Constant term of FFT of signal should be zero: (-8.59893907551e-07+0j)\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYEAAAD7CAYAAACMlyg3AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8VNX5+PHPmUkmO4SQhbAFEkLYt4RNVBINbrTWryJY\nq9baGpdW27q2X7Xab/21xWql1g20tSq1iKi1amUJEhBFlEQ2gQAJO4FshOzrnN8fc0MCJCRkkty5\nk+f9es1rZu69J/OczMx95txz7zlKa40QQoieyWZ2AEIIIcwjSUAIIXowSQJCCNGDSRIQQogeTJKA\nEEL0YJIEhBCiB5MkIIQQPVinJAGl1KRzrJujlEpVSj10rmVCCCG6n9tJQCmVCrzTyrpJAFrrdKBE\nKTWppWXuxiCEEKJjfNz9A1rrdKVUbiur5wGrjMe5QCrQt4VlWa39/fDwcD1kyJAOxVZRUUFQUFCH\nynoKq9fB6vGD1METWD1+6N46ZGZmFmqtI9qzrdtJoA2hQHGz531bWdaqIUOGsGnTpg69eEZGBsnJ\nyR0q6ymsXgerxw9SB09g9fihe+uglDrQ3m27Ogl0iFIqDUgDiIqKIiMjo0N/p7y8vMNlPYXV62D1\n+EHq4AmsHj94bh26OgmUAGHG41CgyHjc0rJTtNaLgEUASUlJuqPZU349mM/q8YPUwRNYPX7w3Dp0\nSRJQSoVqrUuAt4EkY3EskG48bmmZEEKIbtYZZwfNAZKM+0arAbTWWcY2qUCJ1jqrpWXuxiCEEKJj\nOuPsoGXAsjOWJTZ7vKiFMmctE0II0f3kimEhhOjBJAmINhVUFnDr8lsprCo0OxQhRCfzyFNERdcq\nKq9hZ14Zu46VkltYwcnKOkqqaqmoacDf10aAr51gf19iwgIZGh7EuuKFZB3P4qUtL/HYtMfMDl8I\n0YkkCfQAdQ1Ovt5XzKe78vk0O5/cgopT6/oE+hIW5CA00EGIvw819U6KKmrJLaxgTfVtqGP1p7Zd\nmr2UpdlLcdgcZN6caUZVhBCdTJKAFztUXMmSrw+ydNNhCspqcNhtTI0N44bJgxjdvzcj+oXQN9iv\n1fJHypbz/758ii/z1lKna0D7Ulc6mor82fzg1S+5duJAZo+Lxt/X3o21EkJ0JkkCXmjXsVIWrNrD\nih3HUEByQiRzkwZyUXwEQX7tf8sHhETRL7g39boWh91BXUMd3xsXS3T9JN7LOsL972xh/vJdpF0c\ny41TBxPokI+TEFYj31ovsq+wgmdWZvPxtjyCHT78NHkY3586mAGhAR3+m8XVxcxNmMv1w6/nnd3v\nUFhVyC9Sh/PzS+P5fG8RL6zZy5Mf7+SFNXu599J4bp4Wg49dzjcQwiokCXiB6roGXsrI4aWMHHzs\niruT47j9olhCAx1u/+0FKQtOPX502qOnHiuluDA+nAvjw8k8cIJnV+3mtx/uYMlXh3ji6tFMjzvn\nuIBCCA8hP9ksqPkpmxtzi7jqL5/xl9V7uHJsPzIeTObBy0d0SgJor8SYPrz54ym8fFMi5TX1fP+V\nL3nwnS2U19S3XVgIYSppCVjQy1tfJut4Fnd+9Acysy5hUJ9A3rhtChcPb9fw4V1CKcUVY/qRnBDB\nc6v38PLaHDbuK+bmeCfJpkUlhGiLJAELSVycSG1D7ann2ZUrCR6xkgqbg4uHe8Ypm/6+dh66YgSX\njIjkl0s38/uN1dT03sNPU4ahlDI7PCHEGeRwkIUsv3Y5UyNSwekLgK/yY/bQ2ayYs8LkyM6WNCSM\n/957EVOj7Ty9cjd3/zOLCjk8JITHkSRgIWt31vD5njJQ9fjaHNTrWoIcQYQHhJsdWotC/H25Y5wf\nj1w1khXfHuO6l77gUHGl2WEJIZqRJGABWmueXpHNA+9sISyklv8ZNod/zX6LuQlzKao6a04ej6KU\n4vaLY/nHj6aQd7Kaa1/6gp15pWaHJYQwSJ+Ah3NqzWMfbGfxlwe5YfIgfnfNP/A1zsNvfsqmp7t4\neATL7pzOLX//irkLN/C3H05mytCwtgsKIbqUtAQ8WH2Dk79tq2Xxlwe5Y2Ysf7h27KkEYEXxUSEs\nu+sCIkL8uPlvG/l013GzQxKix7PuHsXL1Tc4+fmSzXx+tJ77Zw3nV1eM8IqzawaEBvDOHdMZHhXC\nnW9mkZGdb3ZIQvRokgQ8kNOpeejdrXy8LY95CQ7uuTTeKxJAo77Bfiz+8VSGRQZzx5uZfL5X5ikQ\nwiydMsewUipVKfVQC+smKaW0UirHuC00ls837tPcfX1vo7Xmtx9+y3tZR7hv1nCuHOprdkhdoneg\nL4t/MpUhfYP4yeub+GpfsdkhCdEjuZUElFKTALTW6UBJ4/NmwrTWSmsdB1wPzDeWpymlcoBcd17f\nGz2zcjevbzjA7RcN5Z5LhpkdTpcKC3Kw+CdT6R/qz49f/5rsY2VmhyREj+NuS2AeUGI8zgVSm680\nkkOjJK11407/dq113Bnre7y3Nh7k+TV7uWHyIP73qpFedQioNREhfrx+2xQCfO386LWvOF5abXZI\nQvQo7iaBUKB5O77FoSOVUqnA0maLYls7hNRTrdtdwGMfbCc5IYInrxnTIxJAo4F9AnntR5M5WVXH\nra99zb4TR2VOYyG6idJad7yw6xj/Qq11lrGjn6W1friF7ea3thxYdWaLwOgrSAOIiopKXLJkSYfi\nKy8vJzg4uENlu9PhMif/b2MVff0Vj0wLIMCnKQFYpQ6tOZ/4txfW82xmDZExH1AR8CUzgmcwr++8\nLo6wbVZ/D8D6dbB6/NC9dUhJScnUWie1Z1t3LxYrARqv+AkFWrt89VRfgbGDL9ZaLzO2jz1zY631\nImARQFJSkk5OTu5QcBkZGXS0bHcprqjlkb+uJzjAj7d/OuOsCWCsUIdzOZ/471+cSOCIWsqN5+vL\n17O+fD0Ou4PMm8wbIM/q7wFYvw5Wjx88tw7uHg56m6adeCyQDqCUCm3cQCl15k5+U+N2QJzxvEdq\ncGru/dc3FJTX8OotSW7NAOYNll+7nKuGXoUd11wIpwbIu87zBsgTwlu4lQS01llw6ph/SeNzYPUZ\nm+aeUWauUmoOkNOsTI/z51XZrN9byO++N5rxg0LbLuDlIgIjCPINwkkdSvtS66ylts7hsQPkCeEN\n3B47yDh0c+ayxGaPc4E72irT06z89hgvrMnhhsmDmDd5sNnheIzGOY0vHXA1d//nRdbszaF0Rh29\n/L3zegkhzCYDyJngYFEl9y/dwriBvXni6tFmh+NRms9p/Op3/sDchRv49bvbeP7GiT3qjCkhuosM\nG9HN6hqc3LvkG1Dw4g8m4e9rNzskj5UY04cHLkvg4215vPXVQbPDEcIrSRLoZgvSd7P5UAl/vHYc\nA/sEmh2Ox7vj4lguHh7Bbz/cIfMQCNEFJAl0oy/2FvJiRg7zkgYxe1y02eFYgs2m+PPc8YQG+PKz\nt7KormswOyQhvIokgW5SUlnLL5duZmh4EI9fPcrscCwlPNiPP8+dQE5BBU8tzzY7HCG8iiSBbvLE\nf76lqLyW526YSKBD+uPP14Xx4dx6wRD+/vk+vsiR4SSE6CySBLrB8u3H+Pfmo9xzSTxjBvQ2OxzL\neviKEcSGB/HgO1spra4zOxwhvIIkgS5WXFHLo//exuj+vbg7Jc7scCwtwGHnmbnjyTtZxf99uMPs\ncITwCpIEuthvPtjOyao6npk73tLzA3uKiYP7cHfyMJZlHpapKYXoBLJX6kIrvj3GR1vzuPeSeEb0\n62V2OF7jnkuHERcRxCPvb6eipt7scISwNEkCXaSsuo7HP/iWEf1CuDNZDgN1Jj8fO/OvG8eRkiqe\nWbnb7HCEsDRJAl3kmZW7OV5WzR+vGyeHgbpA0pAwbp4Ww2tf7OObgyfMDkcIy5K9UxfYfKiE1zfs\n55ZpMUyQ0UG7zENXJNCvlz+/encbdQ1Os8MRwpIkCXSyugYnv35vG1Eh/jxweYLZ4Xi1EH9f/u97\nY8g+XsY/Pt9vdjhCWJIkgU72xoYD7Mwr5YmrRxMiwx93uVmjorh0RCQL0ndz7KRMUi/E+ZIk0IkK\ny2tYkL6bmcMjuHx0lNnh9BiPf3c09U7Nkx/LtQNCnC9JAp2koLKA696/iaqGEn7z3VEy9n03Gtw3\nkLuTh/HR1jw+3ytDSghxPiQJdJLff/EcRQ3ZjB29kbiIYLPD6XHumBlLTN9AHvtgO7X10kksRHu5\nnQSUUnOUUqlKqYdaWT/fuE9rbxkrSVycyNjXx5J+5N8opdlbs4qxr48lcXFi24VFp/H3tfP4d0eR\nW1DB4i8PmB2OEJbhVhJQSk0C0FqnAyWNz8+QppTKwZhsvp1lLGP5tcsZG5qMdro6gf3t/sweOpsV\n160wObKeJyUhkoviw/nL6j2UVNaaHY4QluBuS2AeUGI8zgVSW9jmdq11nLHTb28Zywiw92HH4VqU\nqsdhd1DTUEOQI4jwgHCzQ+txlFI8MnskZdV1/GX1HrPDEcIS3E0CoUBxs+d9W9gm9oxDP+0pYxkv\nrNlLtfMklw74Hm9d9RZzE+ZSVFVkdlg91oh+vZg3eTBvbjhAbkG52eEI4fGU1rrjhZVaCCzUWmcp\npVKBWVrrh1vZdj6wCri+rTJG/0EaQFRUVOKSJUs6FF95eTnBwV3XSVtY5eRX66qYGu3D7eP8uuQ1\nuroOXc2M+E/WaB5eV8nIvnZ+Psnf7b9n9fcArF8Hq8cP3VuHlJSUTK11Unu2dXeKqxIgzHgcCpz2\nE9jYmRdrrZcZ62LbKgOgtV4ELAJISkrSycnJHQouIyODjpZtj/uWbsZuz+NPt1xM/9CALnmNrq5D\nVzMr/iN+e3lqeTaOQWO4IM69Q3NWfw/A+nWwevzguXVw93DQ27h27Bj36QBKqcYBczY1LgPijOct\nlrGaHUdLef+bI9w6Y0iXJQDRcbfNGMqA0ACe/GgnDc6Ot3aF8HZuJQGtdRaAcVinpPE5sLrZ+rlK\nqTlAjtY66xxlLGX+8l308vfl7pnDzA5FtMDf187DV45gR14p72YeNjscITyW2zOeG4duzlyW2Mb6\ns5ZZyRd7C1m7u4BHrhpJ70AZH8hTfXdcNH9fv49n03dz9YT++PvazQ5JCI8jVwyfJ6dT84dPdjEg\nNICbp8eYHY44B6UUD12RQN7JarmATIhWSBI4Tx9vy2PbkZPcN2u4/LK0gAviwrlwWDgvrNlLWXWd\n2eEI4XEkCZyH2nonf1qRzYh+IVwzcYDZ4Yh2evDyBE5U1vG39fvMDkUIjyNJ4Dws+fogB4srefjK\nEdhtMkqoVYwfFMoVo/vx6mf7KK6Q4SSEaE6SQDtV1zXw/Kd7mTIkjOThEWaHI87T/ZcNp7K2npcy\n9podihAeRZJAO/1z40Hyy2q477LhMleABcVHhfA/Ewfy+oYD5J2sMjscITyGJIF2aPwFeUFcX6bF\nWnqoox7tF6nxaK15brW0BoRoJEmgHd7ccIDC8lrumzXc7FCEGwaFBXLjlMEs3XSI/YUVZocjhEeQ\nJNCG8pp6Xl6bw8XDI0gaEtZ2AeHRfnrJMHxsiufXSGtACJAk0KbXv9jPico6fpkab3YoohNEhvjz\ng6kxvP/NEQ4USWtACEkC51BWXceidblcMiKSiYP7mB2O6CR3zozFx6Z4QVoDQkgSOJe/r9/Pyao6\nfpkqfQHeJLKXP9+fMph3s45wsKjS7HCEMJUkgVacrKrj1fW5XDYqirEDe5sdjuhkdyXHYZfWgBCS\nBFrz5ob9lFXXc++l0hfgjaJ6+XPjlMG8m3WYQ8XSGhA9lySBFlTW1vO39ftISYhgzABpBXirO2fG\nYVOKP3/6Nbcuv5XCqkKzQxKi20kSaMFbGw9yorKOn10iE8Z4s369/blhyiCWH15M1vEsXtryktkh\nCdHt3J5UxttU1zWwaF0u02P7khgj1wV4s8TFidQ21OLTBzSwNHspS7OX4rA7yLwp0+zwhOgW0hI4\nw7LMw+SX1UgroAdYfu1yrhp6FXYcAPjZ/Zg9dDYrrlthcmRCdB+3k4BSao5SKlUp9VAr69OM2/xm\ny+Y3rnP39TtTXYOTl9fmMGFQKBfEyRhB3i4iMIIg3yCc1KGdPtQ01BLkCCI8INzs0IToNm4lAaXU\nJACtdTpQ0vi82fpUIN2YUzjWeA6QppTKAXLdef3O9p/NRzl8oop7LhkmI4X2EMXVxcxNmMtUvydw\nnpzGsfICs0MSolu52ycwD1hlPM4FUoGsZutjjdsiY32ssfx2rfUyN1+7UzU4NS9k7GVkdC8uGRFp\ndjiimyxIWQBA9pAyLl/gQ/woOSVY9CxKa93xwkotBBZqrbOMX/mztNYPt7LtKuBhY9uHcCWLSVrr\np1rYNg1IA4iKikpcsmRJh+IrLy8nODi4Xdt+dayeFzfXcPd4P6ZEe05/+fnUwRNZKf4FmdXsKWng\nmZmB+Ps0tQStVIfWWL0OVo8furcOKSkpmVrrpHZtrLXu8A1YiGtHDq5WwPxWtpvU0jpgPpB6rtdI\nTEzUHbVmzZp2bed0OvUVC9bplKfX6PoGZ4dfryu0tw6eykrxb9pfrGMe/ki/si7ntOVWqkNrrF4H\nq8evdffWAdik27kfd7djuARoPI8yFChqZbtUbbQQjE7iOcbyIpoOEZlmTXY+O/NKuTt5mMwd3IMl\nxvRhWmwYr362j5r6BrPDEaJbuJsE3qZpJx4LpAMopUIbN1BKpWnjkI9xyGhT43ZAnPHcNFpr/vrp\nXgb2CeB7E/qbGYrwAHcnD+NYaTX//uaI2aEI0S3cSgJa6yw4tXMvaXwOrG62fL5SKkcpdaJZmblG\nayCnWRlTfJlbzDcHS7hjZhy+drlsoqe7KD6cMQN68fLaXBqcHe8vE8Iq3O4B1a7TP89clmjcpwNn\nDcTfUhmzLFqXQ3iwg+sTB5odivAASinuTh7G3f/M4pPteXxnnLQOhXfr0T99s4+VsSa7gB9OH4K/\nr93scISHuHx0P2LDg3hxTU7jCQxCeK0enQQWrcslwNfOzdNjzA5FeBC7TXHnzDh25JWydrdcPCa8\nW49NAnknq/hg8xHmTR5EaKDD7HCEh7lm4gCie/vzYkaO2aEI0aV6bBJ47fP9aODHFw41OxThgRw+\nNn584VC+2ldMTomcLiq8V49MAqXVdby18SCzx0YzKCzQ7HCEh7phymBC/H34ZF+d2aEI0WV6ZBJ4\na+NBymvqSbvY9OvUhAcL9vPhpmkxZB5vYH9hhdnhCNElelwSqKlv4O/r93HhsHCZOlK06UcXDMGu\n4NX1HjXgrRCdpsclgQ82HyW/rEZaAaJdInv5M72/D+9sOkxReY3Z4QjR6bw2CdTl59PnmT9TX9B0\nip/TqXllXS4jo3txUbxMHCLa58qhvtTUO3nzywNmhyJEp/PaJFD47FP47t1LwbNPQWkelBfw2dbd\n5OXnc/cF0Shnfae+Xl1+Pvtvuvm0pCO8Q/9gG5eOiOSNDQeoqpUzhXqyLvueO51QWwmVxVB6FE4c\ngOJ9nfsarfCcgfM7ya7xE9A1rma7Akre+4iS9z5C2TQz5+ax3R/42LjZfMEvGBwhxn1w071/LwiK\ngMBw131Q32bPw8HH77TXLXzxJaoyMyl44UWin3i8u6stutjtF8dyw6IvWZZ1mJunycWFPdU5v+c1\n5VB+HMryXPdVJ6CqBKpLoOoEow/thf1Pu5bXlEF9tetWVw0NLRxqVDZ4/ESX18nrkkDcqpXkP/UU\nZStXomvrUA4fQpKGU3/lZJ7YfIQrR0UwNaYXNNRBbbnrVtN4XwbVJ+HkEdd9ZSG01mIICIPeA9j1\nfDG6vmlogZIlSyhZsgTlcDBi65ZuqrXoalOHhjF+YG9e/SyXG6cMliHHe5hd48aja2tPPT/1Pbcr\nRtwZAGXHobas5cI+ARAQSoDTAYEDIDTG9WPTxx98A06/9/F3/cC0+4LqnqFsvC4J+EZGYgsORtfV\no319oa4e2+Bx/LFmFl84inhwziXg185qa+3K4hVFroRQUQAVha5b2VE4eYS4WxvIX11I2UEbusGG\nsjsJGVhN1KRCeG4ShMWefQsdDD5ylbKVKKVIuziOn76Vxaodx7hiTLTZIYnOVn0SCvdC0R4o3AMn\n9rkOy5QcJO7KQvK/6UXZEX/je64JibUTNas/9BsAw2ZBSBSERENwlOsWGAb+oeDrD8CmjAySk5PN\nrWMLvC4JANQXFhF6ww3kxMUSl5NL6ZFjLK85xt3JcQS1NwEAKAUBfVw3hrW4iS9gK38CvX8pyuGD\nrqvHFj8Dn0sToDjXdTu4wdXSOPV3bdB7EEQkGLcRxi0B/ELcqrvoOleM6cfgsEAWrsvl8tH9UEpa\nA5bjdELJfijY3bSzL9rruq/Ib9pO2V0/1vrEQMIV+E6NwVa5BX0oE+XwRdfVYUuci8+91j/065VJ\nYNDzfwUgOyOD6Jtu4vn3t+GbeZgfXjCkS16vMen0mTeXE28vdXUapTb7cGjtaj00JoXi3KYPXu7a\n048H9hrYLDEkEFJaDbVTwNF0ZXNdfj5H7rufgc/+GZ+IiC6pkzib3ab4yUVD+c0H37LpwAkmDwlr\nu5DoUq1+F7SG8nzI3+G6HTfuC3ZBXWXTdoF9oW88DL/MdR8e77rvM+Ss1nr90nvO/p57Aa9MAs0V\nltfwTuZhrps0gMgQ/y55jcakAxD9+G/O3kApCI5w3QZPPX1dQz2UHHB9OAt2QUG2637T36G+ikSA\nrIehbxz0GwtRYyh8f5d0Qpvk+sRBPLtqNwvX5koS8ACnOmqf+i3RNyQ17ezzd0Bls9lugyIgchRM\n+iFEjXL9yOo7zHXIpp3a/J5blNcngTe+2E9dg5OfXOShF4fZfVw7+L5xMGJ203KnE0oOsH3124yJ\nAI5tY9cTG9ENG09tcqpzysfGiMWPQNQY14db+hu6TIDDzs3Th/Dc6j3szS9nWGSw2SH1LLUVru/C\n7J+g65pO1y35cDUlH65G2TUjfjHA9V2KHNV0C5YWc2vcTgLGNJElwKTGuYTbWt9Wmc5SU69548sD\nzBoZRVyExb6sNhuEDaUwYhoYnUlxl+ST/4ffU/bpGnRNLcrXRkicg6jRx+DfdxnlfF2/dKInQP8J\nED0eIkef6pwS7rtlegwL1+bw6me5/PG6cWaH472qT0LeVsjbwsgdK2H7Q1C4G9DEXWUjf1sEZQd9\n0PUa5fAlJHkGUY88AVFRZkduKW4lAaXUJHBNI6mUilVKTWo+Z3BL6xvXtVamM607Uk9JZR13zPTQ\nVsB58o2MxNY71HXqq58furYW24Rr8PnNo1CUA8e3GV+azbDjA8h63VXQ5gORI12JIXo89J8IUaNd\np6WJ8xYe7Md1iQNZtukw9102vMsOM/YolcWuz23eFtft6GbX2TmGUEdfGDoVxlwL0ePxjR6P7ZlF\n6H1Lm74LYf3wkQRw3txtCcwDVhmPc4FUIKuN9X3bKNMp8sqO82H1X5kQcxeJMd5z7LbFTmibHSKG\nu25jrnNtqLWrr6HxC5W3GXZ9DN+86Vqv7K5DR/0nNCWHfmNP64BuJB3RZ7v9olj+9dVB3vjiAA9c\nnmB2OB6jXZ+VsuNNO/vGHf/JQ03rQ2Ncn8eJNxmfzXFs2LTjrNMr64ta+C6I8+ZuEggFips979uO\n9W2V6RSPrl2A028/UZGfAVd2xUuYot2dU0q5znDoMwRGfc+1TGs4edj1xWtMDLtXwOZ/GmVsEJ7Q\ndBgpegL0GytXQ7dgaHgQl42K4s0vD3DX+Z567MVO+6w8/hsoPWJ81rY03cqPNRXoOwwGTYEptxs/\nRMa10lm746wl3tpR292UOxNpK6UWAgu11llKqVRgltb64XOtx5UEWi1jlEsD0gCioqISlyxZ0u6Y\nfnnwl9Trs6/y9VE+PDv42fOvpMnKy8sJDu7C/gyt8aspIrg8h5Ay1y24PAe/2hPsWhqNdp59Lrz2\n8SG/2RfwXLo8/m7QWh32nmjgyY3V/GCEg1lDfE2IrP26+n2I/Nk9qPqzv3fKphkxNw+NjcrAgZSF\nxFIeHEdZSBzlwUNp8GnfpE7e/DnqCikpKZla66T2bOvuz5cSoDFthwJF7Vx/rjJorRcBiwCSkpL0\n+Vxlt7JyJU99/SdWHVhNg67F3+7PpYMv5YHJDxAeYL2RQzPMusqwNI+4y9eR/+LrlH2z39X51ng1\n9IRSRm2+x/jlNt5oNYyDkH6eE38naq0OycDyY1+w9ng1v735YnzsnjseY6e+Dw31UJht/LJ3ddzW\nXVNC/le2pitqfSBkTBRRt10DI2egokYT5AgkyBPiN4mn1sHdJPA20JhtYoF0AKVUqNa6pLX1rSzr\nFBGBEYQ4gnHqOnzwoaahhiBHkCUTgKl6ReM7fR62FTvRm/Y3db6Nugyf745q2gHs/LCpTHCUqzkf\nPa6pae9GS9MK0i6OJe3NTP67/RhXj+9vdjidr67Kde79sS1N7/nxb5sucPQJgH5j8Z0+F9uJPPSh\nzSg/X3RtHbYRKfhc9gtz4xdtcisJGId0kozDOiXNzvJZDSS2tr6VMp2muLqYuQlzGVI6hP299lNY\nVdjZL9FjtNgRfeEvmzaoLoXj20/9IuTYVsj5FLTrHO4ZPkFwYJIrIUQZ52xHjGixA7qRlTqiU0dG\nERsRxKJ1OXx3XLRHDiVxam6N0aNb/39q7RrCOH+HayefvwOObXNdvGi8l/j3dr2Pjcfvo8e7junb\nXAOd1X/inVfUmmFDThHr9xZwV/Iwgru4v8ntv24cujlzWWIb689a1pkWpCwAXM2vm6bd1JUv5fXa\n7Hzz7wUxF7hujeqqXTuRvC3kZy1nQF3hqSugXRSEDW26kKcxOYTFgd3HUh3RNpsi7aJYfvXeNr7I\nKWLGMM9rcRa++JJrbo3G/2f1Scjf2bSzb7zKtrqkqVBIf9dpxAlXNR3uC41xnXDQCumo7TwvrNnL\n7uNl3HtpfJe/lpzSIDqfrz8MmAQDJrGnfCgDkpPB2QAn9hs7np2Q/61r55P9X9BOgLM6oq0yLPc1\nEwfwzKrdvLw2x3OSgNauuTWM4Y8Vzf6fRmct4JpLI2oUjP4f104/cpTrmpLzGE5BdK7tR06yfm8h\nD18xAj/IixCmAAAaWklEQVSfrh9OWpKA6B42e9PwGKOublpeV+3qZDy+g7jRmeQv+YyyXaXoBpp1\nRB+DP49yHXpoHOCrb5zrl2noINMvevP3tXPrBUP404psdhwtZVT/Xt3zwo0DE5YccCXY4n2uK2qL\n9kDhXuKurDh9+GMfCBkVTtQtl0H8ZNfOv/egc/66F93vlc9yCXLYuXHq4G55PUkCwly+/qeOL/tO\n+D62b59A71jaNFzviBR8rhrdNOzv1qVQU3r63wiOcg37GxrTNPxvrwGus5VCol0TANlaP3OnM/og\nbpoaw4tr9rJoXQ4Lbph4zm3b/Xr1NVB2zLjluc65P3HA2Om7xrmnruL0Mr0HuZLlhO/jGz4cW+3X\n6EPr0b52qG/ANioVn+881qE6iq53+EQlH23N47YZQ+gd0D2nHUsSEB6lxY7omQ82bdA4RHBxrmsn\nWNJsp3j4a/j2/aaOzEY2Hwjud/qkH43zRASEUvj6p64+iKefJPqRh11zOvgEuAb3a6fegb7cMGUw\n//hiPw9cnsDAPs06vrU2phGsgpoyCp9+2vV6v3uA6JsvapqGsLLYmJ7Q2OlXFZ/9Qo4QV5ILi4W4\nFFfi6xPTdO84/STM+je2nja3hnTWera/rd+HAn40Y2i3vaYkAeFR2jUsd0iU6xYz/ez1DfWuWd9K\n85rmem38NV1+rGmSn6oSdr0deXofxAcrKflgZdMxc2UH3wAu0HbICnG1WnwCzjh80nQK7MMNmjk+\npQS8rMHR4OoIr6tyJQBa6PNY+RUlK79yvd4PSl1JKTjSdZV3zPSmlsyp+2jXNudx+ObMuTWE5zpZ\nWcfbXx/i6vH96R/afYc4JQkI72L3MQ4NtXE8VWvibttH/p+epmzt565RWR2+hCTFEzVnMgTZjZ14\nNQUHcxkQGda0Qz/z2gdjp+wA6qtK+bLMyazhQ3H4B7r6K3wDwdefuOma/Le/oGzTbmMQQAchKRcT\n9atfuaYoFD3a4o0HqKxt4PaLu3fAS0kComdSCt9BsdjCIk8flXXwOHyu+tVpm+7JyHCd4dQOPnml\n/PQvn/FgWAI/TTl9SlJfwLahFL1hR9PrhYbjIwmgx6uua+C1z/czc3gEI6O76cQCg+de5y5EN2js\ngxjy9hJCb7iB+kL3LiwcGd2LmcMjeO3zfVTXNZy1vrNfT3iH9785QmF5DXd0cysApCUgeriuuMDp\njotjufHVjbyXdeSs0/zkgipxJqdT88pnuYwZ0IvpcV0yqPI5SUtAiE42Pa4vYwf05tXPcmlwevfY\nScJ96TuPk1tQQdrFcaYMOyJJQIhOppTijpmx5BZWsGrHcbPDER5u0bpcBvYJ4KoxZ4/C2x0kCQjR\nBa4Y3Y9BYQEsXJeDO3N2CO+WeaCYTQdO8JMLh5o2FLkkASG6gI/dxu0XxfLNwRI2HThhdjjCQ728\nNpfQQF/mTh5kWgySBIToItcnDqJPoC8L1+aYHYrwQLuPl7Fqx3F+OH0IgQ7zztGRJCBEFwlw2Lll\n+hDSd+azN7/M7HCEh3k5I4dAh2vwQTNJEhCiC90yPQY/HxuL1uWaHYrwIIeKK/lgy1G+P2UwfYIc\npsYiSUCILtQ32I+5SYN4/5sjHC+tNjsc4SEWrcvFpuAnF3XfQHGtkSQgRBe7/aJYGpyaV6Q1IICC\nshqWbjrEdZMGEt3b3LkwQJKAEF1ucN9Arh7fn39uPEhxRa3Z4QiT/f3zfdQ1OLljZpzZoQCdkASU\nUnOUUqlKqYdaWZ9m3OY3Wza/cZ27ry+EFdydMoyqugZe+3yf2aEIE52squPNDQe4cmw0Q8OD2i7Q\nDdxKAkqpSQBa63SgpPF5s/WpQLoxsXys8RwgTSmVA0j7WPQIw6NCuGJ0P/7xxX5Kq+vMDkeYZPGX\nByivqecuD2kFACh3rmY0ftGv0lqnGzv4SVrrp5qtTwPQWi8yts0xHs/RWi87x99NA9IAoqKiEpcs\nWdKh+MrLywkODu5QWU9h9TpYPX7ovDrsP9nAExuq+W58FYeCFnNbxG30snfPsMFWfx+sHj9A8cly\nHs9UDO1l574k/y59rZSUlEytdVJ7tnX3CoVQoPkceKcNgWe0ABpNAt42Hse2lDTOKLcIICkpSSe3\ncyz3M2VkZNDRsp7C6nWwevzQuXVYU/QV6yr/hvbJZUvQFh6b1j3z/Vr9fbB6/ACPvr6KstpaHr1u\nClOGhpkdzindcpmacZgoS2udBdC441dKzVJKpRqHk4TwaomLE6m110KI6/nS7KUszV6Kw+4g86ZM\nc4MTXaquwckn++pIiunjUQkA2pEEWum8zW3sBwAaaxQKFLXyZ1K11g83+3vFxuGgIqD7Z1EQwgTL\nr13O05ue5pPcdLSqxc/uT+rgS3lg8gNmhya62Aebj1JUrflTiuf0BTRqMwmccUjnTG8DjcedYoF0\nAKVUqNa6xHic1uyXfyqwiaYO4ThgYcdCF8JaIgIjCPINAlWHdvpQSw1BjiDCA8LNDk10oQan5sWM\nvQwKsZGSEGl2OGdx6+ygxsM7xs69pPE5sLrZ8vlKqRyl1IlmZeYqpebg6ijOauFPC+GViquLmZsw\nl4FVD+NbOYPCSple0tt9tPUouQUVXB3na8qkMW1xu0+gpZaC1jrRuE8H+rSnjBA9wYKUBQCsDDlG\n2pshJE8bb3JEois1ODXPrd5DQlQIiVFnzzntCeSKYSFMkDoyihH9QnhxTQ5OmYLSa328LY+cggru\nuXQYNg9sBYAkASFMYbMp7k4Zxp78clZ8e8zscEQXcDo1f129h/jIYK4aE212OK2SJCCESWaPjSYu\nIogF6XukNeCF/rs9jz355dx7aTw2m2e2AkCSgBCmsdsUP08dTvbxMj7elmd2OKITOY2+gGGRwVw1\n1nNbASBJQAhTzR4bzfCoYBak76ZBWgNe45Ptx9h9vJx7LhmG3YNbASBJQAhT2W2KX6QOJ6eggv9s\nOWJ2OKITNLYC4iKC+M64/maH0yZJAkKY7IrR/RgZ3Yu/pO+hvsFpdjjCTZ9sP0b28TLuvTTe41sB\nIElACNPZbIpfpsazv6iS97+R1oCV1Tc4eWZVNvGRwZZoBYAkASE8wqxRUYwd0JvnPt1DnbQGLOu9\nb46QW1DB/ZclWKIVAJIEhPAISinumzWcQ8VVLMs8bHY4ogNq6hv4S/oexg/szeWjo8wOp90kCQjh\nIZITIpgwKJS/rt5DdZ1nDjEgWvfWxoMcKaniwctHeOQYQa2RJCCEh1BK8eDlCRw9Wc3iLw+YHY44\nDxU19Tz/6V4uiOvLhfHWGhVWkoAQHmTGsHAuig/n+TV7OVklcxFbxWuf76OoopYHLk8wO5TzJklA\nCA/zqytHcLKqjpcycswORbRDSWUtC9flMmtUFJMGnzVosseTJCCEhxndvzfXTBjAa5/vI+9kldnh\niDa8tDaH8pp6HrjMeq0AkCQghEe6b9ZwtIZnV+02OxRxDkdLqvjH5/u5ZsIAEvqFmB1Oh0gSEMID\nDQoL5ObpMSzLPMzu42VmhyNa8acV2QDcf9lwkyPpOLeTgFJqjlIqVSn1UCvr5xv3ae0tI4SAn6UM\nI8jhw1PLd5kdimjB5kMlvP/NEX5y0VAG9gk0O5wOcysJKKUmwalpJEsan58hTSmVgzG5fDvLCNHj\n9QlycGdyHOk789mYW2R2OKIZrTVPfrSD8GAHdyUPMzsct7jbEpgHlBiPc4HUFra5XWsdZ+z021tG\nCAHcNmMo/Xr58+THO2XiGQ/yyfZjbDpwgvsvSyDYz+2p2k3lbvShQHGz531b2CZWKZUKTNJaP9We\nMsahozSAqKgoMjIyOhRceXl5h8t6CqvXwerxg/l1+N4QzcKtJ/ndW+nMHOjbob9hdh3c5Unx1zk1\nj39WxcBgRVRFDhkZue0q50l1aK7LU5ix40cpNctIBu0pswhYBJCUlKSTk5M79NoZGRl0tKynsHod\nrB4/mF+HmVqTeXID/9lXwS+um0HvgPNPBGbXwV2eEn9BZQE3/ednFNb+D2/emnpeVwd7Sh3O1Obh\nIKVUWgu3xp15CRBmPA4FilooO8d4WgTEtlVGCHE6pRRPXD2a4spa/pK+x+xwerQFm17kaPVOhsZ/\nYbnhIVrTZkvA+FXemreBJONxLJAOoJQK1VqXAJswOoSBOGChseysMkKI1o0Z0JsbJg/mjQ37+f6U\nQcRHWfOcdKtKXJxIbUOt64mCAtYw9vWxOOwOMm/KNDc4N7nVMay1zgIwWgYljc+B1c3WzzVaAzla\n66xzlBFCnMMDlw0n0GHntx/uQGvpJO5Oy69dzvSoVLTTdSjO3+7P7KGzWXHdCpMjc5/bfQIttRS0\n1oltrD9X60II0YK+wX78ctZwfvvhDlbuOM7lo/uZHVKPEeYfzvZDNSjfehw2BzUNNQQ5gggPsP4h\nIbliWAgLuWlaDMOjgnny4x0y50A3+tdXBymuLmZK+Gzemv0WcxPmUlTlHd2Z1j7BVYgextdu44nv\njubGVzfy/Kd7LTl0sdUUldfwpxXZTIz+Ba/OnopSikenPWp2WJ1GWgJCWMwFw8K5duIAXl6bI+MK\ndYM/fLKLipp6/u97oy01Y1h7SRIQwoIemT2SEH8ffv3eNrmSuAut31PIsszDpF0c67VnZEkSEMKC\n+gb78b9XjSTzwAn+9fVBs8PxSlW1Dfz6/a3Ehgdx76XxZofTZSQJCGFRcxIHMj22L3/8ZBfHS6vN\nDsfr/HlVNoeKq/j9tWPx97WbHU6XkSQghEUppfj9tWOprXfyv+9tk2sHOtHWwyX8bf0+vj9lMNNi\nWxoSzXtIEhDCwoaGB/HQFSNYvSufd7OOmB2OV6iua+DBd7YSEeLHr68aYXY4XU6SgBAW96MLhjBl\nSBi//fBbmZO4Ezy7ajfZx8v443Xj6OXfsVFbrUSSgBAWZ7MpnpozjvoGza/elcNC7vhqXzGLPsvl\nxqmDSUmINDucbiFJQAgvMCQ8iF9dOYK1uwtYvFHOFuqI8pp67n9nM4P6BPLIVSPNDqfbSBIQwkvc\nPC2Gi4dH8ORHO8g+JheRna/ffbiDwyeqeGbueIIsPlvY+ZAkIISXsNkUz1w/nhB/H+791zcyttB5\n+GDzEd7edIi7ZsYxeUhY2wW8iCQBIbxIRIgfT18/nuzjZfz+vzvNDscS9hdW8L/vbSMxpg/3zRpu\ndjjdTpKAEF4mOSGSH184lDc2HGD59mMUVBaw4NgCCqsKzQ7N49TUN/Czf2XhY7fx3Pcn4mPvebvE\nnldjIXqAh65IYPzA3jzwzhae2vhXcmtyeWnLS2aH5XH+8N9dbD9SytPXj2dAaIDZ4Zii5/R+CNGD\n+PnYOdz756iQWpYbJwstzV7K0uylXjElYmd4N/Mw//hiPz++cCizRkWZHY5ppCUghJdaft1yJodf\n6pVTIrpr6+ESfv3+Ni6I68uvr/T+q4LPxe0koJSao5RKVUo91MK6SUoprZTKMW4LjeXzjfs0d19f\nCNGyiMAIhoSFoWz1aKcP1V40JaI7CspquOPNTCKC/Xj+xkk9sh+gObdqr5SaBKC1TgdKGp83E6a1\nVlrrOOB6YL6xPE0plQPkuvP6QohzK64uZl7CXGLL76G2eCrfHuvZ4wtV1zVw1+JMTlTWsvDmRMKC\nHGaHZDp3+wTmAauMx7lAKpDVuNJIDo2Smk0wf7vWepmbry2EaMOClAUArKxYw8vZCXyzqZTNiSVM\nGBRqcmTdz+nU3P/OFjYdOMHzN05kzIDeZofkEZQ744wYh3cWaq2zlFKpwCyt9cMtbJcKbNJalxjP\nH8KVLCZprZ9qYfs0IA0gKioqccmSJR2Kr7y8nODg4A6V9RRWr4PV4wfvqYPTEcTvNlRR06B5bFoA\nEYHWOQzSGe/Bkl01LN9fz7wEB1cO7f6B4brzc5SSkpKptU5q18Za6w7fgIW4duTgagXMb2W7VpcD\nqed6jcTERN1Ra9as6XBZT2H1Olg9fq29qw57jpfqsY8v1zOf+lQfL60yN6jz4O578Pf1uTrm4Y/0\nb/69TTudzs4J6jx15+cI14/udu3H2/wpoJRKa+GWaqwuARqvsQ4Filr5M6f6Cozyc4ynRUBsu7KV\nEMJtwyJDeO1HU8gvq+HmV7+ipLLW7JC63NKvD/HbD3dw2agofvNd75ws3h1t9gnopuP4LXkbaGxy\nxALpAEqpUN106OfMnfwmmjqE43C1JoQQ3SQxpg+v3JLEj177mh/+/Sv+efs0gr10wLT3vznMw+9t\n5aL4cJ77/kTsNkkAZ3LroKDWOgtOHfMvaXwOrD5j09wzysw1WgM5zcoIIbrJjGHhPH/jRLYfLeWW\nv23kZFWd2SF1uo+2HuX+pVuYNrQvr9yS5NXzBLvD7Z4hrfUirXV68xaD1jqx2eNcrfUdLZRZplvo\nFBZCdI/LRvfjhRsnsu3ISW585UuKymvMDqnTLN10iHv/9Q2JMX34262SAM7FOqcHCCE63RVjonnl\nliT25pczb9GXHC+tNjskty1cm8NDy7YyY1g4//jRFAId3nmoq7NIEhCih0tOiOT126aQV1LFNS98\nzrdHT5odUoc4nZo//Hcnf/hkF7PHRfPqD5N61OQwHSVJQAjBtNi+LL1zOgDXv7yBVTuOmxxR+xRU\nFnDr8lvZfyKPtDc3sXBdLjdNG8xzN0zEz0cOAbWHJAEhBACj+/fmg5/OID4ymLQ3N/H8p3twOj17\n0vqXt75M1vEs5i79P9ZkF/Dbq0fzu++NkbOAzoO0lYQQp0T28uftO6bz0LKtPL1yN1/mFvPneeOJ\nDPE3O7TTJC5OpLah6RqHKv/1BCas57lcBz+8QIbJPh/SEhBCnMbf185fbpjA/OvGsulAMVcu+IzV\nOz3r8NDSKz8kUk07NUy2n91PhsnuIEkCQoizKKWYN3kwH/7sQiJC/Pjx65u4+5+Zpp89pLXmo61H\n+cHCnRwt1ihbPQ67g9qGWhkmu4PkcJAQolXxUSH852cXsmhdDs99upd1uwv5RWo8N02L6fZz77OP\nlfH4f7bzZW4xo6J7MTLej7iwuVw//Hre2f2OzKHcQZIEhBDn5PCx8bNL4vnOuP785j/f8uTHO3nl\ns1x+dkk885IG4fDp2gMK2cfKeHlLNV+tWEevAF9+d80YbpwyGLvtolPbPDrt0S6NwZtJEhBCtMuQ\n8CDeuG0KX+QU8ueVu3ns39v56+o93DB5EPOmDGZAaAAFlQU8uO5Bnp75tFuHZuobnHy2p5C3vjrI\nqh3H8bfDTy6K5a6ZcfSRiWA6lSQBIcR5uSAunOl39mXdnkL+8fk+/rpmL8+v2cuMYeHosHfZWprF\nS1te4rFpj53X362uayDzwAnW7Mrn35uPUlheQ59AX35+aTzx+gjfuWxkF9WoZ5MkIIQ4b0opZg6P\nYObwCA4VV/LdDy9kM3VQ6lq/NHspS7OXYseX/zfxY/oE+RIa4CDIz05NvZOq2gZOVtWxv6iC3IIK\ndh0rJetgCbX1TnxsiktGRHJd4kBSEiJx+NjIyDhqboW9mCQBIYRbBoUFsur6FTy96WlWH1hNjbMG\nGw58qsZy4vAV/HTnuQcKDnLYiYsM5pZpMVwwrC+Th4QR4t/9M3/1VJIEhBBuiwiMIMg3iFpnLQ67\ng7qGOq6ZEMd9P5rDgaJKSqpqOVlZR0VtAwG+dgIcNoL9fBnSN5CIED+Z6MVEkgSEEJ2iuLqYuQmn\nn7IZ5OfDqP69zA5NnIMkASFEp1iQsuDUYzll0zo65QRfpdSkc6ybo5RKVUo9dK5lQgghup/bScCY\nWvKdVtZNAtBapwMlSqlJLS1zNwYhhBAd0xnTS6bTbA7hM8wDSozHuUBqK8uEEEKYoKsHkAsFips9\n79vKMiGEECbwyI5hpVQakAYQFRVFRkZGh/5OeXl5h8t6CqvXwerxg9TBE1g9fvDcOrSZBIwd8ply\njcNAbSkBwozHoUCR8bilZadorRcBiwCSkpJ0cnJyO17qbBkZGXS0rKeweh2sHj9IHTyB1eMHz61D\nm0nA2CGfF6VUqNa6BHgbSDIWxwKNiaOlZUIIIbqZ24eDlFJzgCSl1Byt9TJj8WogUWudpZRKMs4g\nKtFaZxllzlrWmszMzEKl1IEOhhcOWH2QcavXwerxg9TBE1g9fujeOsS0d0OltWdPJO0OpdQmrXVS\n21t6LqvXwerxg9TBE1g9fvDcOsj0kkII0YNJEhBCiB7M25PAeXdqeyCr18Hq8YPUwRNYPX7w0Dp4\ndZ+AEEKIc/P2loAlnTmekgy4JzpKBm4UbfHKJGDlD/uZA/JZbcA9pVSacZvfbJml3g8j1lQr1wFO\nfZZmGY+t9jmab9ynNVtmqffAGDBzjnEafeMyj6uD1yUBq33Yz9TCgHyWGXDP2OmkGxcYxhofdku9\nH0YdrjfineRFI99a5nNkSFNK5WB8Fyz6HvzauHYq1pM/R16XBLDeh70tVhpwL5am/3eu8dxS74fW\nOl1rfYfxNNa4mNFSdQDXTvOMoV2s9DkCuF1rHdesDpZ6D4xf/18DaK2f8uTPkTcmAat92L2G1npR\ns2FGJgGbsOj7YTTXG5OBFesQ1vYmHi32jMMmVnsPJgN9jRaAR9fBG5OAt2ltED6PZTRzs9oaEsST\naa2fAu5QSoWaHcv5aqEVABb7HBm/ntNx7Ug94hdzBxQ1GypnTlsbm8Ujh5J2k6U+7O3Q2iB8nixV\na/2w8dhS70ez47ZZuJrsaVisDrh+RcfiijnMqJNlPkdGZ3CxcTy9CFe8VnsPimjq2yvB1TLwyDp4\nY0vgbVwfGvDwD3tLmg/IB6d2Ro0dlm0OuGc2pVSa8Su6MWarvR+pnP5FzcViddBaL2s2mGOoscxK\nn6NNNP2P44znlnoPgGU0xRuKq3/AI+vglReLGb8kcnF17HnkVXreqNnprcW4dqTXa63TrfR+GId/\n5hpPExs7ia1UB2/Q2BrA9f9+qtkyy7wHzeowubFl7Il18MokIIQQon288XCQEEKIdpIkIIQQPZgk\nASGE6MEkCQghRA8mSUAIIXowSQJCCNGDSRIQQoge7P8Dpe7Ql0ldd5gAAAAASUVORK5CYII=\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "t, x, e, amps, phases = ms.hb_time(duff_osc_ss, sp.array([[0,1,-1],[.1,-.1,0]]), .1, eqform='first_order', 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": 8, "metadata": { "scrolled": true }, "outputs": [ { "data": { "text/plain": [ "[]" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD7CAYAAAB68m/qAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHdBJREFUeJzt3Xl8VOW9BvDnzb6RDFlYAoQQFmWHZMLitbbWYGsV1wiy\nI0jQq7el1YvX6r22tuoFa63SeiUsAiKLUG2ttipBS7UWyCJLWLMAAUJCtpmsk8zy3j9yAjGETCAz\nc+aceb6fD585M+9J5vcyOU/enPOec4SUEkREpH1+ahdARESuwUAnItIJBjoRkU4w0ImIdIKBTkSk\nEwx0IiKdYKATEekEA52ISCcY6EREOhHgyTeLjY2ViYmJnnxLIiLNy83NrZRSxjlbz6OBnpiYiJyc\nHE++JRGR5gkhznRnPe5yISLSCQY6EZFOMNCJiHSCgU5EpBMMdCIinWCgExHpRLcCXQiR3EVbuhAi\nTQix3HVlERHph93hmTvDOQ10IUQagB1XaUsGACllFgBTV8FPROSLjpfVYtpre5B/3uz293Ia6EpY\nF1+leSYAk7JcDCDNRXUREWne8bJazF6zD43NdoQHu/88zp7uQzcAqG73PKaH34+ISBfawjzI3w9b\nM6ZgSGy429+TB0WJiFzsRFkdZq/Zh0B/4bEwB3oe6CYA0cqyAUBVxxWEEBlCiBwhRE5FRUUP346I\nyLudKKvDrDV7EegvsC1jqsfCHLjOQBdCGJTF7QCSlOUkAFkd15VSZkopjVJKY1yc04uFERFplpph\nDnRvlks6AKPy2GY3AEgp85R10gCY2p4TEfma1t0srWG+dYnndrO05/Swq5RyJ4CdHV5Labec6Ya6\niIg0oy3M/f1awzwpLkKVOnhQlIioBwrKL4f5tgz1whxgoBMRXbfCi/WYtWYf/PxaZ7OoGeYAA52I\n6LoUVdRj1pq9AICtS6ZgqMphDjDQiYiu2anKBsxesxcOh8SWJZMxrI/6YQ54+J6iRERad6aqAbMy\n98Jqbw3zEX17qV3SJRyhExF109nqRszK3AuLzY7Niyfjxn6Rapf0LQx0IqJuOFfTiIcy96KhpTXM\nR8V7V5gDDHQiIqdKTU2YtWYvai1WbF48GWMGRKldUqcY6EREXbhgbg1zU0NrmI8d6J1hDvCgKBHR\nVZXXWjB7zT5U1bdg0+JJGD/I4PyLVMQROhFRJy7WWjBrzV5crLVg46JUJCf0VrskpzhCJyLqoKKu\nGbPX7kOZ2YKNiyYhZXC08y/yAhyhExG1U1XfjDlr9+J8TRPWL0xFaqI2whxgoBMRXVLd0II5a/eh\npLoR6xYaMSVJW3fVZKATEQEwNbaG+anKBqydn4qbhsaqXdI14z50IvJ55kYr5qzdh6KKeqydb8TN\nw7UX5gBH6ETk48xNVsxbvw8F5fVYPS8Ft4zQ7q0yGehE5LPqLFYsWL8fxy7U4v/mJuPWG/qoXVKP\nMNCJyCc1ttiwaEM28s+b8fvZybhtZF+1S+oxBjoR+RyL1Y5HNuYg90wNfvfQBPxgdD+1S3IJHhQl\nIp/SbLPj0c25+FdxFX47YzzuGhevdkkuwxE6EfkMq92BJ7Z8g7+fqMDL943FfRMHql2SSzHQicgn\n2OwOLNt2ALuOluOFe0bjoUkJapfkcgx0ItI9u0PiP3cewseHL+DZH43E/KmJapfkFgx0ItI1h0Pi\n2Q8O44NvzuOp20dgyS1JapfkNgx0ItItKSV++Zcj2JZ9Fv/x/WF44vvD1S7JrRjoRKRLUkq89Ndj\n2PivM8i4JQk/mzZC7ZLcjoFORLr0210nsebLU1gwdTCeueNGCCHULsntGOhEpDurdhdg1eeFmDVp\nEJ6fPtonwhxgoBORzmT+owiv7jqJ+ycOwIv3joWfn2+EOcBAJyId2fj1abz01+O4a1x/rEwf51Nh\nDjDQiUgntu0vwfMfHsHto/ritZkTEODve/Hmez0mIt15P+8cnvngML53QxxWzZ6IQB8Mc4CBTkQa\n99GhUjy14yBuGhqDt+amIDjAX+2SVOP0aotCiHQAJgDJUsqVXbQnSSkzXV8iEVHnPjtShp9sOwDj\n4GismW9ESKDvhjngZIQuhEgGACllFgBT2/MO7cVKe3HHdiIid/myoAJPbPkGYwdEYf3DqQgL4tXA\nne1ymYnW0TcAFANI62SdFcpjkpQyz1WFERFdTfbpamRsysXQPhHY+PAkRAQzzAHngW4AUN3ueUz7\nRiXAi4UQNR3Wu0QIkSGEyBFC5FRUVPSoWCKiw+fMWPR2NvobQvDO4kmICgtUuySv0aODokIIA1pH\n8C8DWCOEuOIyZlLKTCmlUUppjIvT7t20iUh9J8vrMH/9PkSGBmLz4smIjQhWuySv4uzvFBOAaGXZ\nAKCqQ3sGgJellCYhRDGAdABXHDglIuqpM1UNmLt2HwL9/bBlyWTEG0LVLsnrOBuhbwfQNupOApAF\nXBqZf4uUcicu728nInKZUlMTZq/ZB6vdgc2PTMbgmHC1S/JKXY7QpZR5QgijECINgKndQc/dAFKk\nlCuFEMuV0Xk0py0SkatV1jdj7tp9qG2yYsuSKRjRt5faJXktp4eGOwtpKWVKu2XuYiEitzA3WjFv\n3X6UmpvwzuLJGDswSu2SvBrPFCUir1TfbMOCt/ej6GI91sw3IjUx2vkX+ThO3iQir2Ox2vHIxmwc\nPm/Gm3OS8Z3hnCHXHRyhE5FXabE58NjmXOw7VY1XHxyPH4zup3ZJmsFAJyKvYbM78NPtB/DFiQq8\neO9Y3DtxgNolaQoDnYi8gsMh8V/vH8bHhy/g2R+NxOzJCWqXpDkMdCJSnZQSL3x0FDtzz+HHtw3H\nkluuOOmcuoGBTkSqe23XSWz4+jQW3zwEP00brnY5msVAJyJVrf/qFN74vBAzjAPx3J0jIYRv3QfU\nlRjoRKSaP+aewwsfHcUPRvfFS/eNZZj3EAOdiFSx62g5lv/xEG4aGoPXH5rokzd1djX+DxKRx+0t\nrsLjW/IwJj4Smbx1nMsw0InIo/LPm/HIxhwkRIfhbd5tyKUY6ETkMcUV9Viwfj+iQgPxzuJJiA4P\nUrskXWGgE5FHXDA3Yd66/QCAdxZPQv8o3qDC1fi3DhG5XXVDC+at2w9zkxXbMqYgKS5C7ZJ0iSN0\nInKr+mYbHn57P85WN2LtAiPGDOA1zd2FI3Qicptmmx1L38lBfmktVs9NwZSkGLVL0jWO0InILewO\niWXbDuCfhVV4JX0c0kb1Vbsk3WOgE5HLSSnxP3/Ox9/yy/Dfd43C/ckD1S7JJzDQicjlVn1eiHf3\nleDR7w7F4puHqF2Oz2CgE5FLbdtfgt/uOon7kwfg6R/eoHY5PoWBTkQus+toOX7+wWF8d0QcVjww\njhfb8jAGOhG5RO6ZajyxJQ9jB0ThzTnJCOTFtjyO/+NE1GOFF+uwaEMO+keFYP3CVITz+iyqYKAT\nUY+UmS2Yv24/Av39sGnRZMREBKtdks9ioBPRdTM3WrFg/X7UWmzY8HAqEmLC1C7JpzHQiei6WKx2\nLNmUg+LKeqyel8JT+r0Ad3QR0TVrOwt0/+lqvDFrIv5tWKzaJRE4QieiaySlxPMf5uOTI61ngd49\nPl7tkkjBQCeia/KHLwqxeW8Jlt6SxLNAvQwDnYi67f28c/jNZydx38QBePqHN6pdDnXAQCeibvm6\nsBJP//EQpibFYMUD4+Dnx7NAvQ0DnYicOlleh6Wbc5EYE4635qUgKIDR4Y2cznIRQqQDMAFIllKu\n7KQ9GUASAEgpd7q8QiJS1cVaCx5+Oxshgf54++FURIUGql0SXUWXv2aVsIaUMguAqe15B88oQZ50\nlXYi0qiGZhsWbcxGTWML3l6YioG9eeKQN3P2d9NMtI7OAaAYQFr7RmX0ng0AUsqVUso8l1dIRKqw\n2R14YksejpbW4g+zk3nikAY4C3QDgOp2zzveEDAVQIwQIlkIsbyzbyCEyBBC5AghcioqKnpQKhF5\nipQS//PhEXxxogK/uncMbr2xj9olUTe44shGVdvIXBmxf4uUMlNKaZRSGuPi4lzwdkTkbm/tKcaW\nfSV47HtDMWfyYLXLoW5yFugmANHKsgFAVYf2KrTuimlbN9V1pRGRGj48WIoVnxzH9PHx+M/becch\nLXEW6NuhzGBRHrMAQAhhUF7b2a7dAGV/OhFp077iKjz13kFMGhKN3zzIueZa02Wgt9uVkgbA1O6g\n526lvRits1/SAcRw2iKRdp2qbMDSzbkYGB2KzHkpCA7wV7skukZO56FLKTM7eS2lk3aGOZFGmRpb\nsGhDNvyEwIaFk2AIC1K7JLoOPN2LyMe12Bx4dHMuztc0IXNeCm9SoWG8HjqRD5NS4rk/Hcbe4mr8\nbuYEGBOjnX8ReS2O0Il8WOY/ivFezjn8+PvDcO/EAWqXQz3EQCfyUZ8eKcP/fnIcd47rj2VpI9Qu\nh1yAgU7kg/LPm7Fs2wGMH2jAqw+O5/REnWCgE/mYMrMFizdmIzo8CJnzUxASyOmJesFAJ/IhjS02\nLN6YjXqLDWsXGNGnV4jaJZELcZYLkY9wOCSWbTuAYxdqsW5BKkb2j1S7JHIxjtCJfMSKT4/js6Pl\n+O+7RvHqiTrFQCfyAduzS7B6TzHmTknAwpsS1S6H3ISBTqRz/yqqwrMf5OM7w2Pxi+mjIQRntOgV\nA51Ix85WN+Kxd3ORGBuOP8xJRoA/N3k946dLpFMNzTYs2ZQDh0Ni7XwjIkN4c2e94ywXIh1yOCR+\n9t4BnCyvw4aHJyExNlztksgDOEIn0qE3Pi/Ap0fK8fMfjcQtI3jrR1/BQCfSmU/yL+B3WQW4P3kA\nFt88RO1yyIMY6EQ6crysFj977yDGDzLgpfvGckaLj2GgE+lEdUMLlmzKQURwADLn8RotvogHRYl0\nwGp34PF381Be24ztGVPQN5LXaPFFHKET6cCLHx/Dv4qr8PJ9YzExobfa5ZBKGOhEGrc9uwQbvj6N\nxTcPwQMpA9Uuh1TEQCfSsNwz1XjuT62n9T9zx41ql0MqY6ATadQFcxOWvpOHeEMoVs2ayNP6iQdF\nibTIYrUjY1Mumlps2LpkMgxhQWqXRF6AgU6kMVJKPP3HQ8gvNSNznhHD+/ZSuyTyEvwbjUhjMv9R\njD8fKMWT00Zg2qi+apdDXoSBTqQhX5y4iP/95DjuHNsfj986TO1yyMsw0Ik0oqiiHj/e+g1G9ovE\nKw+O42n9dAUGOpEG1FqsWLIpB4H+fsicn4KwIB7+oisx0Im8nN0h8ZOt36CkqhH/NycZA3uHqV0S\neSn+mifycq98egJfnKjAr+8dg8lJMWqXQ16MI3QiL/bnA+fx1p4izJ6cgLlTBqtdDnk5p4EuhEgX\nQqQJIZY7Wa/LdiK6NofPmbF85yFMSozGL6aPVrsc0oAuA10IkQwAUsosAKa2552slwZgmuvLI/JN\nF+ssyHgnB7ERwXhzbjKCAvjHNDnn7KdkJgCTslwMIM295RBRs82OxzbnoaaxBavnpSA2Iljtkkgj\nnAW6AUB1u+dXHJERQiQrI3gi6iEpJZ7/8xHknqnBK+njMWZAlNolkYa44u+46K4ahRAZQogcIURO\nRUWFC96OSL/e2XsG27LP4vFbh2L6+Hi1yyGNcRboJlwObAOAqvaN3RmdSykzpZRGKaUxLi7u+isl\n0rmviyrxy78cRdrIPnhy2g1ql0Ma5Gwe+nYARmU5CUAWAAghDFJKE4AkIUQSWkM/Wgn4PLdVS6RT\nZ6sb8fi7eRgSG47XZk6Anx9P66dr1+UIvS2clVkspnZhvVtp3yml3Km8ZnBblUQ61tBsw5JNObA7\nJNbMN6JXSKDaJZFGOT1TVEqZ2clrKZ2sc8V6RNQ1h0PiqR0HcbK8DhsenoQhseFql0QaxsmtRCp6\nfXcB/pZfhmfuGIlbRvAYE/UMA51IJR8fuoDXdxcgPWUgHvnOELXLIR1goBOpIP+8GU/uOIDkBANe\nvG8Mr21OLsFAJ/KwirpmZGzKQe+wILw1LwXBAf5ql0Q6wcvnEnlQs82ORzfnorqxBTsfvQl9eoWo\nXRLpCAOdyEOklHjug3zknqnBH2Yn87R+cjnuciHykHVfncKO3HP48W3Dcee4/mqXQzrEQCfygD0n\nK/DSX4/hjjH9sOy24WqXQzrFQCdys6KKejyxJQ839IvEqzPG87R+chsGOpEbVTe0YPGGbAT5+2HN\n/BSEBfGwFbkPf7qI3MRitSNjUw5KzRZsXTIFA3uHqV0S6RxH6ERu4HBILN95CDlnavDbGeORMri3\n2iWRD2CgE7nBa1kn8eHBUiz/4Q24axxvVEGewUAncrH3cs5i1eeFeCh1EB777lC1yyEfwkAncqGv\nCyvx8/cP4+ZhsfjVvbxGC3kWA53IRQov1mHp5lwMiQ3Hm3OTEejPzYs8iz9xRC5QXmvBgvXZCA7w\nx/qFqYjkXYdIBQx0oh4yN1oxf91+mBpbsH6hEYOiOT2R1MF56EQ90NRix+KN2SiurMfbCydh3EDe\nWpfUw0Anuk42uwNPbMlDbkkNVs2aiJuHx6pdEvk47nIhug5SSjzz/mHsPn4RL9w9mnPNySsw0Imu\nkZQSv/74GHbknsNPbhuOeVMT1S6JCAADneiaSCmx4pMTWPfVKSy8KRHL0ngpXPIeDHSia/BaVgHe\n2lOEOZMT8Pz0UTxxiLwKA52om1btLsAbuwswwzgQv7qHZ4GS92GgE3XD6j1FeHXXSdw/cQBevn8c\nb1JBXonTFom6IKXEG7sL8VrWSdw1rj9Wpo+DP8OcvBQDnegqpJR48eNjWPvVKdyfPAArHxiHAF6f\nhbwYA52oE3aHxLMfHMa27LNYMHUwnp8+mrtZyOsx0Ik6aLE58LP3DuCjQxfwxK3D8OTtI3gAlDSB\ngU7UTq3FisffzcOXBZV45o4bsZQ3qCANYaATKc5WN2LRhmycqmzAygfGYUbqILVLIromDHQiAHkl\nNcjYlIMWmwObFk/CTUN5oS3SHgY6+by/HCzFkzsOon9UCLYvTcXQuAi1SyK6Lk4DXQiRDsAEIFlK\nubKT9gxlcaiU8mkX10fkNla7A7/57ARW7ylGamJvrJ5nRHR4kNplEV23LifVCiGSAUBKmQXA1Pa8\nXXsagCwpZSaAJOU5kdcrM1swe81erN5TjLlTErD5kckMc9I8ZyP0mQB2KcvFANIA5LVrT1L+ZSrt\nSa4ukMjVviyowLJtB9BkteP1hybgngkD1C6JyCWcBboBQHW75zHtG5WReZtkANs7fgNll0wGACQk\nJFxflUQu0GJz4PefF2DVF4UY3icCb85JxrA+vdQui8hlXHJQVNkVkyelzOvYpoR+JgAYjUbpivcj\nulZHS2vx1I6DOHqhFukpA/HCPaMRFsQ5AaQvzn6iTQCilWUDgKqrrJfGA6Lkjax2B976exHe+LwA\nUaGByJyXgttH91O7LCK3cBbo2wEYleUkAFkAIIQwSClNynJG2+wXIUSacgCVSHX558145v3DOHze\njOnj4/HLu0fzwCfpWpezXNp2oSizV0ztdqnsbvf6CiFEkRCixq2VEnVTTUMLnv3gMKb//iuUmprw\n5pxkrJo1kWFOuud0J2KHA59tr6Uoj1kAeruhLqJrZndIbN1fgt98dgJ1FhsWTE3ET6eNQFRooNql\nEXkEjwqR5kkp8fcTFXjl0xM4eqEWU5Ki8Yu7R+PGfpFql0bkUQx00rSviyrx6mcnkXumBgnRYVg1\nayLuGtefl7sln8RAJ82RUiL7dA1e330S/yysQr/IELx43xjMMA5CIO8oRD6MgU6aYXdIfJJfhswv\ni3HwrAkx4UF47s6RmDtlMEIC/dUuj0h1DHTyerUWK97PPYd1/zyFs9VNGBwThl/dMxrpKYMQGsQg\nJ2rDQCevJKXEgbMmbN1fgr8cvIAmqx0TEwx49kcjMW1UP/jz/p5EV2Cgk1epqGvGR4dK8V7OORy7\nUIuwIH/cMyEeD01KwPiBUTzYSdQFBjqprtZixaf5ZfjwYCn+WVgJhwRGx0fi1/eOwT0T4tErhPPI\nibqDgU6qqKhrxu5j5fjsaDm+KqxEi82BQdGh+PfvDcPdE+Ixoi+vgkh0rRjo5BEOh8SR0lp8WViB\n3ccuIq+kBlICAwyhmDM5AdPHx2PiIAN3qRD1AAOd3OZcTSO+KqjEl4WV+LqwEjWNVgCtu1OW3TYC\n00b1xcj+vRjiRC7CQCeXcDgkiirqkXOmBjmna5B7phqnqxoBAH0jg/H9G/vi5uEx+LdhsejTK0Tl\naon0iYFO16Wyvhn5583IP29G7pka5JWYYG5qHYHHhAcheXBvzJ+aiO8Mj8WwPhEchRN5AAOduiSl\nxAWzpTW8S2txtNSM/PO1KKu1XFpnWJ8I3DGmH1IG94YxMRqJMWEMcCIVMNAJQOsuk1JzEwou1qPo\nYj0KyutRWFGPgvI61FpsAAA/AQyNi8DUoTEYHR+J0fFRGBUfycvTEnkJBroPsTskLpibUFLdiLPV\njThT1YiS6kacrmpA0cUGNFntl9aNCQ/CsD4RmD4+Hjf064UxA6Iwsl8kT7Un8mIMdB2xWO0or7Xg\ngtmCMnPr43lTI0qqm3C2uhHnahphtV++T3eAn8CA3qFIiA7DQ5OiMbxPLwzrE4FhfSJ4dx8iDWKg\nezkpJRpa7Kiqb0ZlfTMq61taH+tacLHucnCX1VpQ3dByxddHhQYiIToMo+Ij8cMx/ZAQHYbB0WEY\nFB2G/lEhCODlZol0g4HuQW3hbGpsgbnJCnOjFeYmK0xNVpjalhtbA7uivuVSiFusjk6/X++wQPSP\nCkX/qBBMTDCgf1QI+inP+0WFoF9kCMKD+RET+Qpu7d3gcEhYbHbUW2yob7ahodmO+ua25cuPDc02\n1F1abl2n1mL9VnjbHPKq7xPoL2AIC0JsRDBiI4KQFBuO2IggxEQEX3otVlmODg9CUABH10R0mSYC\n3dTYgjqLDVa7A1a7VB5bl212B1rsDtjaXndIWG0O2BwOtCjtbes22xxottphsdphsTrQ1LZsc8Bi\ntStt7V5X2lpsnY+QOxICCA8KQHiwP8KDAxARHIBeIQGIjwpFVFggDKGBiAoNhCEsEFGhQe2WWx9D\nA/053Y+IrpsmAv3lvx7H9pyzLvlewQF+CAn0R0ig8hjgj5Agf4QE+MEQFvTt1wP9lDZ/hAa1hbQ/\nwoMCEBHSGthtwR0eHICwQH/48TrdRKQSTQR6unEgjIm9ERTghwA/PwT6CwT6+yHQ3w8Bl5a//Rig\nLAcpywF+rcsMXCLSK00EempiNFITo9Uug4jIq/GoGhGRTjDQiYh0goFORKQTDHQiIp1goBMR6QQD\nnYhIJxjoREQ6wUAnItIJIeXVLxbl8jcTogLAmXYvxQKo9FgBnqG3PumtP4D++qS3/gD661NP+zNY\nShnnbCWPBvoVby5EjpTSqFoBbqC3PumtP4D++qS3/gD665On+sNdLkREOsFAJyLSCbUDPVPl93cH\nvfVJb/0B9NcnvfUH0F+fPNIfVfehExGR66g9QicvI4RI7qItXQiRJoRY7smaesJJf1Yojxmeq4jI\nfTwW6M7CQGth0Y3+aC4shBBpAHZcpS0ZAKSUWQBMXQWlt+iqP4oMIUQRgGIPldRjQogM5d+Kq7Rr\nbTty1h9NbkfKP49/Rh4JdGdhoLWw6Ga9mgsLpT9Xq3cmAJOyXAwgzSNF9YCT/gDAEinlUGU9r6f8\ngsqSUmYCSFKet2/X2nbUZX8UmtqOlD48qHwGyZ7OOk+N0J2FgdbCojv1aiosusEAoLrd8xi1CnGh\nJC2NZgEk4fLPWrHyvD2tbUfO+gNobDuSUmZJKZcqT5OklHkdVnHrZ+SpQHcWBloLi+7Uq7Ww8DlS\nypVKUMRcZXToVaSUmcpoFgCSAeR0WEVT21E3+gNodDtS6l3aSZNbPyMeFHUTrYVFN5gAtN3Y1QCg\nSsVaekzZb5uuPK1C56NDr6T8mZ7XyehPk7rqj1a3IynlSgBLhRAGT76vpwLdWRhoLSy6rFfLYdFR\nux/I7bjcjyQAmvgTuKN2/cnB5T4MReejQ2+VJqV8upPXtbYdtem0P1rcjoQQ7febFwPoeDDXrZ+R\npwK90zDQcFg4648mw0LZeIztNiIA2A0AbaMnZZRk0sLosBv9maG0FWmhP0BryCmjv7bPQsvbkbP+\naHE7SsO3A7sY8Nxn5LETi5RpR8VoPVCQqbyWK6VMuVq7N+tmf6qV9pXqVUp60W4aZjVaQ+NBKWWW\nVreja+iPZrYjJbhnKE9T2g6Qeuoz4pmiREQ6wYOiREQ6wUAnItIJBjoRkU4w0ImIdIKBTkSkEwx0\nIiKdYKATEenE/wOVAHOiFDeWZQAAAABJRU5ErkJggg==\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(ms.duff_osc, 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": [] }, { "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" ] } ], "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.2" }, "toc": { "nav_menu": {}, "number_sections": true, "sideBar": true, "skip_h1_title": false, "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 }