Skip to content

Instantly share code, notes, and snippets.

@genkuroki

genkuroki/2017-06-17-232248octave.ipynb

Last active July 25, 2017 13:09
Show Gist options
  • Save genkuroki/1024efda1e084a31cbda0cf5fa3d35f7 to your computer and use it in GitHub Desktop.
Save genkuroki/1024efda1e084a31cbda0cf5fa3d35f7 to your computer and use it in GitHub Desktop.
Download ZIP
2017-06-17-232248octave.ipynb Comparison of continuous-time and discrete-time open Toda lattices
Display the source blob
Display the rendered blob
Raw
2017-06-17-232248octave.ipynb
{
"cells": [
{
"cell_type": "markdown",
"metadata": {
"collapsed": false
},
"source": [
"#オープン戸田格子とその離散化の解の一致\n",
"\n",
"2002年11月にやった計算の紹介\n",
"\n",
"黒木 玄"
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": false
},
"source": [
"##函数の定義\n",
"\n",
"toda_open(x) はLax形式でのオープン戸田格子の定義.\n",
"\n",
"sqreshape(x) は正方行列をベクトル化したものを正方行列に戻すための函数.\n",
"\n",
"##オープン戸田格子の定義\n",
"\n",
"オープン戸田格子の運動方程式は $x_1,\\ldots,x_n$ に関する次の微分方程式で定義される:\n",
"\\begin{align*}\n",
"& \\ddot x_1 = e^{x_2-x_1},\\\\\n",
"& \\ddot x_2 = e^{x_3-x_2}-e^{x_2-x_1},\\\\\n",
"& \\quad\\cdots\\cdots\\cdots\\cdots \\\\\n",
"& \\ddot x_{n-1} = e^{x_n-x_{n-1}}-e^{x_{n-1}-x_{n-2}},\\\\\n",
"& \\ddot x_n = -e^{x_n-x_{n-1}}.\n",
"\\end{align*}\n",
"これのLax形式について説明しよう.\n",
"簡単のため $4\\times4$ で説明する.\n",
"$$\n",
"p_i=\\dot x_i,\\quad q_i=e^{-(x_i-x_{i+1})/2}\n",
"$$\n",
"とおき, 3重対角行列 $L$ を次のように定める:\n",
"\\begin{align*}\n",
"L=\n",
"\\begin{bmatrix}\n",
"p_1 & q_1 & 0 & 0 \\\\\n",
"q_1 & p_2 & q_2 & 0 \\\\\n",
"0 & q_2 & p_3 & q_3 \\\\\n",
"0 & 0 & q_3 & p_4 \\\\\n",
"\\end{bmatrix}.\n",
"\\end{align*}\n",
"行列 $M=M(L)$ を次のように定める:\n",
"\\begin{align*}\n",
"M\n",
"&=\n",
"\\frac12\\text{($L$ の対角部分)}\n",
"+\\text{($L$ の上三角部分)}\n",
"\\\\ &=\n",
"\\begin{bmatrix}\n",
"p_1/2 & q_1 & 0 & 0 \\\\\n",
"0 & p_2/2 & q_2 & 0 \\\\\n",
"0 & 0 & p_3/2 & q_3 \\\\\n",
"0 & 0 & 0 & p_4/2 \\\\\n",
"\\end{bmatrix}.\n",
"\\end{align*}\n",
"ただし, 上三角部分に対角部分は含まれないものとする.\n",
"\n",
"このとき次の微分方程式はオープン戸田格子と同値である:\n",
"$$\n",
"\\frac{dL}{dt}=ML-LM.\n",
"$$\n",
"これをオープン戸田格子のLax形式と呼ぶ.\n",
"\n",
"函数 toda_open(x) ではLax形式でオープン戸田格子の微分方程式を定義している.\n",
"\n",
"`triu(ones(n,n),1).* L` は対角成分を含まない $L$ の上三角部分になる."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": false
},
"outputs": [
],
"source": [
"## open Toda lattice\n",
"##\n",
"## Example: To integrate the open Toda lattice from t = 0 to t = 1,\n",
"##\n",
"## solution = lsode('toda_open', vec(L0), 0:1)\n",
"## L1 = sqreshape(sol(2,:))\n",
"\n",
"function xdot = toda_open(x)\n",
"\n",
" L = sqreshape(x);\n",
" M = diag(diag(L))/2 + triu(ones(rows(L),columns(L)), 1).* L;\n",
" Ldot = M * L - L * M;\n",
" xdot = reshape(Ldot, rows(x), columns(x));\n",
"\n",
"endfunction"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"collapsed": false
},
"outputs": [
],
"source": [
"## Reshape rectangulat matrix to square matrix\n",
"##\n",
"## retval = sqreshape(x)\n",
"##\n",
"## x = matrix\n",
"## retval = square matrix\n",
"\n",
"function retval = sqreshape(x)\n",
"\n",
" size = prod(size(x));\n",
" n = ceil(sqrt(size));\n",
" retval = reshape(x, n, n);\n",
"\n",
"endfunction"
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": false
},
"source": [
"##戸田格子の初期条件 L0 を CCC と定義\n",
"\n",
"等間隔で質点が並んでおり、運動量は -2, -1, 1, 2."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"CCC =\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
" -2 1 0 0\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
" 1 -1 1 0\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
" 0 1 1 1\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
" 0 0 1 2\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n"
]
}
],
"source": [
"CCC = [\n",
" -2, 1, 0, 0;\n",
" 1, -1, 1, 0;\n",
" 0, 1, 1, 1;\n",
" 0, 0, 1, 2\n",
"];\n",
"CCC"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"L0 =\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
" -2 1 0 0\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
" 1 -1 1 0\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
" 0 1 1 1\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
" 0 0 1 2\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n"
]
}
],
"source": [
"L0=CCC"
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": false
},
"source": [
"##X0はオープン戸田格子の離散化の初期条件"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"X0 =\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
" 0.26161 0.38561 0.43429 0.24631\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
" 0.38561 1.08150 1.93478 1.41954\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
" 0.43429 1.93478 5.93631 5.94709\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
" 0.24631 1.41954 5.94709 10.46386\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n"
]
}
],
"source": [
"X0 = expm(L0)"
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": false
},
"source": [
"##L0とX0の固有値を計算"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"ans =\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
" -2.70493\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
" -0.82667\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
" 0.82667\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
" 2.70493\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n"
]
}
],
"source": [
"eig(L0)"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"ans =\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
" 14.953210\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
" 2.285684\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
" 0.066875\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
" 0.437506\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n"
]
}
],
"source": [
"eig(X0)"
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": false
},
"source": [
"##オープン戸田格子の離散化の時刻を1進める"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"X1 =\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
" 1.78285 2.78489 2.17629 0.94154\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
" 2.78489 5.95483 5.55793 2.88354\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
" 2.17629 5.55793 6.18306 4.02311\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
" 0.94154 2.88354 4.02311 3.82254\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n"
]
}
],
"source": [
"X1 = chol(X0) * chol(X0).'"
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": false
},
"source": [
"##オープン戸田格子の連続時間を1進める"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"sol =\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
" Columns 1 through 7:\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
" -2.00000 1.00000 0.00000 0.00000 1.00000 -1.00000 1.00000\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
" -0.52601 1.40051 0.00000 0.00000 1.40051 0.04909 1.94886\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
" Columns 8 through 14:\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
" 0.00000 0.00000 1.00000 1.00000 1.00000 0.00000 0.00000\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
" 0.00000 0.00000 1.94886 -0.04909 1.40051 0.00000 0.00000\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
" Columns 15 and 16:\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
" 1.00000 2.00000\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
" 1.40051 0.52601\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n"
]
}
],
"source": [
"sol = lsode('toda_open', vec(L0), 0:1)"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"L1 =\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
" -0.52601 1.40051 0.00000 0.00000\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
" 1.40051 0.04909 1.94886 0.00000\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
" 0.00000 1.94886 -0.04909 1.40051\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
" 0.00000 0.00000 1.40051 0.52601\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n"
]
}
],
"source": [
"L1 = sqreshape(sol(2,:))"
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": false
},
"source": [
"##連続時間版と離散時間版の比較\n",
"\n",
"expm(L1)とX1は一致する.\n",
"\n",
"このようにLax形式の微分方程式で定義されたオープン戸田格子の時間発展の整数時刻での解の値のexponentialは、コレスキー分解で定義した離散オープン戸田格子の値に一致する."
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"ans =\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
" 1.78285 2.78489 2.17629 0.94154\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
" 2.78489 5.95483 5.55793 2.88354\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
" 2.17629 5.55793 6.18306 4.02311\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
" 0.94154 2.88354 4.02311 3.82254\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n"
]
}
],
"source": [
"expm(L1)"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"X1 =\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
" 1.78285 2.78489 2.17629 0.94154\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
" 2.78489 5.95483 5.55793 2.88354\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
" 2.17629 5.55793 6.18306 4.02311\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
" 0.94154 2.88354 4.02311 3.82254\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n"
]
}
],
"source": [
"X1"
]
},
{
"cell_type": "code",
"execution_count": 0,
"metadata": {
"collapsed": false
},
"outputs": [
],
"source": [
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Octave",
"language": "octave",
"name": "octave"
},
"language_info": {
"file_extension": ".m",
"help_links": [
{
"text": "GNU Octave",
"url": "https://www.gnu.org/software/octave/support.html"
},
{
"text": "Octave Kernel",
"url": "https://github.com/Calysto/octave_kernel"
},
{
"text": "MetaKernel Magics",
"url": "https://github.com/calysto/metakernel/blob/master/metakernel/magics/README.md"
}
],
"mimetype": "text/x-octave",
"name": "octave",
"version": "4.0.0"
}
},
"nbformat": 4,
"nbformat_minor": 0
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment