{ "cells": [ { "cell_type": "markdown", "metadata": { "toc": "true" }, "source": [ "

Table of Contents

\n", "
" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "# Coupled Nonlinear Oscillators - an example of mathematical modeling -\n", "\n", "## Modeling" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "The harmonic oscillator is given by\n", "$$\n", " \\frac{dx}{dt} = \\omega y,\\quad \\frac{dy}{dt} = -\\omega x.\n", "$$\n", "$x$ and $y$ represent the position and the velocity in the Newton equation of a mass under elastic force.\n", "\n", "If we introduce a a complex variable $\\psi = x- iy$, the equation leads to\n", "$$\n", " \\frac{d\\psi}{dt} = (i\\omega - \\gamma)\\psi,\n", "$$\n", "where we add a dumping term. Its solution is easily obtained as\n", "$$\n", " \\psi (t) = \\psi_0 e^{i\\omega t -\\gamma t}\n", "$$\n", "\n", "Consider a nonlinear oscillater with a third-power term which yields a stable amplitude $|\\psi|=1$:\n", "$$\n", " \\frac{d\\psi}{dt} = i\\omega \\psi + (1- |\\psi|^2)\\psi\\ .\n", "$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "Then, we introduce a mutual force entraining each other between two oscillators with :\n", "$$\n", " \\frac{d\\psi_1}{dt} = i\\omega_1 \\psi_1 + (1- |\\psi_1|^2)\\psi_1 + D(\\psi_2 - \\psi_1)\n", "$$\n", "$$\n", " \\frac{d\\psi_2}{dt} = i\\omega_2 \\psi_2 + (1- |\\psi_2|^2)\\psi_2 + D(\\psi_1 - \\psi_2)\n", "$$\n", "\n", "![image of coupled nonlinear oscillator](http://8tops.yamanashi.ac.jp/~toyoki/lectures/simulationMethods/images/2ctdgl_s.png)\n", "\n", "Our interest is on how these two oscillators' behavior depends on the difference of frequencies $\\delta\\omega = |\\omega_1 - \\omega_2 |$ and the interaction strength $D$. An analytical consideration leads to the following phase diagram which means the stationary state is divided into three regions.\n", "\n", "For example, this is a simplest model to describe how the cells of heart muscle could oscillate synchronously even if their own frequence is slightly different from others.\n", "\n" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "## まずは描画してみる\n", "\n", "下のサブスライドにmatplotlibのアニメーション関数を使ったプログラムを掲載してあるので、相図を参考にパラメータを変えて実行してみてほしい。\n", "\n", "パラメータは２つである。\n", " - $\\Delta \\omega$ ２つの振動子の角振動数の違い(一方は$\\omega=1$に固定）\n", " - 差が大きいほど同期して動きにくくなる。\n", " - $D$ 振動子間の結合係数\n", " - 相手を自分のペースに引きずり込む効果の大きさを表す。\n", "\n", "相図（再掲）\n", "\n", "$D$と$\\Delta\\omega$の大きさにより振る舞いは3つの領域に分けられる。(このような図を相図(phase diagram)という。）\n", "\n", "\n" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ " - matplotlibのanimationをインポート\n", " - 微分方程式の右辺関数を定義" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": true, "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "# coding: UTF-8\n", "\"\"\"\n", " TDGL model of complex variables as a model of reaction-diffusion phenomena described in chap. 15 of Ohta's textbook.\n", "\n", "\"\"\"\n", "\n", "import numpy as np\n", "from scipy import integrate\n", "from matplotlib import pyplot as plt\n", "# from mpl_toolkits.mplot3d import Axes3D\n", "from matplotlib import animation\n", "\n", "def oscillator_deriv(x, t0, omega1=2., omega2= 1, D=1.5):\n", " \"\"\"Compute the time-derivative of a coupled oscillator system. array of sets consisting of real part and imaginary part \"\"\"\n", " ampFac1 = 1.0 - (x[0]*x[0] + x[1]*x[1])\n", " ampFac2 = 1.0 - (x[2]*x[2] + x[3]*x[3])\n", " return [ampFac1*x[0] - omega1* x[1] + D*(x[2]- x[0]),\n", " ampFac1*x[1] + omega1* x[0] + D*(x[3]- x[1]),\n", " ampFac2*x[2] - omega2* x[3] + D*(x[0]- x[2]),\n", " ampFac2*x[3] + omega2* x[2] + D*(x[1]- x[3])]" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "アニメーションのための関数（メイン文から呼び出すと繰り返し実行される）" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": true, "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "def animate(i):\n", " global start_t, x0, data4plot, t4plot, plot1, plot2\n", "\n", " start_t = start_t + width_ode\n", " # integration\n", " x_t = integrate.odeint(oscillator_deriv, x0, t, args=(omega1, omega2, D))\n", " # update data for plotting\n", " data4plot = np.vstack((data4plot, x_t)) # append new columns\n", " data4plot = data4plot[step:] # delete \"step\" columns at the head\n", " # new values of the horizontal axis\n", " t4plot = np.linspace(start_t, start_t+width_plot, int(step*width_plot/width_ode))\n", " \n", " plt.xlim(min(t4plot), max(t4plot))\n", " plot1.set_xdata(t4plot)\n", " plot1.set_ydata(data4plot[:,0])\n", " plot2.set_xdata(t4plot)\n", " plot2.set_ydata(data4plot[:,2])\n", " x0 = x_t[-1] # initial value for the next calculation of ode\n", " return plot1, plot2" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "$D$, $\\Delta\\omega$をセットし、一枚目の図を作成、その後、アニメーション関数を実行" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "# parameters of the differential equation\n", "omega1 = 1.0\n", "deltaOmega = 2.5\n", "D = 1.4\n", "omega2 = omega1 + deltaOmega\n", "\n", "# setting for animation\n", "step = 5\n", "width_ode = 0.05\n", "width_plot = 10\n", "t = np.linspace(0,width_ode,step)\n", "x0 = [0.2, 0.2, 0.2, 0.2 ] # starting vector\n", "\n", "# Create new Figure and an Axes which fills it.\n", "fig = plt.figure()\n", "\n", "t4plot = np.linspace(0, width_plot, int(step*width_plot/width_ode))\n", "data4plot = integrate.odeint(oscillator_deriv, x0, t4plot,\n", " args=(omega1, omega2, D)) # ode calc. for the first plot\n", "plot1, = plt.plot([], [], 'b-')\n", "plot2, = plt.plot([], [], 'r-')\n", "plt.ylim(-1.1, 1.1)\n", "\n", "x0 = data4plot[-1]\n", "start_t =0\n", "\n", "anim = animation.FuncAnimation(fig, animate,frames=1000, interval=20)\n", "# anim.save('oscillator_xyz_t.mp4', fps=15, extra_args=['-vcodec', 'libx264']) #アニメーションを動画ファイルにする記述\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Drifting Oscillation （ドリフトモード）\n", "\n", "$D$が小さくて、振動数の差が大きいので、影響し合いながらも、それぞれの振動数に基づく動きが主である。\n", "青は$\\omega=3.5$, 赤は$\\omega=1.0$の振動子である。\n", "\n", "" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Arrest of Oscillation (振動の停止）\n", "\n", "$\\Delta\\omega >2$では、$D=1$よりも大きいと振動が減衰していく。お互いの動きを牽制しあって動けなくなる。\n", "\n", "" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Sycronized Oscillation (同期振動）\n", "\n", "お互いの中間的な振動数で同期して動く。ただし、境界に近いところでは、振幅が小さくなることにも注目してほしい。\n", "\n", "" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "## GUI for the simulation of coupled nonlinear oscillation\n", "\n", "The following is the simulation program implementing gui by using \"tkinter\", which is a somewhat classical but easy-to-use gui library. I could not run \"animation\" in matplotlib in the Tk frame, so I decided to introduce an button to update the system.\n", "The Tk create its own window and make the code run in the created window (\"inline\" statement of matplotlib does not work).\n", "\n", "[note] The library name is \"Tkinter\" for python2 and \"tkinter\" for python3." ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": true, "slideshow": { "slide_type": "subslide" } }, "outputs": [], "source": [ "# coding: UTF-8\n", "from Tkinter import *\n", "import numpy as np\n", "from scipy import integrate\n", "from matplotlib.backends.backend_tkagg import FigureCanvasTkAgg\n", "from matplotlib.figure import Figure\n", "from matplotlib.lines import Line2D\n", "\n", "class TwoTDGL:\n", " \"\"\"\n", " class definition\n", " TDGL model of complex variables as a model of reaction-diffusion phenomena described in chap. 15 of Ohta's textbook.\n", " \"\"\"\n", " pause = False\n", " \n", " def oscillator_deriv(self, x, t0, omega1, omega2, D):\n", " \"\"\"Compute the time-derivative of a coupled oscillator system. array of sets consisting of real part and imaginary part \"\"\"\n", " ampFac1 = 1.0 - (x[0]*x[0] + x[1]*x[1])\n", " ampFac2 = 1.0 - (x[2]*x[2] + x[3]*x[3])\n", " return [ampFac1*x[0] - omega1* x[1] + D*(x[2]- x[0]),\n", " ampFac1*x[1] + omega1* x[0] + D*(x[3]- x[1]),\n", " ampFac2*x[2] - omega2* x[3] + D*(x[0]- x[2]),\n", " ampFac2*x[3] + omega2* x[2] + D*(x[1]- x[3])]\n", " \n", " def __init__(self, D, delta, master):\n", " # create a container\n", " self.frame = Frame(master)\n", " self.fig = Figure()\n", " self.ax = self.fig.add_subplot(111)\n", " self.ax.set_ylim(-1.1, 1.1)\n", " self.ax.set_title(\"$\\omega_1=$\" + str(omega1) +\", $\\Delta\\omega$=\" + str(delta) + \", $D=$\" + str(D))\n", " self.line1 = Line2D([],[], color='red')\n", " self.line2 = Line2D([],[], color='blue')\n", " self.ax.add_line(self.line1)\n", " self.ax.add_line(self.line2)\n", " self.canvas = FigureCanvasTkAgg(self.fig, master=master)\n", " self.canvas.get_tk_widget().pack()\n", " self.frame.pack(side=\"top\")\n", " self.start_t = 0\n", " self.x0 = [0.2, 0.2, 0.2, 0.2 ] # starting vector\n", "\n", " def run(self, D, delta):\n", " # set parameters of the differential equation and run\n", " self.omega1 = 1.0\n", " self.omega2 = omega1 + delta\n", " self.D = D\n", " self.step = 50\n", " self.width_ode = 0.5\n", " self.width_plot = 10\n", " self.update()\n", " \n", " def reset(self, D, delta):\n", " self.start_t =0\n", " setParams(D, delta)\n", " \n", " def update(self):\n", " self.t4plot = np.linspace(self.start_t, self.start_t + self.width_plot,\n", " int(self.step*self.width_plot/self.width_ode))\n", " self.data4plot = integrate.odeint(self.oscillator_deriv, self.x0, self.t4plot,\n", " args=(self.omega1, self.omega2,self.D)) # ode calc. for the first plot\n", " self.x0 = self.data4plot[-1]\n", " self.ax.set_xlim(self.start_t, self.start_t + self.width_plot)\n", " self.line1.set_data(self.t4plot,self.data4plot[:,0])\n", " self.line2.set_xdata(self.t4plot)\n", " self.line2.set_ydata(self.data4plot[:,2])\n", " self.im = 0\n", " self.canvas.show()\n", " self.start_t += self.width_plot\n", "\n", " def setParams(self, D, delta):\n", " self.D = D\n", " self.delta= delta\n", " self.omega1 = 1.0\n", " self.omega2 = self.omega1 + delta\n", " self.ax.set_title(\"$\\omega_1=$\" + str(omega1) +\", $\\Delta\\omega$=\" + str(delta) + \", $D=$\" + str(D))\n", " \n", "# end of CTDGL class definition\n", "\n", "class controlFrame() :\n", " \"\"\"\n", " control panel for starting, restarting and changing parameter values\n", " \"\"\"\n", " def setupControlBox (self) :\n", " \n", " # control buttons\n", " self.controlBox.StartBtn = Button(self.controlBox,text=u\"start\", width=30)\n", " self.controlBox.StartBtn.bind(\"\", self.startApp)\n", "\n", " # sliders for parameter setting\n", " self.controlBox.set_D = Scale(self.controlBox, from_=0.0, to=4.0, resolution=0.01, orient=HORIZONTAL, label=\"Interaction strength: D\", length=500, command=self.changeParams)\n", " self.controlBox.set_delta = Scale(self.controlBox, from_=0.0, to=5.0, resolution=0.01, orient=HORIZONTAL, label=\"freq. difference: delta_omega\", length=500, command=self.changeParams)\n", " self.controlBox.set_D.set(self.D) # initial value\n", " self.controlBox.set_delta.set(self.delta)\n", "\n", " # pack widgets in the control frame\n", " self.controlBox.StartBtn.pack()\n", " self.controlBox.set_D.pack()\n", " self.controlBox.set_delta.pack()\n", " return\n", " \n", " def startApp (self, ev) :\n", " D = self.controlBox.set_D.get()\n", " delta = self.controlBox.set_delta.get()\n", " self.simulationFrame.run(D,delta)\n", " self.controlBox.StartBtn.config(text=u\"continue\")\n", "\n", " def changeParams(self, ev):\n", " self.simulationFrame.setParams(self.controlBox.set_D.get(),\n", " self.controlBox.set_delta.get())\n", " \n", " # constructor of control panel\n", " def __init__(self, master=None, simulationFrame=None, D=0.5, delta=2.5):\n", " self.D = D\n", " self.delta =delta\n", " self.controlBox = Frame(master)\n", " self.setupControlBox()\n", " self.controlBox.pack(side=\"bottom\")\n", " self.simulationFrame=simulationFrame\n", " \n", "# end of controlBox class\n", "\n", "if __name__ == \"__main__\" :\n", " # Create new Figure and an Axes which fills it.\n", " omega1 = 1.0\n", " delta_ini = 2.5\n", " D_ini = 1.4\n", " root = Tk()\n", " simF = TwoTDGL(D=D_ini, delta=delta_ini, master=root)\n", " app = controlFrame(master=root, simulationFrame=simF, D=D_ini, delta=delta_ini)\n", " root.mainloop()\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Inline GUI for setting parameters and controling runs\n", "\n", "\n", "More useful code using \"ipywidgets\" is shown in the other page. [Here](http://toyoki-lab.ee.yamanashi.ac.jp/~toyoki/lectures/simulationMethods/coupledNonlinearOscillationEx.html) The nootebook can be download from [here](http://toyoki-lab.ee.yamanashi.ac.jp/~toyoki/lectures/simulationMethods/coupledNonlinearOscillationEx.ipynb)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "## Excitability (興奮性）を示す非線形方程式\n", "\n", "神経膜の興奮性（閾値を超える刺激に対して平衡値から大きくずれ、ゆっくり元に戻る現象）を表す方程式としてホジキン・ハックスリー(Hodgkin-Huxley)方程式やそれを簡単化した南雲・フィッツヒュー(FitzHugh-Nagumo)方程式が有名である。\n", "\n", "以下は、それを簡単化したもので太田のテキスト（太田隆夫、非線形の物理学、裳華房）に載っている。時間があれば数値計算して振る舞いを調べてみてほしい。\n", "\n", "$$\n", "\\begin{eqnarray}\n", "\\frac{du_1}{dt} & = & f(u_1) - v_1 + I_1 + k(u_2 -u_1)\\\\\n", "\\frac{dv_1}{dt} & = & b_1u_1 -\\gamma_1v_1 + k(v_2 -v_1)\\\\\n", "\\frac{du_2}{dt} & = & f(u_2) - v_2 + I_2 + k(u_1 -u_2)\\\\\n", "\\frac{dv_2}{dt} & = & b_2 u_2 - \\gamma_2 v_2 + k(v_1 - v_2)\n", "\\end{eqnarray}\n", "$$\n", "where \n", "$$\n", "f(u) = 10 u(u-1)(0.5-u)\n", "$$\n" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "## Extension to oscillator chain described by one-dimensional differential equation\n", "\n", "次の図のように振動子が１列に並んだ系を考えよう。\n", "\n", "![1d TDGLのイラスト](http://8tops.yamanashi.ac.jp/~toyoki/lectures/simulationMethods/images/1dTDGL.png)\n", "\n", "\n", "$\\psi_n$に対する相互作用項は\n", "\n", "$$\n", " -D (\\psi_n - \\psi_{n-1}) + D(\\psi_{n+1} - \\psi_{n})\n", "$$\n", "\n", "のように書かれる。振動子間の「距離」$h$を無限小にとり、$\\psi_n$の$n$を連続変数$x$として$\\psi (x)$について方程式を考える。上記の相互作用は\n", "\n", "$$\n", " Dh^2 \\frac{1}{h} \\left(\\frac{\\psi_{n+1} - \\psi_{n}}{h} - \\frac{\\psi_n - \\psi_{n-1}}{h}\\right) \\stackrel{h\\to 0}{\\rightarrow}\n", " D\\frac{\\partial^2 \\psi (x)}{\\partial x^2}\n", "$$\n", "\n", "のように２階微分として表される。(ただし、相互作用係数を新たに$D$とした。）\n", "\n", "逆に、数値計算するときは、偏微分方程式の２階微分は、単純差分として\n", "\n", "$$\n", " \\nabla^2 \\psi \\approx \\psi_{n+1} + \\psi_{n-1} - 2\\psi_{n}\n", "$$\n", "を計算する。\n", "\n", "単純なpythonスクリプトにすると\n", "\n", "```python\n", "for i in range(N):\n", " diff[i] = psi[i+1] + psi[i-1] - 2.*psi[i]\n", "\n", "```\n", "\n", "のように書くことができるが、先週試したように非効率である。かわりに、\n", "```python\n", "diff = np.roll(psi,1) + np.roll(psi,-1) - 2.*psi\n", "```\n", "のようにnumpyの関数を利用すると高速に計算できる。roll()は配列を循環的に一つずらす（シフト）する関数である。つまり、上記のようにすると、周期的境界条件$\\psi_{N} = \\psi_{0}$を設定したことを意味する。固定周期境界や自由境界の場合は、端の設定を条件に合わせてセットする必要がある。\n", "\n", "![1d TDGLのイラスト](http://8tops.yamanashi.ac.jp/~toyoki/lectures/simulationMethods/images/partialDiff_illust.png)\n", "\n", "もし自由境界で解きたいなら\n", "$$\n", "\\psi_0 = \\psi_1,\\ \\psi_{N-1} = \\psi_{N-2}\n", "$$\n", "のように、固定境界にするには、$\\psi_0$と$\\psi_{N-1}$を一定の値に保つようセットする。\n" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([1, 2, 3, 4, 5, 6, 7, 8, 9, 0])" ] }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" } ], "source": [ "import numpy as np\n", "a = np.arange(10)\n", "np.roll(a,-1)\n" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[0.00905979 0.9580148 0.60695763 0.53766746 0.3271431 0.08697423\n", " 0.91899813 0.18152794 0.38802286 0.88343062]\n" ] }, { "data": { "text/plain": [ "array([ 1.82332584, -1.30001217, 0.281767 , -0.14123419, -0.02964451,\n", " 1.07219278, -1.56949409, 0.94396511, 0.28891284, -1.3697786 ])" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "a = np.random.rand(10)\n", "print(a)\n", "np.roll(a,1) + np.roll(a,-1) - 2.*a" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### プログラム例\n", "\n", "説明は省略するが、少し一般化し係数を複素数とした方程式\n", "$$\n", "\\frac{\\partial \\psi}{\\partial t} = (1 + ic_0)\\psi - (1+ic_2)|\\psi|^2 \\psi + (1+ic_1)\\frac{\\partial^2 \\psi}{\\partial x^2}\n", "$$\n", "を数値計算し、アニメーション表示するプログラム例を掲げる。\n", "\n", "境界条件により振る舞いが異なるので、それをメイン文でセットできるようにしてある。（独立したプログラムとするには、コマンドライン引数やGUIでパラメータをセットできるようにするのがよい。余裕があれば試してみてほしい。）\n", "\n", "あまりかっこうのよいプログラムではないが、パラメータはメイン文でセットし、animate関数のなかではglobal宣言して値をうけとっている。（グローバル変数を関数内で扱う場合の例でもある。）\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "'''\n", "Calculation and Animation for One-dimensional complex TDGL equation\n", "(solved by simple Euler method.)\n", "\n", "'''\n", "%matplotlib auto\n", "import numpy as np\n", "from matplotlib import pyplot as plt\n", "from matplotlib import animation\n", "import sys\n", "import math\n", "\n", "\n", "size = 400\n", "amp = np.zeros(size) # used in the computing process as amplitudes\n", "diff = np.array(amp, dtype=np.complex128) # used as the Laplacian term\n", "\n", "x = np.random.rand(size)\n", "x = np.array(x, dtype=np.complex128)\n", "y = np.array(x, dtype=np.complex128)\n", "site_list = range(size)\n", "\n", "def updateField(x_in, x_out, c0i =0.0 + 0.0j, ci=1 + 0.6j, Di= 1 -1.0j, boundary=\"f\"): \n", " DT = 0.002 # time step\n", " amp = np.absolute(x_in)\n", " amp = np.power(amp,2)\n", " diff[1:-1] = x_in[0:-2,] + x_in[2:] - 2.0* x_in[1:-1]\n", " x_out[1:-1] = x_in[1:-1] + DT*(c0i - ci* amp[1:-1])*x_in[1:-1] + DT* Di * diff[1:-1]\n", " # periodic boundary condition\n", " if boundary==\"periodic\":\n", " x_out[0] = x_out[-2]\n", " x_out[-1] = x_out[1]\n", " elif boundary == \"fixed\":\n", " x_out[0] =0.2\n", " x_out[-1] = 0.2\n", " else:\n", " x_out[0] = 0.2\n", " x_out[-1] = x_out[-2]\n", " \n", "fig = plt.figure(figsize=(7,4))\n", "ax = fig.add_subplot(1,1,1)\n", "ax.set_ylim(-2.0, 2.0)\n", "plt_real_x, = ax.plot(site_list, np.real(x))\n", "\"\"\"\n", "animation function\n", "\"\"\"\n", "def animate(i):\n", " global x, y, c0i, ci, Di, size, boundary\n", " \n", " # integration\n", " for count in range(20):\n", " updateField(x, y, c0i, ci, Di, boundary)\n", " updateField(y, x, c0i, ci, Di, boundary)\n", " \n", " # update data for plotting\n", " # crgb = (np.angle(x) + math.pi)/2.0/math.pi\n", " plt_real_x.set_data(site_list, np.real(x))\n", " return [plt_real_x]\n", "\n", "\n", "if __name__ == '__main__':\n", " # parameters\n", " c0 = -0.5\n", " c1 = 0.0\n", " c2 = 0.6\n", " initial_condition = 1 # 1 = random, otherwise = [0]\n", " boundary = \"free\" # \"free\", \"fixed\" or \"periodic\"\n", " \n", " c0i = 1.0 + c0*1.0j\n", " Di = 1.0 + c1*1.0j\n", " ci = 1.0 + c2*1.0j\n", " plt.title(\"$c_0=\" + str(c0) + \", c_1=\" + str(c1) + \", c_2=\"\n", " + str(c2) + \"$, boudary=\" + boundary)\n", " # initial configuration\n", " if initial_condition == 1:\n", " for n in range(size):\n", " x[n] = 0.01*(np.random.rand() -0.5 + (np.random.random() - 0.5)* 1.0j)\n", " else:\n", " for n in range(size):\n", " x[n] = 0.0 + 0.0j\n", " \n", " anim = animation.FuncAnimation(fig, animate,\n", " frames=1000, interval=10)\n", " # anim.save('ctdgl2d.c1_1.0_c2_0.6.mp4', fps=15, extra_args=['-vcodec', 'libx264'])\n", " # anim.save('oscillator_xyz_t.gif', writer='imagemagick', fps=30)\n", " plt.show()\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "パラメータの設定例\n", "\n", "- c0=0.5, c1 =0., c2=0.6, initial_condition=0, boundary=\"free\": \"励起波\"\n", "- c0=0.5, c1 =0., c2=0.6, initial_condition=0, boundary=\"free\": \"左右の波の「干渉」\"\n", "- c0=0.5, c1=2.0, c2=-0.2, initial_condition=1, boundary=\"periodic\": 長期的には「孤立波」的な振る舞い\n", "\n", "メイン文のなかのパラメータをこのような値にセットして振る舞いを確かめてみてほしい。\n", "\n", "古い資料だが、次のページにもう少し詳しい事例が紹介されている。\n", "\n", "https://toyoki-lab.ee.yamanashi.ac.jp/~toyoki/lectures/complexSys/ctdgl_pde/ctdgl_pde.html\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "anaconda-cloud": {}, "celltoolbar": "Slideshow", "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.11" }, "latex_envs": { "bibliofile": "biblio.bib", "cite_by": "apalike", "current_citInitial": 1, "eqLabelWithNumbers": true, "eqNumInitial": 0 }, "toc": { "base_numbering": 1, "nav_menu": { "height": "155px", "width": "252px" }, "number_sections": true, "sideBar": true, "skip_h1_title": true, "title_cell": "Table of Contents", "title_sidebar": "Contents", "toc_cell": true, "toc_position": { "height": "701px", "left": "0px", "right": "1130px", "top": "106px", "width": "212px" }, "toc_section_display": "block", "toc_window_display": true } }, "nbformat": 4, "nbformat_minor": 1 }