{
"cells": [
{
"cell_type": "markdown",
"metadata": {
"toc": "true"
},
"source": [
"# Table of Contents\n",
" "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"
\n",
"\n",
"\n",
"描画(CG)\n",
"
\n",
"
\n",
"\n",
"file:/~/python/doing_math_with_python/cg.ipynb\n",
"
\n",
"cc by Shigeto R. Nishitani 2017 \n",
"
"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# CG(ComputerGraphics)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## listplot, pointplot \n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"リスト構造にある離散的なデータはlistplotで表示してくれる．listplotは受け取ったlistの要素をy値に，1から始まる添字をx値にして，デフォルトで\n",
"は線でグラフを書く．\n",
"```maple\n",
"> T:=[seq(exp(-i),i=0..5)]; \n",
"> listplot(T);\n",
"```\n",
"$$\n",
"T\\, := \\,[1,\\exp(-1),\\exp(-2),\\exp(-3),\\exp(-4),\\exp(-5)]\n",
"$$\n",
"\n",
"|{{attach_view(MapleCGplot2d1.png,)}}|\n",
"|:----|\n",
"\n",
"\n",
"以下のようにoptionをつけるとpointで描く．\n",
"```maple\n",
"> listplot(T,style=point):\n",
"```\n",
"\n",
"それぞれの値の横軸xが1,2,3,..では不都合なときには，2次元のlistlist構造を用意し，[x[i],y[i]]を入れてpointplot関数で表示する\n",
"．\n",
"```maple\n",
"> T:=[seq([i/2,exp(-i/2)],i=0..6)]; \n",
"> pointplot(T,symbol=circle,symbolsize=20);\n",
"```\n",
"$$\n",
"T\\, := \\,[[0,1],[1/2,\\exp(-1/2)],[1,\\exp(-1)],[3/2,\\exp(-3/2)], \\notag \\\\\n",
"[2,\\exp(-2)],[5/2,\\exp(-5/2)],[3,\\exp(-3)]]\n",
"$$\n",
"\n",
"|{{attach_view(MapleCGplot2d2.png,)}}|\n",
"|:----|\n",
"\n",
"\n",
"\n",
"listplotのように線でつなぎたい時には，以下のようにoptionをつける．\n",
"```maple\n",
"> pointplot(T,connect=true):\n",
"```\n",
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 写像の表示 \n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"ある行列によって点を移動させる写像の様子を示すスクリプトを通して，plottoolsが提供するdisk, arrowの使い方を示す．先ず描画に必要なライブラリーパッケージ(plotsおよびplottools)をwithで読み込んでおく．\n",
"```maple\n",
"> restart; with(plots):with(plottools):\n",
"```\n",
"$$\n",
"$$\n",
"行列$A= \\left[ \\begin {array}{cc} 1&2\\\\ 2&1\\end {array} \\right]$\n",
"によって点$a_0$(1, 2)が$a_1$(5, 4)に移動するとする(LinearAlgebra参照）．\n",
"```maple\n",
"> with(LinearAlgebra): A:=Matrix([[1,2],[2,1]]): a0:=Vector([1,2]): a1:=A.a0;\n",
"```\n",
"$$\n",
"{\\it a1}\\, := \\, \\left[ \\begin {array}{c} 5\\\\ 4\\end {array} \\right]\n",
"$$\n",
"ベクトルが位置座標を意味するようにlistへ変換(convert)する．\n",
"```maple\n",
"> p0:=convert(a0,list):p1:=convert(a1,list):\n",
"```\n",
"位置p0に円(disk)を半径0.2,赤色で描く．同じように位置p1に半径0.2，青色でdiskを描く．\n",
"```maple\n",
"> point1:=[disk(p0,0.2,color=red),disk(p1,0.2,color=blue)]:\n",
"```\n",
"もう一つ，p0からp1に向かう矢印(arrow)を適当な大きさで描く．後ろの数字をいじると線の幅や矢印の大きさが変わる．\n",
"```maple\n",
"> line1:=arrow(p0,p1,0.05,0.3,0.1):\n",
"```\n",
"これらをまとめて表示(display)．このとき，表示範囲を0..6,0..6とする．\n",
"```maple\n",
"> display(point1,line1,view=[0..6,0..6],gridlines=true);\n",
"```\n",
"\n",
"|{{attach_view(MapleCGplot2d3.png,)}}|\n",
"|:----|\n",
"\n",
"$a_0$(1, 2)の赤点が，$a_1$(5, 4)の青点へ移動していることを示している．\n",
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 回転写像 \n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"次に原点周りでの回転の様子を示す．回転の行列．\n",
"```maple\n",
"> Matrix([[cos(theta),sin(theta)],[-sin(theta),cos(theta)]]);\n",
"```\n",
"$$\n",
"\\left[ \\begin {array}{cc} \\cos \\left( \\theta \\right) &\\sin \\left( \\theta \\right) \\\\ -\\sin \\left( \\theta \\right) &\\cos \\left( \\theta \\right) \\end {array} \\right]\n",
"$$\n",
"これを関数のように定義している．\n",
"```maple\n",
"> A:=t->Matrix([[cos(t),sin(t)],[-sin(t),cos(t)]]);\n",
"```\n",
"$$\n",
"A\\, := \\,t\\mapsto \\left[ \\begin {array}{cc} \\cos \\left( t \\right) &\\sin \\left( t \\right) \\\\ -\\sin \\left( t \\right) &\\cos \\left( t \\right) \\end {array} \\right]\n",
"$$\n",
"tに回転角(Pi/3)を入れている．\n",
"```maple\n",
"> a0:=Vector([3,0]);\n",
"> a1:=A(Pi/3).a0;\n",
"```\n",
"$$\n",
"{\\it a0}\\, := \\, \\left[ \\begin {array}{c} 3\\\\ 0\\end {array} \\right] \\notag \\\\\n",
"{\\it a1}\\, := \\, \\left[ \\begin {array}{c} 3/2\\\\ -3/2\\,\\sqrt {3}\\end {array} \\right] \\notag\n",
"$$\n",
"表示の仕方は，前節と同じ．\n",
"```maple\n",
"> p0:=convert(a0,list):p1:=convert(a1,list):\n",
"> point1:=[disk(p0,0.2,color=red),disk(p1,0.2,color=blue)]:\n",
"> line1:=arrow(p0,p1,0.05,0.3,0.1):\n",
"> display(point1,line1,view=[-4..4,-4..4],gridlines=true);\n",
"```\n",
"\n",
"|{{attach_view(MapleCGplot2d4.png,)}}|\n",
"|:----|\n",
"\n",
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 平行投影図の作成 \n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"もう少し複雑な対象物として，ここでは立方体の表示を考える．まず３次元座標を打ち込む．\n",
"```maple\n",
"> restart; with(plots): with(plottools): \n",
" p:=[[0,0,0],[1,0,0],[1,1,0],[0,1,0],\n",
" [0,0,1],[1,0,1],[1,1,1],[0,1,1]]:\n",
"```\n",
"$$\n",
"$$\n",
"次にこれをpointplot3dで簡便に表示．\n",
"```maple\n",
"> points:= { seq(p[i],i=1..8) }:\n",
"> pointplot3d(points,symbol=circle,symbolsize=40,color=black);\n",
"```\n",
"\n",
"|{{attach_view(MapleCGplot3d5.png,)}}|\n",
"|:----|\n",
"\n",
"もうすこし見やすいように頂点を結んでおく．たとえば，p[1]とp[2]との間を線で結ぶには，\n",
"\n",
"```maple\n",
"> line(p[1],p[2]);\n",
"```\n",
"とする．それをseqで複数の点間に対して施す．\n",
"\n",
"```maple\n",
"> ll:=[[1,2],[2,3],[3,4],[4,1],[1,5],[2,6],[3,7],[4,8],\n",
"> [5,6],[6,7],[7,8],[8,5]]:\n",
"> lines:=[seq(line(p[ll[i][1]],p[ll[i][2]]),i=1..nops(ll))]:\n",
"> display(lines,scaling=constrained,color=black);\n",
"```\n",
"\n",
"|{{attach_view(MapleCGplot3d6.png,)}}|\n",
"|:----|\n",
"\n",
"\n",
"```maple\n",
"> l3:=display(lines,scaling=constrained,color=black):\n",
"> p3:=pointplot3d(p,symbol=circle,symbolsize=40,color=black):\n",
"> display([p3,l3],scaling=constrained,color=black);\n",
"```\n",
"\n",
"|{{attach_view(MapleCGplot3d7.png,)}}|\n",
"|:----|\n",
"\n",
"\n",
"\n",
"画像をつまんでぐるぐる回してみよ．Mapleではこんな操作は簡単にできるが，よく見ればわかるように，この3次元表示では透視図ではなく，平行投影図といわれるものを書いている．\n",
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 透視図 \n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"透視図のもっとも簡単な変換法は\n",
"```maple\n",
"> proj2d:=proc(x,z) \n",
" local x1,y1; \n",
" x1:=x[1]*z/(z-x[3]); \n",
" y1:=x[2]*z/(z-x[3]);\n",
" return [x1,y1]; \n",
" end proc:\n",
"```\n",
"$$\n",
"$$\n",
"zに視点の距離を入れて，xで座標を受け取って変換した結果を[x1,y1]として返している．この関数を前の表示と組み合わせれば透視図の描画ができる．\n",
"```maple\n",
"> z_p:=-8:\n",
" lines:=[seq(line(proj2d(p[ll[i][1]],z_p), proj2d(p[ll[i][2]],z_p)), i=1..nops(ll))]:\n",
" display(lines);\n",
"```\n",
"\n",
"|{{attach_view(MapleCGplot2d8.png,)}}|\n",
"|:----|\n",
"\n",
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Mapleの描画関数の覚書 \n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"mapleにはいくつかの描画レベルに合わせた関数が用意されている．どのような目的にどの関数（パッケージ）を使うかの選択指針として，それぞれがどのような意図で作ら\n",
"れ，それらの依存関係は以下の通り．\n",
"\n",
"- 描画の下位関数
\n",
"
\n",
"plot[structure]にあるPLOT,PLOT3Dデータ構造が一番下でCURVES,POINTS,POLYGONS,TEXTデータを元に絵を描く．\n",
"\n",
"- plottoolsパッケージ
\n",
"
\n",
"PLOTよりもう少し上位で，グラフィックスの基本形状を生成してくれる関数群．arc, arrow, circle, curve, line, point,sphereなどの関数があり，PLOT構造を吐く．表示にはplots[display]を使う．\n",
"\n",
"- plotsパッケージ
\n",
"
\n",
"簡単にグラフを書くための道具．たとえばlistplotは，listデータを簡易に表示する事を目的としている．\n",
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# 動画(Animation)\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## matplotlibでanimation\n",
"\n",
"以下は，[勾配降下法の結果を matplotlib アニメーション出力する](https://qiita.com/QUANON/items/e4a4847f338eb3251e53)\n",
"にあるのをそのまま写した．shift-returnでdirectoryにgradient_descent.gifが作られる．macではこれをfinderで選んで\n",
"spaceを押すとquick viewでanimationが見れる．"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"総ステップ数: 10792\n",
"最小値 (勾配降下法): -0.331667951428822\n",
"最小値 (f1x = 0 の解): -0.3333333333333333\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAGECAYAAADDQ9xjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd0XOWB/vHvq24Vy5ZVXOTeC65ywZgSOqYZCNh0YmoS\nsrDLLgu/lBOW7Cab7CakQIgxLVRjMMUOJWAMtnGVe+9F1SqW1evM+/tDclZxJGtkz8yd8nzO0fFo\n2n1850qPbnuvsdYiIiLSlginA4iISOBSSYiISLtUEiIi0i6VhIiItEslISIi7VJJiIhIu1QSIiLS\nLpWEiIi0SyUhIiLtivLFm6amptoBAwb44q1FRMQLNmzYUGKtTevoeT4piQEDBpCdne2LtxYRES8w\nxhzx5Hna3CQiIu1SSYiISLtUEiIi0i6VhIiItEslISIi7VJJiIhIuzw6BNYYcxioBFxAk7U2y5eh\nREQkMHTmPIlvWWtLfJZEREQCjjY3iYhIuzwtCQt8YYzZYIx5wJeBREQkcHi6uWmGtTbPGJMOfG6M\n2W2tXd76CS3l8QBAv379vBxTRESc4NGahLU2r+XfIuB9YEobz5lnrc2y1malpXU4ZpSIiASBDkvC\nGJNgjEk6eRu4HNju62AiIuI8T9YkMoCVxpgtwDrgL9baT30bS0RE2rP6QCkVdY1+mVaH+ySstQeB\ncX7IIiIiHahpaGLuK+u5aVIffjbrHJ9PT4fAiogEkc93HqO20cW1Y3v7ZXoqCRGRILJ4SwEZXWOZ\nPCDFL9NTSYiIBInymka+3lvENWN7ExFh/DJNlYSISJD4bEchjS7LdeP8s6kJVBIiIkFj8dZ8+qXE\nMzYz2W/TVEmIiASBkqp6vtlfwrXjemGMfzY1gUpCRCQofLytALeF68b18et0VRIiIkFg8ZZ8hmUk\nMrxnkl+nq5IQEQlweSdqWX+4zK87rE9SSYiIBLi/bM0H4Bo/nUDXmkpCRCTALd5SwNjMZAakJvh9\n2ioJEZEAdqikmm155Y5sagKVhIhIQFu8pXlT09VjezkyfZWEiEiAstby0ZZ8pgxIoVdyF0cyqCRE\nRALUzoIK9hdVce14ZzY1gUpCRCRgfbg5n6gIw9XnOLOpCVQSIiIByeW2fLQ5nwuHpZGSEONYDpWE\niEgAWnfoOIUVdVw/wb/DcJxKJSEiEoA+3JxHfEwkl43McDSHSkJEJMDUN7n4eFsBV4zuSZeYSEez\nqCRERALMst3FVNQ1cb2DRzWdpJIQEQkwH27OIzUxhhlDUp2OopIQEQkkFXWNLN3dfB3rqEjnf0U7\nn0BERP7m0+2FNDS5A2JTE6gkREQCyoeb8+jfI57xfbs5HQVQSYiIBIxjFXWsOlDK9eP7+PU61qej\nkhARCRCLt+RjLQGzqQlUEiIiAePDzfmc0yeZwWmJTkf5G5WEiEgAOFBcxba88oBaiwCVhIhIQPhw\nUx4RBseuQNcelYSIiMOstXywOZ/pg1NJ7xrndJy/o5IQEXHYppwTHD1eE3CbmkAlISLiuA835RET\nFcEVY3o6HeUfqCRERBzU0ORm8dYCLhuVQde4aKfj/AOVhIiIg77eW8zx6gZumujsxYXao5IQEXHQ\noo25pCbGcP7QNKejtEklISLikPKaRpbuKuK6cX2IDoARX9sSmKlERMLAkm35NLjc3Bigm5pAJSEi\n4phFG/MYlpHI6N5dnY7SLpWEiIgDDpdUs+FIGTdOzAyYEV/bopIQEXHAok15GAOzxgfupiZQSYiI\n+J21lvc35XLe4FR6JgfWMBynUkmIiPhZ9pEyco7XBvQO65NUEiIifrZoYy7xMZFcMTrwhuE4lUpC\nRMSP6hpdLNlawJWje5IQG+V0nA6pJERE/GjpriIq65q4cWKm01E8opIQEfGjRRtz6dk1jnMH93A6\nikdUEiIiflJSVc9Xe4uZNaEPkRGBe25EayoJERE/+WhzPi63DYqjmk5SSYiI+Mn7m/IY06crwzKS\nnI7iMZWEiIgf7D1Wyba8cm6YEBw7rE9SSYiI+MHC7ByiIgyzAvA61qfjcUkYYyKNMZuMMUt8GUhE\nJNQ0uty8vymfi0ek0yMx1uk4ndKZNYlHgF2+CiIiEqq+3lNMSVU9N2f1dTpKp3lUEsaYTOBqYL5v\n44iIhJ6FG3JITYzhouGBeYnS0/F0TeIZ4HHA7cMsIiIhp7SqnqW7irhhQuBeovR0OkxsjLkGKLLW\nbujgeQ8YY7KNMdnFxcVeCygiEsw+2JxPk9vy7UnBt6kJPFuTOA+4zhhzGHgbuNgY8/qpT7LWzrPW\nZllrs9LSgm+VSkTE26y1LMzOYWxmMsN7Bs+5Ea11WBLW2iettZnW2gHAHOBLa+0dPk8mIhLkduRX\nsLuwkpsnBde5Ea0F3wYyEZEgsTA7h5ioCK4bFzzDcJyqU4OZW2u/Ar7ySRIRkRBS3+Tiwy35XD4q\ng+T4aKfjnDGtSYiI+MDSXUWcqGkMynMjWlNJiIj4wMLsHHolxzFjSKrTUc6KSkJExMuOVdTx9d5i\nbpwYPNeNaI9KQkTEyxZtzMNtCdpzI1pTSYiIeJG1loUbcpg8oDsDUxOcjnPWVBIiIl608egJDhZX\nc3MIrEWASkJExKve3ZBDl+hIZo7t5XQUr1BJiIh4SU1DE4u3FDDznF4kxnbqNLSApZIQEfGSJVsL\nqKpvYs6U0NjUBCoJERGvWbA+h8FpCWT17+50FK9RSYiIeMG+Y5VsOFLGnMn9MCa4z41oTSUhIuIF\nC9bnEB1puGFi8A7m1xaVhIjIWapvcvHexlwuG5VBamKs03G8SiUhInKWPt95jLKaRmZP7ud0FK9T\nSYiInKUF63Po061L0A/m1xaVhIjIWcg5XsOKfSXcktU36Afza4tKQkTkLLyTnYMxcHNW8F6i9HRU\nEiIiZ6jJ5WZhdi4XDkujd7cuTsfxCZWEiMgZWr6vmMKKOuZMDp0zrE+lkhAROUNvrcshNTGGS0Zm\nOB3FZ1QSIiJnoKiiji93F3HTpEyiI0P3V2no/s9ERHzo3Y25uNyW2Vmhu6kJVBIiIp1mrWXB+hym\nDExhUFqi03F8SiUhItJJqw+WcqS0JqR3WJ+kkhAR6aS31uXQNS6KmeeExtXnTkclISLSCSVV9Xy6\nvYCbJmUSFx3pdByfU0mIiHTCwuxcGl2W26eG3mB+bVFJiIh4yO22vLXuKFMHpjAkPcnpOH6hkhAR\n8dDK/SUcPV7D7dP6Ox3Fb1QSIiIeemPtEVISYrhidOieYX0qlYSIiAcKy+v4YlcRN2dlEhsV+jus\nT1JJiIh4YMH6HFxuy21TwmOH9UkqCRGRDjS53Ly9/ijnD02lf48Ep+P4lUpCRKQDy/YUU1Bex+1T\nw2eH9UkqCRGRDryx9ggZXWO5ZGS601H8TiUhInIaOcdr+HpvMbMn9wvpIcHbE37/YxGRTnh7/VEM\nhMVgfm1RSYiItKOhyc2C9blcPCIjZK9h3RGVhIhIOz7feYySqvqwGaepLSoJEZF2vLH2CH26deGC\nYWlOR3GMSkJEpA0HiqtYdaCU26b2IzLCOB3HMSoJEZE2vLb6CDGREcwO0x3WJ6kkREROkV9WyzvZ\nOVw6KoPUxFin4zhKJSEi0srXe4q47tmV1DS42JJTxtd7ipyO5CiVhIhIi/KaRp5avIOK2ia6REcS\nHxPJfyzZSXlNo9PRHKOSEBFpUVRZR12jmwaXm5SEGOKio7C2+f5wpZIQEWmRnhRHeW0jEQa6xUdT\n3+TCmOb7w1WU0wFERAJFTWMTtY0uusZFU9vQXBA/uWYUyfHRTkdzjEpCRKTFm2uPYoE37ptKTFQE\n6UlxYV0QoJIQEQGgvsnFW+uOcsmIdEb3SXY6TsDQPgkREeCTbYWUVDVw17kDnI4SUFQSIiLAq6sP\nMzA1gRlDUp2OElA6LAljTJwxZp0xZosxZocx5il/BBMR8ZdtueVsOnqCO6f1JyKMx2lqiyf7JOqB\ni621VcaYaGClMeYTa+0aH2cTEfGLP68+THxMJDdNynQ6SsDpcE3CNqtq+Ta65cv6NJWIiJ+UVTfw\n0ZZ8Zk3oQ3KX8D6SqS0e7ZMwxkQaYzYDRcDn1tq1vo0lIuIf72TnUN/k5q5z+zsdJSB5VBLWWpe1\ndjyQCUwxxow59TnGmAeMMdnGmOzi4mJv5xQR8boml5s/rz7C1IEpjOjZ1ek4AalTRzdZa08Ay4Ar\n23hsnrU2y1qblZYWvldxEpHg8fnOY+SdqGXujIFORwlYnhzdlGaM6dZyuwtwGbDb18FERHztpW8O\n0TelC5eOzHA6SsDyZE2iF7DMGLMVWE/zPoklvo0lIuJbW3NPsP5wGfdMHxjWlyftSIeHwFprtwIT\n/JBFRMRvXv7mMImxUdySpcNeT0dnXItI2DlWUceSrfl8e1ImSXE67PV0VBIiEnZeX3OEJrflnukD\nnI4S8FQSIhJW6hpdvLH2KJeMyGBAaoLTcQKeSkJEwsqHm/M4Xt3A3BkDnI4SFFQSIhI2rLW8/M1h\nRvRM4txBPZyOExRUEiISNlYfKGV3YSVzZwzEGB326gmVhIiEjZe+OUSPhBiuG9fb6ShBQyUhImHh\nUEk1S3cXcfu0/sRFRzodJ2ioJEQkLLy66jBREYY7pvVzOkpQUUmISMgrr23knewcrh3Xm/SkOKfj\nBBWVhIiEvAXrj1LT4GLueRrttbNUEiIS0hqa3Ly08jDnDurBmD7JTscJOioJEQlpS7bmU1hRxwMX\nDHI6SlBSSYhIyLLWMm/5QYamJ3LRcF0M7UyoJEQkZK3cX8Luwkruv2CQTp47QyoJEQlZ85YfJD0p\nluvH6+S5M6WSEJGQtDO/ghX7SrjnvAHERunkuTOlkhCRkPTCioPEx0Ry+5T+TkcJaioJEQk5+Sdq\nWbwlnzmT+5EcryvPnQ2VhIiEnFdWHcaCrhnhBSoJEQkpFXWNvLn2KFef04vM7vFOxwl6KgkRCSlv\nrztKVX2TTp7zEpWEiISMk0NwTB+sITi8RSUhIiHj5BAc92stwmtUEiISEk4OwTEsI5GLhmkIDm9R\nSYhISPhqTzG7Cyt54ILBGoLDi1QSIhISnvtqP326ddEQHF6mkhCRoLfu0HHWHy7j/vMHEh2pX2ve\npLkpIkHvua/20yMhhtmTdf1qb1NJiEhQ25Ffzld7ipk7YyBdYjSQn7epJEQkqD331QESY6O4Y5oG\n8vMFlYSIBK1DJdV8sq2AO6b1J7mLBvLzBZWEiAStP319gOjICO6dMdDpKCFLJSEiQamgvJb3NuZy\nS1Zf0pJinY4TslQSIhKU5q84hNuigfx8TCUhIkGnrLqBN9ce5fpxvembouHAfUklISJB5+VVh6lt\ndPHQRYOdjhLyVBIiElSq6pt4ddVhLhuVwbCMJKfjhDyVhIgElbfWHqW8tpHvaS3CL1QSIhI0ahtc\n/Gn5QaYP7sGEft2djhMWVBIiEjTeXHeUkqp6HrlkqNNRwoZKQkSCQl2ji+e/PsC0QSlMHdTD6Thh\nQyUhIkHh7XVHKa6s55FLhjkdJayoJEQk4NU1uvjj1weYMiCFaYNSnI4TVlQSIhLwFmbncKyinkcu\nHapLk/qZSkJEAlp9k4vnvjrApP7dmT5Y+yL8TSUhIgHt3Q25FJTX8cglWotwgkpCRAJWQ5Ob55Yd\nYEK/bpw/NNXpOGFJJSEiAWvRxlzyTtTyT1qLcIxKQkQCUqPLzR+W7WdcZjIXDUtzOk7YUkmISEB6\nf1MeuWVai3CaSkJEAk6Ty82zy/Yzpk9XLh6R7nScsNZhSRhj+hpjlhljdhpjdhhjHvFHMBEJX4s2\n5XGktIZ/ulhrEU6L8uA5TcBj1tqNxpgkYIMx5nNr7U4fZxORMFTf5OK3X+xjbGYyl43KcDpO2Otw\nTcJaW2Ct3dhyuxLYBfTxdTARCU/vrM8h70Qtj10+XGsRAaBT+ySMMQOACcBaX4QRkfBW1+ji91/u\nZ8qAFC7QeREBweOSMMYkAu8Bj1prK9p4/AFjTLYxJru4uNibGUUkTLy2+ghFlfU8dvkwrUUECI9K\nwhgTTXNBvGGtXdTWc6y186y1WdbarLQ0HdMsIp1TVd/EH78+wPlDU3W9iADiydFNBngR2GWt/bXv\nI4lIOHp55SGOVzfw2OXDnY4irXiyJnEecCdwsTFmc8vXTB/nEpEwUl7TyLwVB7l0ZAbj+3ZzOo60\n0uEhsNbalYA2DoqIz7yw4iCVdU38y2W66lyg0RnXIuKokqp6XvrmENeM7cWo3l2djiOnUEmIiKOe\n/+oAdY0uHr1UaxGBSCUhIo4pLK/jtTVHuGFCJkPSE52OI21QSYiIY/6wbB8ut+XRS4c6HUXaoZIQ\nEUccKqnm7XU5zJ7cl74p8U7HkXaoJETEEf/z1z1ER0bwiNYiAppKQkT8bkvOCf6ytYD7zx9IelKc\n03HkNFQSIuJX1lp+/skueiTE8MCFg52OIx1QSYiIX321p5g1B4/zT5cMJTHWk0vaiJNUEiLiNy63\n5b8/3U3/HvHcOqWf03HEAyoJEfGb9zflsbuwkn+7YjgxUfr1Ewz0KYmIX9Q1uvj1X/cwNjOZmWN6\nOR1HPKSSEBG/eHXVYfLL63jiqhFERGjM0GChkhARnztR08Czy/Zz0fA0pg/WZUmDiUpCRHzuj18d\noLK+iX+/coTTUaSTVBIi4lN5J2p5edVhbpjQh5G9NBR4sFFJiIhP/c9newB0WdIgpZIQEZ/ZnHOC\n9zflcd+MgfTp1sXpOHIGVBIi4hPWWv5j8Q5SE2P53reGOB1HzpBKQkR8YvHWAjYePcHjVwzX8BtB\nTCUhIl5X1+jiFx/vYnTvrtw0KdPpOHIWVBIi4nXzVxwkv7yOH18zikidOBfUVBIi4lXHKup47qsD\nXDm6J9MG9XA6jpwllYSIeNWvPttDk8vy5EydOBcKVBIi4jXbcst5d0Mu35kxgP49EpyOI16gkhAR\nr7DW8vSSnfRIiOFhHfIaMlQSIuIVn2wvZN3h4zx2+XCS4qKdjiNeopIQkbNW1+jivz7exYieScye\n3NfpOOJFKgkROWvzlh8kt6yWn+iQ15CjkhCRs5JzvIZnl+3n6rG9mD5E14oINSoJETkrP/vLTiKM\n4YczRzodRXxAJSEiZ+zrvcV8tuMYP7hkCL01ymtIUkmIyBmpb3Lx0492MDA1gXtnDHQ6jviIhmYU\nkTPy4spDHCqp5tW5U4iNinQ6jviI1iREpNMKymv5/dL9XD4qgwuHpTkdR3xIJSEinfaff9mF21p+\nfM0op6OIj6kkRKRTVuwrZsnWAr530RD6psQ7HUd8TCUhIh6ra3Tx4w+2MzA1gQcvHOR0HPED7bgW\nEY89t2w/h0treOO+qcRFa2d1ONCahIh4ZH9RFX/8+gCzxvfmPJ1ZHTZUEiLSIWstP/pgG12iI/nh\n1dpZHU5UEiLSoUUb81hz8DhPXDWStKRYp+OIH6kkROS0yqob+M+PdzGxXzfmaBjwsKOSEJHT+sUn\nuymvbeQ/bziHCA0DHnZUEiLSrjUHS1mQncN9MwYysldXp+OIA1QSItKmukYXT7y3lX4p8Tx66TCn\n44hDdJ6EiLTpN1/s5XBpDW/eN5UuMTonIlxpTUJE/sG23HLmrzjEnMl9dbW5MKeSEJG/0+hy8/h7\nW+mREMOTutpc2NPmJhH5O/OWH2RXQQV/unMSyV2inY4jDtOahIj8zf6iKn67dB8zz+nJFaN7Oh1H\nAoBKQkQAcLstTy7aSpfoSH563Win40iA6LAkjDEvGWOKjDHb/RFIRJzx0jeHWH+4jB9fM4r0pDin\n40iA8GRN4hXgSh/nEBEHHSiu4lef7eGSEencNLGP03EkgHRYEtba5cBxP2QREQc0udw89s4W4qIj\n+fmN52CMht6Q/xNw+yT+uqOQV1cddjqGSNiYt+Igm3NO8PSsMaR31WYm+XteKwljzAPGmGxjTHZx\ncfEZv89HW/J5eslOtueVeyuaiLRjT2Elz3zefDTTtWN7OR1HApDXSsJaO89am2WtzUpLSzvj9/nZ\nrDH0SIzhnxdspq7R5a14InKKRpebf3lnM0lxUTx9/RhtZpI2Bdzmpm7xMfzy2+PYV9S8I01EfOPZ\nZfvZkV/Bf95wDj0SdSEhaZsnh8C+BawGhhtjco0x9/o61IXD0rjr3P68uPIQq/aX+HpyImFnS84J\n/vDlfmaN782VY3TSnLTPk6ObbrXW9rLWRltrM621L/oj2JNXjWRQagL/unAL5bWN/pikSFiorm/i\n0QWbSU+K5anrxzgdRwJcwG1uOqlLTCS/mT2eY5X1PPXRDqfjiISMn/1lJ4dLq/n17PEam0k6FLAl\nATCubzce/tYQFm3K4+NtBU7HEQl6n+0o5K11OTx4wWCmDerhdBwJAgFdEgAPXzyEcZnJPLloG/kn\nap2OIxK0iirqeOK9rYzp05V/uUxXmhPPBHxJREdG8MycCTS53Dy6YDMut3U6kkjQcbstjy3cQm2j\ni2dmTyAmKuB/9CVABMWSMjA1gadnjWHdoeM8u2y/03FEgs4rqw6zYl8JP7p6FEPSE52OI0EkKEoC\n4MaJmdwwoQ/PfLGX7MMaSkrEUzvyy/nFp7u5ZEQ6t0/t53QcCTJBUxIA/3H9aPqmxPPI25spr9Fh\nsSIdqapv4uE3N9E9Pppf3TxOZ1VLpwVVSSTFRfO7ORM4VlHHE4u2Yq32T4i0x1rLj97fxpHSan43\nZwIpCTFOR5IgFFQlAc2Hxf7bFcP5ZHvzoXwi0raFG3L5YHM+/3zpMKbqcFc5Q0FXEgD3nz+I84em\n8tTiHezI12ixIqfad6ySn3y4nfOG9OB73xridBwJYkFZEhERht/MHk+3+Gi+98ZGKuq0f0LkpNoG\nF99/cyOJsVH8ZvZ4IiO0H0LOXFCWBEBqYizP3jaR3LJaHl+o/RMiJ/30ox3sK6riN7PH61rVctaC\ntiQAsgak8MSVI/h0RyEvrjzkdBwRx72zPocF2Tl876LBnD/0zK/rInJSUJcEwH3nD+TyURn84pPd\nbDii8yckfG3PK+dHH25nxpBU/uWy4U7HkRAR9CVhjOFXN4+jd7cuPPzmJkqr6p2OJOJ3ZdUNPPT6\nBlITYvjtHO2HEO8J+pIASO4SzXO3T6S0uoFH3t5Mk8vtdCQRv3G5LY8s2ExRRT3P3TFJV5kTrwqJ\nkgAY0yeZn10/hpX7S/ilLnsqYeS3S/exfG8xP71uNOP7dnM6joSYKKcDeNMtk/uyPb+cecsPMrp3\nV64f38fpSCI+tXTXMX63dB83T8rk1il9nY4jIShk1iRO+vE1o5gyMIXH393K9jydaCeh60BxFY8u\n2Mzo3l15etYYjcskPhFyJREdGcFzt08kJSGGB1/boB3ZEpLKaxq5/9VsoiMjeP6OScRFRzodSUJU\nyJUENJ9oN+/OLEqq6vn+mxtp1I5sCSFNLjcPv7WRnLIanr9jEn1T4p2OJCEsJEsC4JzMZH5+4zms\nOXicny3Z6XQcEa/5r493s2JfCU9fP4YpA1OcjiMhLqR2XJ/qxomZ7MyvYP7KQwxKS+Tu6QOcjiRy\nVt5Zn8NL3xzinukDmDNFFxAS3wvpkgB4cuZIDpdW89TiHfTrEc+3hqc7HUnkjGQfPs4PP9jG+UNT\n+dHVI52OI2EiZDc3nRQZYfjtnAmM6NmVH7y5iT2FlU5HEum0o6U1PPjaBjK7x/OHWycSFRnyP7oS\nIMJiSUuIjeLFe7KIj4lk7ivrKaqsczqSiMfKqhu45+V1uKzlxbuzSI6PdjqShJGwKAmAXsldePHu\nyRyvbuCBP2+grtHldCSRDtU1unjgtWxyT9Tywl1ZDEpLdDqShJmwKQloPuLpmTnj2ZJ7gn9esBmX\nW9egkMDldlseW7iF9YfL+N+bxzF5gI5kEv8Lq5IAuGJ0T344cySfbC/kqcU7dLEiCVj//dlu/rK1\ngCeuGsG143o7HUfCVMgf3dSW+84fRFFlPfOWHySjaxzf1zWAJcC8vuYIf/r6ILdP7ceDFwxyOo6E\nsbAsCYAnrhxBUUUdv/psD2mJsdwyWYOjSWD4eFsBP/lwOxePSOep60ZrTCZxVNiWRESE4ZffHkdp\ndQNPvr+NHokxXDIyw+lYEuZW7Cvmkbc3MaFfd569TYe6ivPCegmMiYrgj3dMYlSvrnz/zY1sOFLm\ndCQJY5uOlvHgaxsYnJbIS3dPpkuMBu0T54V1SQAkxkbx8ncmk9E1ju+8vI4d+RpeXPxv77FKvvPK\nelITY/nz3Ck6F0ICRtiXBDSPGvv6vVNJjI3izhfXse+YzsoW/8k5XsOdL64lJjKC1++dSnrXOKcj\nifyNSqJF35R43rh/GpERhtvnr+VwSbXTkSQMFJbXcceLa6ltcPHne6fQr4eG/ZbAopJoZWBqAm/c\nN5VGl5vb568l70St05EkhB2rqOPWF9ZQWtXAK3OnMKJnV6cjifwDlcQphmUk8dq9U6moa+T2F9ZQ\nVKFxnsT7iiqbC6Kooo5X505mYr/uTkcSaZNKog1j+iTzynemUFRZz5wX1nBMRSFeVFxZz20vrKWw\nvI5X5k5hUn8NtyGBSyXRjkn9u/Pq3CkcK69j9p9Wk69NT+IFpVX13D5/DbllNbx0z2SNxyQBTyVx\nGpMHpPDne6dSWtXA7HmryTle43QkCWJFlXXc9sJajpTW8NLdk5k2qIfTkUQ6pJLowKT+3Xn9vqmU\n1zQyZ94ajpaqKKTzcstquOX51Rw93rwGMX1IqtORRDyikvDAuL7dePP+aVQ3NHHLn1ZzsLjK6UgS\nRA6VVHPL86sprW7g9fumcJ4KQoKISsJDY/ok89b902h0ubn5+dVsz9OZ2dKx3YUV3Pz8auqa3Lx1\n/zTtpJago5LohJG9uvLOQ+cSFx3JnHlrWLW/xOlIEsC25Jxgzrw1REbAOw9OY0yfZKcjiXSaSqKT\nBqcl8t53p9OnWxfueXk9H28rcDqSBKBle4q49YU1JMVFsfDB6QxJT3I6ksgZUUmcgZ7Jcbzz4LmM\nzUzm+29u5LU1R5yOJAFkwfqj3PdqNgNTE3jvoekaakOCmkriDCXHR/PavVO5eHg6P/5gO//z2R7c\numZ2WLMD9EwFAAAN+ElEQVTW8uvP9/Lv721jxpBUFjx4rgbrk6CnkjgLXWIief7OSczO6ssflu3n\nB29voq7R5XQscUCjy82/vbuV3y3dxy1Zmcy/O4vE2LC9ppeEEC3FZyk6MoJf3HQOg9IS+MWnu8kt\nq+WFuyaRnqS/IMPFiZoGHn5zEyv3l/DopUN55JKhuuSohAytSXiBMYYHLxzM83dMYm9hJbP+8A07\n8yucjiV+sPdYJdc/+w3rDh3nV98ey6OXDlNBSEhRSXjRFaN7svChc3FZy7efX8VnOwqdjiQ+9Ncd\nhdzw7DdU17t464Fp3JzV1+lIIl6nkvCyMX2S+fD7MxiSnsiDr23gvz/dTZPL7XQs8SJrLb9fuo8H\nXtvA4PREFv/gPCb111DfEpo8KgljzJXGmD3GmP3GmCd8HSrYnTxE9tYpffnjVwe466V1lFTVOx1L\nvKC8tpHvvr6R//18L7PG9+adB8+lV3IXp2OJ+EyHJWGMiQSeBa4CRgG3GmNG+TpYsIuLjuTnN47l\nl98eS/aRMq79/Uo2HS1zOpacha25J7jm9yv4fNcxfjhzJL+ZPZ646EinY4n4lCdrElOA/dbag9ba\nBuBt4Hrfxgodt2T1ZdF3pxMZYbjlT6t5aeUhrNX5FMHEWsvL3xzipj+uwuWyvPPgNO6/YJB2UEtY\n8KQk+gA5rb7Pbbnv7xhjHjDGZBtjsouLi72VLySM6ZPMkh/M4MJhafzHkp3c8/J6iip1tbtgUF7b\nyEOvb+CpxTu5YGgaf/mn8zVIn4QVr+24ttbOs9ZmWWuz0tLSvPW2IaNbfAwv3JXF07PGsOZgKVc9\ns4Klu445HUtO45v9JVz1zHKW7irihzNHMv/uLLonxDgdS8SvPCmJPKD1sX2ZLfdJJxljuHNaf5b8\nYAbpXeO499VsfvzBdmobdJZ2IKlpaOInH27n9vlriYuJ5N3vTtfmJQlbnpxxvR4YaowZSHM5zAFu\n82mqEDc0I4kPvj+dX326h/krD7F8XzE/v/Ecpg/WxWictuFIGY+9s5nDpTXMPW8gj185XDunJax1\nuCZhrW0CHgY+A3YB71hrd/g6WKiLjYrkR9eM4s37pwJw2wtreXLRVsprGx1OFp5qGpr4r493cfPz\nq2h0Wd66fxo/uXaUCkLCnvHFkTZZWVk2Ozvb6+8bqmobXPzmi73MX3GQ1MRYnp41hitG93Q6Vtj4\n645Cnlq8k7wTtdw6pS//b+ZIkuKinY4l4lPGmA3W2qwOn6eSCBxbc0/w+Ltb2V1YyaUj0/nR1aMY\nkJrgdKyQlVtWw08/2skXu44xomcSP5s1hqwBOnJJwoNKIkg1uty8uPIQv1+6jwaXm7kzBvKDi4dq\n2Gkvqm1w8dI3h/jDl/sB+OfLhvKd8wYSHalRaiR8qCSCXFFFHb/8bA/vbsglLSmWx68Yzk0TM4mI\n0BE2Z8rltry3MZdf/3UvhRV1XDE6g59cO5o+3TSshoQflUSI2JJzgp8u3sGmoycY0TOJRy8dxhWj\nM3Q4ZidYa/l6bzG/+GQ3uwsrGd+3G/9v5kimDNSmJQlfKokQ4nZblmwr4JnP93KwpJoxfbry2GXD\nuWh4msriNKy1rNxfwu+/3M+6Q8fplxLP41cO5+pzemm+SdhTSYSgJpebDzbn89ule8k5XsuEft34\n/kVDuHhEujZDtWKtZemuIn6/bD9bck6Q0TWWhy4czG1T+xEbpUNaRUAlEdIaXW4WZufyhy/3kV9e\nx6DUBObOGMhNEzPpEhO+vwTrGl18vK2AF1YcYldBBX1TuvDdC4dw06Q+KgeRU6gkwkCjy80n2wuZ\nv+IgW3PL6RYfzR1T+zNnSl8yu8c7Hc9vcstqeGPtURasz+F4dQND0hP53kWDuW5cb6J0xJJIm1QS\nYcRaS/aRMuavOMhfdzYPGnje4FRuzsrkitE9Q/Ks4fomF8v3lrBgfQ5f7m7+P186MoO7zh3AeUN6\naJ+DSAdUEmEq53gN723M5d0NueSW1ZIUF8V143pz9dheTBmQEtR/WbvclrWHSvlocz4fbyugoq6J\n1MQY5kzux21T+9Fbh7KKeEwlEebcbsuag6Us3JDLJ9sLqGt00y0+motHpHP5qJ5cMCyV+JjAP0Gv\nur6JVQdK+WpPEV/sOsaxinoSYiK5YnRPrhvfm/OGpOokOJEzoJKQv6lpaGL53hL+uqOQpbuLKK9t\nJDYqgon9ujNtUA+mDkphfN9uAbFZqqHJzc6CCtYeLOXrvcWsP3ycRpclISaSGUNTuXZcby4ZkRHW\nO+hFvEElIW1qdLlZf+g4X+wqYs3BUnYVVmAtxERFML5vN8b0TmZEryRG9erKkPREnxZHfZOLI6U1\n7CmsZNPRE2zOKWN7fgUNTW4AhmckcdHwNC4cnkZW/xRiorTGIOItKgnxSHlNI+sOH2ftwVLWHylj\nT2EFdY3Nv6QjIwwDesST2T2e3t3i6JXchd7dutCzaxxJcVEkxEYSHxNFQmwUcdERuN3Q6HbT5LI0\nudzUN7kprW7geHU9JVUNlFY1cKyijkMl1RwsqSKvrBZ3y+IXGxXB2MxkJvTrzvi+3ZjYrzs9k+Mc\nnDMioU0lIWfE5bYcKa1md2Eluwoq2HuskvwTdRSU11JS1XDW7x8fE8nA1AQGpiYwKC2RwWkJDE5L\nZHjPJO1bEPEjT0si8Pdcil9FRhgGpSUyKC2Rmef0+rvH6hpdFJbXcayijuqGJqrqXdTUN1Hd4KKu\n0UWEMURHGqIiDFGREcRERdAjIYaUhBhSE2PpkRgTFDvLReT/6CdWPBYXHcmA1ARd40IkjGj9XkRE\n2qWSEBGRdqkkRESkXSoJERFpl0pCRETapZIQEZF2qSRERKRdKgkREWmXSkJERNqlkhARkXapJERE\npF0qCRERaZdKQkRE2qWSEBGRdvnkokPGmGLgyFm8RSpQ4qU4vhIMGSE4cgZDRgiOnMGQEYIjZ6hn\n7G+tTevoST4pibNljMn25IpJTgqGjBAcOYMhIwRHzmDICMGRUxmbaXOTiIi0SyUhIiLtCtSSmOd0\nAA8EQ0YIjpzBkBGCI2cwZITgyKmMBOg+CRERCQyBuiYhIiIBwLGSMMbcbIzZYYxxG2Pa3TtvjLnS\nGLPHGLPfGPNEq/tTjDGfG2P2tfzb3QcZO5yGMWa4MWZzq68KY8yjLY/91BiT1+qxmd7O6GnOlucd\nNsZsa8mS3dnX+zqjMaavMWaZMWZny7LxSKvHfDYv21vGWj1ujDG/a3l8qzFmoqev9SYPct7ekm+b\nMWaVMWZcq8fa/OwdyHiRMaa81ef4E09f68eM/9Yq33ZjjMsYk9LymL/m40vGmCJjzPZ2HvffMmmt\ndeQLGAkMB74Cstp5TiRwABgExABbgFEtj/0SeKLl9hPAf/sgY6em0ZK3kObjjwF+CvyrH+alRzmB\nw0Dq2f4/fZUR6AVMbLmdBOxt9Xn7ZF6ebhlr9ZyZwCeAAaYBaz19rZ9zTge6t9y+6mTO0332DmS8\nCFhyJq/1V8ZTnn8t8KU/52PLdC4AJgLb23ncb8ukY2sS1tpd1to9HTxtCrDfWnvQWtsAvA1c3/LY\n9cCrLbdfBWb5IGZnp3EJcMBaezYnEp6Js50XATEvrbUF1tqNLbcrgV1AHx9kae10y9hJ1wN/ts3W\nAN2MMb08fK3fclprV1lry1q+XQNk+ijLGWf00Wt9mfFW4C0f5Dgta+1y4PhpnuK3ZTLQ90n0AXJa\nfZ/L//3SyLDWFrTcLgQyfDD9zk5jDv+4QP2gZXXwJV9sxmnhaU4LfGGM2WCMeeAMXu+PjAAYYwYA\nE4C1re72xbw83TLW0XM8ea23dHZa99L8l+ZJ7X323uRpxuktn+MnxpjRnXytvzJijIkHrgTea3W3\nP+ajJ/y2TEadzYs7Yoz5AujZxkM/tNZ+6K3pWGutMeaMDtM6XcbOTMMYEwNcBzzZ6u4/Ak/TvGA9\nDfwvMNfBnDOstXnGmHTgc2PM7pa/WDx9vT8yYoxJpPkH81FrbUXL3V6bl6HOGPMtmktiRqu7O/zs\n/WQj0M9aW9WyX+kDYKgDOTxxLfCNtbb1X/SBMh/9xqclYa299CzfIg/o2+r7zJb7AI4ZY3pZawta\nVrOKvJ3RGNOZaVwFbLTWHmv13n+7bYx5AVhyJhm9ldNam9fyb5Ex5n2aV02XE0Dz0hgTTXNBvGGt\nXdTqvb02L09xumWso+dEe/Bab/EkJ8aYscB84CprbenJ+0/z2fs1Y6vSx1r7sTHmOWNMqiev9VfG\nVv5hy4Cf5qMn/LZMBvrmpvXAUGPMwJa/1OcAH7U89hFwd8vtuwGvrZm00plp/MO2y5ZfhifdALR5\npIIXdJjTGJNgjEk6eRu4vFWegJiXxhgDvAjsstb++pTHfDUvT7eMtc5+V8sRJdOA8pZNZ5681ls6\nnJYxph+wCLjTWru31f2n++z9nbFny+eMMWYKzb+DSj15rb8ytmRLBi6k1XLqx/noCf8tk77YM+/J\nF80/6LlAPXAM+Kzl/t7Ax62eN5Pmo1wO0LyZ6uT9PYClwD7gCyDFBxnbnEYbGRNoXtCTT3n9a8A2\nYGvLB9XLR/Oyw5w0H+2wpeVrRyDOS5o3j9iW+bW55Wumr+dlW8sY8BDwUMttAzzb8vg2Wh2N197y\n6aPPuaOc84GyVvMuu6PP3oGMD7dk2ELzzvXp/p6XHWVs+f4e4O1TXufP+fgWUAA00vx78l6nlkmd\ncS0iIu0K9M1NIiLiIJWEiIi0SyUhIiLtUkmIiEi7VBIiItIulYSIiLRLJSEiIu1SSYiISLv+P+7p\n8CgrMKXVAAAAAElFTkSuQmCC\n",
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"%matplotlib inline\n",
"\n",
"from sympy import Derivative, Symbol, sympify, solve\n",
"from numpy import arange\n",
"import matplotlib.pyplot as plt\n",
"import matplotlib.animation as ani\n",
"\n",
"\n",
"def gradient_descent(x0, f1x, x, epsilon=1e-6, step_size=1e-4):\n",
" # f1x = 0 の解を持つか調べる。\n",
" if not solve(f1x):\n",
" return\n",
"\n",
" x_old = x0\n",
" x_new = x_old - step_size * f1x.subs({x: x_old}).evalf()\n",
" X_traversed = []\n",
"\n",
" while abs(x_old - x_new) > epsilon:\n",
" X_traversed.append(x_new)\n",
" x_old = x_new\n",
" x_new = x_old - step_size * f1x.subs({x: x_old}).evalf()\n",
"\n",
" return x_new, X_traversed\n",
"\n",
"\n",
"def draw_graph(f, x):\n",
" X = arange(-1, 1, 0.01)\n",
" Y = [f.subs({x: x_val}) for x_val in X]\n",
"\n",
" plt.plot(X, Y)\n",
"\n",
"\n",
"def draw_frame(i, x, X, Y):\n",
" plt.clf()\n",
" draw_graph(f, x)\n",
" plt.scatter(X[i], Y[i], s=20, alpha=0.8)\n",
"\n",
"\n",
"if __name__ == '__main__':\n",
" x = Symbol('x')\n",
" f = 3 * x ** 2 + 2 * x\n",
" var0 = 0.75 # 勾配降下法の初期値\n",
"\n",
" d = Derivative(f, x).doit()\n",
"\n",
" # gradient_descent() は、勾配降下法で求めた最小値と各ステップでの x の値を返す。\n",
" var_min, X_traversed = gradient_descent(var0, d, x)\n",
"\n",
" print('総ステップ数: {0}'.format(len(X_traversed)))\n",
" print('最小値 (勾配降下法): {0}'.format(var_min))\n",
" print('最小値 (f1x = 0 の解): {0}'.format(float(solve(d)[0])))\n",
"\n",
" X = X_traversed[::100] # (1)\n",
" Y = [f.subs({x: x_val}) for x_val in X]\n",
"\n",
" fig = plt.figure(figsize=(6.5, 6.5))\n",
"\n",
" anim = ani.FuncAnimation(fig, draw_frame, fargs=(x, X, Y), frames=len(X)) # (2)\n",
" anim.save('gradient_descent.gif', writer='imagemagick', fps=10) # (3)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## plotの動画\n",
"\n",
"もう一つ．[ほげおのクラフト](https://hogeocraft.blogspot.jp/2016/12/jupyter-notebook-python3.html)の例．ffmpegでエラーが出た．Mac版をinstallして．．．\n",
"> brew install ffmpeg\n",
"\n",
"再度動かしてみると動いた！！？"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"%matplotlib inline\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"\n",
"from matplotlib import animation\n",
"from IPython.display import HTML"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"dn = np.random.choice([-1,1], size=100)\n",
"swalk = np.cumsum(dn)\n",
"\n",
"# First set up the figure, the axis, and the plot element we want to animate\n",
"fig, ax = plt.subplots()\n",
"\n",
"ax.set_xlim(( 0, 100))\n",
"ax.set_ylim((-10, 10))\n",
"\n",
"line, = ax.plot([], [], lw=2)\n",
"plt.close()\n",
"\n",
"# initialization function: plot the background of each frame\n",
"def init():\n",
" line.set_data([], [])\n",
" return (line,)\n",
"# animation function. This is called sequentially\n",
"def animate(i):\n",
" x = np.arange(i)\n",
" y = swalk[:i]\n",
" line.set_data(x, y)\n",
" return (line,)\n",
"# equivalent to rcParams['animation.html'] = 'html5'\n",
"plt.rcParams['animation.html'] = 'html5'\n",
"\n",
"# call the animator. blit=True means only re-draw the parts that have changed.\n",
"anim = animation.FuncAnimation(fig, animate, init_func=init,\n",
" frames=100, interval=50, blit=True)"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
""
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"anim"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Runge-Kutta4次公式\n",
"2次と同様の考え方で，$h^4$の精度を持つRunge-Kutta４次公式を作ることができる．\n",
"\n",
"微分方程式\n",
"$$\n",
"\\frac{dy}{dx} = f(x,y), \\, {\\textrm where}y(x_0)=y_0\n",
"$$\n",
"\n",
"の数値解は，刻み幅を$h$，$x_n=x_0+nh$として，次の漸化式\n",
"$$\n",
"y_{n+1} = y_n +k (n=0,1,2,\\cdots)\n",
"$$\n",
"で与えられる．ここに，$k$は次で定める．\n",
"$$\n",
"\\begin{aligned}\n",
"k_1 & = hf(x_n,y_n), \\\\\n",
"k_2 & = hf(x_n+\\frac{h}{2}, y_n+\\frac{k_1}{2}), \\\\\n",
"k_3 & = hf(x_n+\\frac{h}{2}, y_n+\\frac{k_2}{2}), \\\\\n",
"k_4 & = hf(x_n+h, y_n+k_3), \\\\\n",
"k & = \\frac{1}{6}(k_1+2k_2+2k_3+k_4)\n",
"\\end{aligned}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 連立方程式にRunge-Kutta4次公式を\n",
"連立微分方程式\n",
"$$\n",
"\\begin{aligned}\n",
"\\frac{dy}{dx} &= f(x,y,z), &\\, \\, where \\, y(x_0)=y_0 \\\\\n",
"\\frac{dz}{dx} &= g(x,y,z), &\\, \\, where \\, z(x_0)=z_0 \\\\\n",
"\\end{aligned}\n",
"$$\n",
"\n",
"の数値解は，刻み幅を$h$，$x_n=x_0+nh$として，次の漸化式\n",
"$$\n",
"\\begin{aligned}\n",
"y_{n+1} & = y_n +k &\\, (n=0,1,2,\\cdots) \\\\\n",
"z_{n+1} & = z_n +l &\\, (n=0,1,2,\\cdots) \\\\\n",
"\\end{aligned}\n",
"$$\n",
"で与えられる．ここに，$k,l$は次で定める．\n",
"$$\n",
"\\begin{aligned}\n",
"k_1 &= hf(x_n,y_n,z_n), \\,\n",
"&l_1 &= hg(x_n,y_n,z_n), \\\\\n",
"k_2 &= hf(x_n+\\frac{h}{2}, y_n+\\frac{k_1}{2}, z_n+\\frac{l_1}{2}), \\,\n",
"&l_2 &= hg(x_n+\\frac{h}{2}, y_n+\\frac{k_1}{2}, z_n+\\frac{l_1}{2}), \\\\\n",
"k_3 &= hf(x_n+\\frac{h}{2}, y_n+\\frac{k_2}{2}, z_n+\\frac{l_2}{2}), \\,\n",
"&l_3 &= hg(x_n+\\frac{h}{2}, y_n+\\frac{k_2}{2}, z_n+\\frac{l_2}{2}), \\\\\n",
"k_4 &= hf(x_n+h, y_n+k_3, z_n+l_3), \\,\n",
"&l_4 &= hg(x_n+h, y_n+k_3, z_n+l_3), \\\\\n",
"k &= \\frac{1}{6}(k_1+2k_2+2k_3+k_4), \\,\n",
"&l &= \\frac{1}{6}(l_1+2l_2+2l_3+l_4)\n",
"\\end{aligned}\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"def runge_kutta4(x0,y0,z0,h):\n",
" k1= h * ff(x0, y0, z0);\n",
" l1= h * gg(x0, y0, z0);\n",
" k2= h * ff(x0 + h/2, y0 + k1/2, z0 + l1/2)\n",
" l2= h * gg(x0 + h/2, y0 + k1/2, z0 + l1/2)\n",
" k3= h * ff(x0 + h/2, y0 + k2/2, z0 + l2/2)\n",
" l3= h * gg(x0 + h/2, y0 + k2/2, z0 + l2/2)\n",
" k4= h * ff(x0 + h, y0 + k3, z0 + l3)\n",
" l4= h * gg(x0 + h, y0 + k3, z0 + l3)\n",
" kk= 1.0/6*(k1 + 2*k2 + 2*k3 + k4)\n",
" ll= 1.0/6*(l1 + 2*l2 + 2*l3 + l4)\n",
" return [kk,ll]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Runge-Kuttaの４次公式をそのままcodingすると上のようになります．これを先ほどのバネ運動の問題に当てはめてみましょう．\n",
"\n",
"先ほど導出した運動方程式の漸化式\n",
"$$\n",
"\\begin{aligned}\n",
"v_{i+1} & = v_i - \\frac{k}{m} x_i\\, dt \\\\\n",
"x_{i+1} & = x_i + v_i\\, dt\n",
"\\end{aligned}\n",
"$$\n",
"\n",
"とRunge-Kutta4次公式を示した連立微分方程式\n",
"$$\n",
"\\begin{aligned}\n",
"\\frac{dy}{dx} &= f(x,y,z), \\, where \\, y(x_0)=y_0 \\\\\n",
"\\frac{dz}{dx} &= g(x,y,z), \\, where \\, z(x_0)=z_0 \\\\\n",
"\\end{aligned}\n",
"$$\n",
"とを比べて，変数の表記の違いと関数$f,g$を具体的に書き下します．\n",
"\n",
"| 4次の公式 | 運動方程式 |\n",
"|:-----:|:-----:|\n",
"|x | t |\n",
"|y | x |\n",
"|z | v |\n",
"|f(x,y,z) | f(t,x,v) = v |\n",
"|g(x,y,z) | g(t,x,v) = -k x |\n",
"\n",
"この変数の書き換えを吸収する中間関数としてEuler3を書き換えます．RungeKutta4の仮引数を上の表に従って置き換えて，数値を渡しています．また，関数$f,g$も定義しておきます．"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"def ff(t,x,v):\n",
" return v\n",
"\n",
"def gg(t,x,v):\n",
" return -k*x\n",
"\n",
"def ode(x0, v0):\n",
" kk,ll = runge_kutta4(0, x0, v0, dt)\n",
" x1 = x0 + kk\n",
" v1 = v0 + ll\n",
" return [x1, v1]\n",
"\n",
"t, dt, k=0.0, 1, 0.01\n",
"tt,xx,vv=[0.0],[0.0],[0.01]\n",
"for i in range(0,500):\n",
" t += dt\n",
" x, v = ode(xx[-1],vv[-1])\n",
" tt.append(t)\n",
" xx.append(x)\n",
" vv.append(v)\n",
"\n",
"#my_plot(xx, vv, tt)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"これを前と同様に走らせると\n",
"発散も収束もすることなく，定常的に振動を繰り返していることが見て取れます．"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAZIAAAD8CAYAAABdCyJkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFHRJREFUeJzt3W+MXNd93vHvU9J0C1m1pIalaZKCKGBjgw4CRh5QAmIb\naSMlJBuYcl8oFNCadQXQAiw3RlukdAUU7jvVtePCqCKBTojSqGNWTaJoYahRKNZI30Qxhw4ti5IZ\nrhgJIrEimRiVUiiQQvvXF3vXHK2GWnLPXfHPfj/AYO4995w7vzkv+HDOnbmbqkKSpIX6O5e6AEnS\nlc0gkSQ1MUgkSU0MEklSE4NEktTEIJEkNeklSJJsTnI0yVSSXWOOfzDJnyR5Pcm/vZCxSW5Isj/J\nse75+j5qlST1qzlIkiwDHgS2ABuAu5NsmNPth8C/Ar50EWN3AQeqagI40O1Lki4zfXwi2QRMVdXx\nqnoD2AdsG+1QVaer6iDwtxcxdhuwt9veC9zZQ62SpJ4t7+Eca4CXRvZPALf2MHZVVU132y8Dq8ad\nIMlOYCfANddc8+EPfvCDF/jSkiSAQ4cO/WVVrVzo+D6CZNFVVSUZey+XqtoN7AYYDAY1HA7f0dok\n6UqX5MWW8X0sbZ0E1o3sr+3aWseeSrIaoHs+3VinJGkR9BEkB4GJJOuTrAC2A5M9jJ0EdnTbO4DH\neqhVktSz5qWtqjqb5D7gCWAZsKeqjiS5tzv+cJL3AUPg7wM/TvI5YENVvTpubHfqB4BHktwDvAjc\n1VqrJKl/uZpuI+81Ekm6eEkOVdVgoeP9ZbskqYlBIklqYpBIkpoYJJKkJgaJJKmJQSJJamKQSJKa\nGCSSpCYGiSSpiUEiSWpikEiSmhgkkqQmBokkqYlBIklqYpBIkpoYJJKkJgaJJKmJQSJJatJLkCTZ\nnORokqkku8YcT5KvdsefTnJL1/6BJIdHHq92f8+dJF9IcnLk2NY+apUk9Wt56wmSLAMeBO4ATgAH\nk0xW1bMj3bYAE93jVuAh4NaqOgpsHDnPSeDRkXFfqaovtdYoSVo8fXwi2QRMVdXxqnoD2Adsm9Nn\nG/D1mvEUcF2S1XP6/CLwfFW92ENNkqR3SB9BsgZ4aWT/RNd2sX22A9+c0/bZbilsT5Lre6hVktSz\ny+Jie5IVwMeB/znS/BBwMzNLX9PAl88zdmeSYZLhmTNnFr1WSdKb9REkJ4F1I/tru7aL6bMF+G5V\nnZptqKpTVfWjqvox8DVmltDeoqp2V9WgqgYrV65seBuSpIXoI0gOAhNJ1nefLLYDk3P6TAKf7L69\ndRvwSlVNjxy/mznLWnOuoXwCeKaHWiVJPWv+1lZVnU1yH/AEsAzYU1VHktzbHX8YeBzYCkwBrwGf\nmh2f5BpmvvH16Tmn/mKSjUABL4w5Lkm6DKSqLnUNvRkMBjUcDi91GZJ0RUlyqKoGCx1/WVxslyRd\nuQwSSVITg0SS1MQgkSQ1MUgkSU0MEklSE4NEktTEIJEkNTFIJElNDBJJUhODRJLUxCCRJDUxSCRJ\nTQwSSVITg0SS1MQgkSQ1MUgkSU0MEklSE4NEktSklyBJsjnJ0SRTSXaNOZ4kX+2OP53klpFjLyT5\nfpLDSYYj7Tck2Z/kWPd8fR+1SpL61RwkSZYBDwJbgA3A3Uk2zOm2BZjoHjuBh+Yc/0dVtXHOH5/f\nBRyoqgngQLcvSbrM9PGJZBMwVVXHq+oNYB+wbU6fbcDXa8ZTwHVJVs9z3m3A3m57L3BnD7VKknrW\nR5CsAV4a2T/RtV1onwKeTHIoyc6RPquqarrbfhlYNe7Fk+xMMkwyPHPmzELfgyRpgS6Hi+0fqaqN\nzCx/fSbJx+Z2qKpiJnDeoqp2V9WgqgYrV65c5FIlSXP1ESQngXUj+2u7tgvqU1Wzz6eBR5lZKgM4\nNbv81T2f7qFWSVLP+giSg8BEkvVJVgDbgck5fSaBT3bf3roNeKWqppNck+RagCTXAL8EPDMyZke3\nvQN4rIdaJUk9W956gqo6m+Q+4AlgGbCnqo4kubc7/jDwOLAVmAJeAz7VDV8FPJpktpbfqao/7I49\nADyS5B7gReCu1lolSf3LzOWHq8NgMKjhcDh/R0nSTyQ5NOfnFxflcrjYLkm6ghkkkqQmBokkqYlB\nIklqYpBIkpoYJJKkJgaJJKmJQSJJamKQSJKaGCSSpCYGiSSpiUEiSWpikEiSmhgkkqQmBokkqYlB\nIklqYpBIkpoYJJKkJr0ESZLNSY4mmUqya8zxJPlqd/zpJLd07euSfDvJs0mOJPm1kTFfSHIyyeHu\nsbWPWiVJ/VreeoIky4AHgTuAE8DBJJNV9exIty3ARPe4FXioez4L/Juq+m6Sa4FDSfaPjP1KVX2p\ntUZJ0uLp4xPJJmCqqo5X1RvAPmDbnD7bgK/XjKeA65KsrqrpqvouQFX9NfAcsKaHmiRJ75A+gmQN\n8NLI/gneGgbz9klyE/BzwJ+ONH+2Wwrbk+T6cS+eZGeSYZLhmTNnFvYOJEkLdllcbE/yHuD3gM9V\n1atd80PAzcBGYBr48rixVbW7qgZVNVi5cuU7Uq8k6Zw+guQksG5kf23XdkF9kryLmRD5RlX9/myH\nqjpVVT+qqh8DX2NmCU2SdJnpI0gOAhNJ1idZAWwHJuf0mQQ+2X176zbglaqaThLgt4Hnquo3Rgck\nWT2y+wngmR5qlST1rPlbW1V1Nsl9wBPAMmBPVR1Jcm93/GHgcWArMAW8BnyqG/7zwD8Hvp/kcNf2\n76vqceCLSTYCBbwAfLq1VklS/1JVl7qG3gwGgxoOh5e6DEm6oiQ5VFWDhY6/LC62S5KuXAaJJKmJ\nQSJJamKQSJKaGCSSpCYGiSSpiUEiSWpikEiSmhgkkqQmBokkqYlBIklqYpBIkpoYJJKkJgaJJKmJ\nQSJJamKQSJKaGCSSpCYGiSSpSS9BkmRzkqNJppLsGnM8Sb7aHX86yS3zjU1yQ5L9SY51z9f3Uask\nqV/NQZJkGfAgsAXYANydZMOcbluAie6xE3joAsbuAg5U1QRwoNuXJF1m+vhEsgmYqqrjVfUGsA/Y\nNqfPNuDrNeMp4Lokq+cZuw3Y223vBe7soVZJUs+W93CONcBLI/sngFsvoM+aecauqqrpbvtlYNW4\nF0+yk5lPObz3ve8lyQLegiRpofoIkkVXVZWkznNsN7AbYDAY1HA4fEdrk6QrXet/wPtY2joJrBvZ\nX9u1XUiftxt7qlv+ons+3UOtkqSe9REkB4GJJOuTrAC2A5Nz+kwCn+y+vXUb8Eq3bPV2YyeBHd32\nDuCxHmqVJPWseWmrqs4muQ94AlgG7KmqI0nu7Y4/DDwObAWmgNeAT73d2O7UDwCPJLkHeBG4q7VW\nSVL/UjX20sMVyWskknTxkhyqqsFCx/vLdklSE4NEktTEIJEkNTFIJElNDBJJUhODRJLUxCCRJDUx\nSCRJTQwSSVITg0SS1MQgkSQ1MUgkSU0MEklSE4NEktTEIJEkNTFIJElNDBJJUhODRJLUpClIktyQ\nZH+SY93z9efptznJ0SRTSXaNtP/nJD9I8nSSR5Nc17XflORvkhzuHg+31ClJWjytn0h2AQeqagI4\n0O2/SZJlwIPAFmADcHeSDd3h/cDPVNXPAn8OfH5k6PNVtbF73NtYpyRpkbQGyTZgb7e9F7hzTJ9N\nwFRVHa+qN4B93Tiq6o+q6mzX7ylgbWM9kqR3WGuQrKqq6W77ZWDVmD5rgJdG9k90bXP9S+B/jeyv\n75a1/jjJR89XQJKdSYZJhmfOnLnI8iVJrZbP1yHJk8D7xhy6f3SnqipJLaSIJPcDZ4FvdE3TwI1V\n9VdJPgz8QZIPVdWrc8dW1W5gN8BgMFjQ60uSFm7eIKmq2893LMmpJKurajrJauD0mG4ngXUj+2u7\nttlz/AvgV4BfrKrqXvN14PVu+1CS54GfBobzviNJ0juqdWlrEtjRbe8AHhvT5yAwkWR9khXA9m4c\nSTYDvw58vKpemx2QZGV3kZ4kNwMTwPHGWiVJi6A1SB4A7khyDLi92yfJ+5M8DtBdTL8PeAJ4Dnik\nqo504/8rcC2wf87XfD8GPJ3kMPC7wL1V9cPGWiVJiyDdatJVYTAY1HDo6pckXYwkh6pqsNDx/rJd\nktTEIJEkNTFIJElNDBJJUhODRJLUxCCRJDUxSCRJTQwSSVITg0SS1MQgkSQ1MUgkSU0MEklSE4NE\nktTEIJEkNTFIJElNDBJJUhODRJLUxCCRJDVpCpIkNyTZn+RY93z9efptTnI0yVSSXSPtX0hysvt7\n7YeTbB059vmu/9Ekv9xSpyRp8bR+ItkFHKiqCeBAt/8mSZYBDwJbgA3A3Uk2jHT5SlVt7B6Pd2M2\nANuBDwGbgd/sziNJusy0Bsk2YG+3vRe4c0yfTcBUVR2vqjeAfd24+c67r6per6q/AKa680iSLjOt\nQbKqqqa77ZeBVWP6rAFeGtk/0bXN+mySp5PsGVkam2/MTyTZmWSYZHjmzJkFvQlJ0sLNGyRJnkzy\nzJjHmz5VVFUBdZGv/xBwM7ARmAa+fJHjqardVTWoqsHKlSsvdrgkqdHy+TpU1e3nO5bkVJLVVTWd\nZDVweky3k8C6kf21XRtVdWrkXF8DvjXfGEnS5aV1aWsS2NFt7wAeG9PnIDCRZH2SFcxcRJ8E6MJn\n1ieAZ0bOuz3Ju5OsByaA7zTWKklaBPN+IpnHA8AjSe4BXgTuAkjyfuC3qmprVZ1Nch/wBLAM2FNV\nR7rxX0yykZklsReATwNU1ZEkjwDPAmeBz1TVjxprlSQtgsxc2rg6DAaDGg6Hl7oMSbqiJDlUVYOF\njveX7ZKkJgaJJKmJQSJJamKQSJKaGCSSpCYGiSSpiUEiSWpikEiSmhgkkqQmBokkqYlBIklqYpBI\nkpoYJJKkJgaJJKmJQSJJamKQSJKaGCSSpCYGiSSpSVOQJLkhyf4kx7rn68/Tb3OSo0mmkuwaaf8f\nSQ53jxeSHO7ab0ryNyPHHm6pU5K0eJY3jt8FHKiqB7qA2AX8u9EOSZYBDwJ3ACeAg0kmq+rZqvrV\nkX5fBl4ZGfp8VW1srE+StMhal7a2AXu77b3AnWP6bAKmqup4Vb0B7OvG/USSAHcB32ysR5L0DmsN\nklVVNd1tvwysGtNnDfDSyP6Jrm3UR4FTVXVspG19t6z1x0k+2linJGmRzLu0leRJ4H1jDt0/ulNV\nlaQWWMfdvPnTyDRwY1X9VZIPA3+Q5ENV9eqY+nYCOwFuvPHGBb68JGmh5g2Sqrr9fMeSnEqyuqqm\nk6wGTo/pdhJYN7K/tmubPcdy4J8CHx55zdeB17vtQ0meB34aGI6pbzewG2AwGCw0yCRJC9S6tDUJ\n7Oi2dwCPjelzEJhIsj7JCmB7N27W7cAPqurEbEOSld1FepLcDEwAxxtrlSQtgtYgeQC4I8kxZgLh\nAYAk70/yOEBVnQXuA54AngMeqaojI+fYzlsvsn8MeLr7OvDvAvdW1Q8ba5UkLYJUXT2rQYPBoIbD\nt6x+SZLeRpJDVTVY6Hh/2S5JamKQSJKaGCSSpCYGiSSpiUEiSWpikEiSmhgkkqQmBokkqYlBIklq\nYpBIkpoYJJKkJgaJJKmJQSJJamKQSJKaGCSSpCYGiSSpiUEiSWpikEiSmhgkkqQmTUGS5IYk+5Mc\n656vP0+/PUlOJ3nmQscn+XySqSRHk/xyS52SpMXT+olkF3CgqiaAA93+OP8N2Hyh45NsALYDH+rG\n/WaSZY21SpIWQWuQbAP2dtt7gTvHdaqq/wP88CLGbwP2VdXrVfUXwBSwqbFWSdIiWN44flVVTXfb\nLwOrehq/BnhqpN+Jru0tkuwEdna7r89dPlvCfgr4y0tdxGXCuTjHuTjHuTjnAy2D5w2SJE8C7xtz\n6P7RnaqqJLXQQhY6vqp2A7sBkgyrarDQGq4mzsU5zsU5zsU5zsU5SYYt4+cNkqq6/W1e/FSS1VU1\nnWQ1cPoiX/98408C60b6re3aJEmXmdZrJJPAjm57B/BYT+Mnge1J3p1kPTABfKexVknSImgNkgeA\nO5IcA27v9kny/iSPz3ZK8k3gT4APJDmR5J63G19VR4BHgGeBPwQ+U1U/uoB6dje+n6uJc3GOc3GO\nc3GOc3FO01ykasGXNSRJ8pftkqQ2BokkqclVEyRJNne3U5lKcr5f2F81xt12ZinecibJuiTfTvJs\nkiNJfq1rX4pz8XeTfCfJ97q5+I9d+5Kbi1lJliX5syTf6vaX5FwkeSHJ95Mcnv2qb69zUVVX/ANY\nBjwP3AysAL4HbLjUdS3ye/4YcAvwzEjbF4Fd3fYu4D912xu6OXk3sL6bq2WX+j30NA+rgVu67WuB\nP+/e71KciwDv6bbfBfwpcNtSnIuROfnXwO8A3+r2l+RcAC8APzWnrbe5uFo+kWwCpqrqeFW9Aexj\n5jYrV60af9uZJXfLmaqarqrvdtt/DTzHzF0QluJcVFX9v273Xd2jWIJzAZBkLfBPgN8aaV6Sc3Ee\nvc3F1RIka4CXRvbPe0uVq9zb3XLmqp+fJDcBP8fM/8SX5Fx0SzmHmflx7/6qWrJzAfwX4NeBH4+0\nLdW5KODJJIe620pBj3PReq8tXaaq2m5Zc6VJ8h7g94DPVdWrSX5ybCnNRc383mpjkuuAR5P8zJzj\nS2IukvwKcLqqDiX5hXF9lspcdD5SVSeT/ENgf5IfjB5snYur5ROJt1SZcaq71QxL6ZYzSd7FTIh8\no6p+v2teknMxq6r+L/BtZv4Mw1Kci58HPp7kBWaWuv9xkv/O0pwLqupk93waeJSZpare5uJqCZKD\nwESS9UlWMPO3TCYvcU2XwpK75UxmPnr8NvBcVf3GyKGlOBcru08iJPl7wB3AD1iCc1FVn6+qtVV1\nEzP/HvzvqvpnLMG5SHJNkmtnt4FfAp6hz7m41N8m6PFbCVuZ+cbO88D9l7qed+D9fhOYBv6WmTXM\ne4B/wMwfCDsGPAncMNL//m5ujgJbLnX9Pc7DR5hZ/30aONw9ti7RufhZ4M+6uXgG+A9d+5Kbiznz\n8guc+9bWkpsLZr7N+r3ucWT238c+58JbpEiSmlwtS1uSpEvEIJEkNTFIJElNDBJJUhODRJLUxCCR\nJDUxSCRJTf4/Mh2D3q1zqosAAAAASUVORK5CYII=\n",
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# %notebook inline\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"import matplotlib.animation as animation\n",
"\n",
"t_t= np.array(tt)\n",
"x_y= np.array(xx)\n",
"v_y= np.array(vv)\n",
"\n",
"fig, ax = plt.subplots()\n",
"ax.set_xlim(( 0, 500))\n",
"ax.set_ylim((-0.1,0.1))\n",
"plt.hlines(0, 0, 500, color='k', linestyle='-', linewidth=1)\n",
"l_x, = ax.plot(x_y, 'b', lw=2)\n",
"l_v, = ax.plot(v_y, 'r', lw=2)\n",
"plt.close\n",
"\n",
"def init():\n",
" l_x.set_data([], [])\n",
" l_v.set_data([], [])\n",
" return (l_x, l_v)\n",
"\n",
"def animate(i):\n",
" x = t_t[:i]\n",
" l_x.set_data(x, x_y[:i])\n",
" l_v.set_data(x, v_y[:i])\n",
" return (l_x, l_v)\n",
"\n",
"plt.rcParams['animation.html'] = 'html5'\n",
"\n",
"anim = animation.FuncAnimation(fig, animate, init_func=init,\n",
" frames=500, interval=50, blit=True)"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"scrolled": true
},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
""
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"anim"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## animate関数 \n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"plotsパッケージにあるanimate関数を使う．構文は以下の通りで，[]内に動画にしたい関数を定義し，tで時間を変えていく．\n",
"```maple\n",
"> with(plots): animate(plot, [sin(x-t),x=0..5*Pi], t=0..10);\n",
"```\n",
"\n",
"|{{attach_view(MapleCGplot2d9.png,)}}|\n",
"|:----|\n",
"\n",
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## リストに貯めて，display表示 \n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"おなじ動作を，display関数でオプションとしてinsequence=trueとしても可能．この場合は第一引数に入れるリスト([])に一連の画像を用意し，コマ\n",
"送りで表示させる．\n",
"```maple\n",
"> tmp:=[]: n:=10: for i from 0 to n do t:=i; tmp:=[op(tmp),\n",
"> plot(sin(x-t)+sin(x+t),x=0..5*Pi)]; end do:\n",
"> display(tmp,insequence=true);\n",
"```\n",
"\n",
"|{{attach_view(MapleCGplot2d10.png,)}}|\n",
"|:----|\n",
"\n",
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 凝った例 \n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"ヘルプにある凝った例．Fという動画のコマを吐く関数を用意する．これを，animate関数から適当な変数を入れて呼び出す．backgroundには動かない絵を指定\n",
"することができる．\n",
"```maple\n",
"> with(plottools,line): F := proc(t) plots[display](\n",
"> line([-2,0],[cos(t)-2,sin(t)],color=blue),\n",
"> line([cos(t)-2,sin(t)],[t,sin(t)],color=blue),\n",
"> plot(sin(x),x=0..t,view=[-3..7,-5..5]) ); end:\n",
"> animate(F,[theta],theta=0..2*Pi, background=plot([cos(t)-2,sin(t),t=0..2*Pi]),\n",
"> scaling=constrained,axes=none);\n",
"```\n",
"$$\n",
"$$\n",
"\n",
"|{{attach_view(MapleCGplot2d11.png,)}}|\n",
"|:----|\n",
"\n",
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## ファイルへの保存 \n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"animationなどのgif形式のplotを外部ファイルへ出力して表示させるには，以下の一連のコマンドのようにする．\n",
"```maple\n",
"> plotsetup(gif,plotoutput=file2): display(tmp,insequence=true);\n",
"> plotsetup(default):\n",
"```\n",
"こいつをquicktimeなどに食わせれば，Maple以外のソフトで動画表示が可能となる．3次元図形の標準規格であるvrmlも同じようにして作成することが可能(?vrml;参照)．\n"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.1"
},
"latex_envs": {
"LaTeX_envs_menu_present": true,
"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": {
"colors": {
"hover_highlight": "#DAA520",
"navigate_num": "#000000",
"navigate_text": "#333333",
"running_highlight": "#FF0000",
"selected_highlight": "#FFD700",
"sidebar_border": "#EEEEEE",
"wrapper_background": "#FFFFFF"
},
"moveMenuLeft": true,
"nav_menu": {
"height": "12px",
"width": "252px"
},
"navigate_menu": true,
"number_sections": true,
"sideBar": true,
"threshold": 4,
"toc_cell": true,
"toc_section_display": "block",
"toc_window_display": true,
"widenNotebook": false
}
},
"nbformat": 4,
"nbformat_minor": 2
}