kOde:常微分方程式

Python sf の kOde(..) 関数を使えば、非線形な常微分方程式でも、その数値解を簡単に求められます。といっても、scipy.integrate package の odeint(..) 関数を簡単に利用するためのインターフェースを被せているだけであり、数値解を求める作業自体は scipy.integrate.odeint(..) 関数が行っています。でも、インターフェースを被せることで、より単純に常微分方程式の数値解を得られるようにしています。多くの常微分方程式問題を one-liner で扱えるようにしています。

scipy package は Matlab に近い機能を実装している Python package です。ここにはMatlab に近い多くの機能が既に実装されています。その多くが Python でありながら Matlab に似せてあります。scipy.integrate.odeint(..) も そのような機能の一つです。でも これは時間に依存して変化する常微分方程式系でも扱えるようにしているぶん、記述が複雑になります。One-liner で扱うには、荷が重くなりがちです。

kOde(..) 関数は、常微分方程式系を時不変なものに限定することで、その扱いを単純にしました。二変数系の三次系ぐらいまでは one-liner で扱えるようにできました。

kOde(..) 関数のソース・コードはユーザー・カスタマイズの好例にもなっています。kOde(..) もソースを公開しています。40 行程度にすぎません。ユーザーの日常計算で、同じような計算が何度も続いて面倒だと感じるようになったら、ユーザー・カスタマイズに挑戦してみてください。小さな Python プログラムを書くだけです。

以下、「大振幅での振り子」「時不変な van der Poll系とリミット・サイクル」「水星の軌道」「特殊相対論での軌道」「時間に依存して変化する van der Poll系」の順序で Python sf の kOde(..), odeint(..) を適用していきます。

これらの常微分方程式系への適用例で Python sf one-liner の便利さが分かってもらえると思います。

kOde(..) ソルバー

Python sf では kOde(関数 f, 初期値、時間) と三つの引数を与えるだけで、次のような




任意の時不変な微分方程式の数値解を得られます。

ここで微分定式系は

  • 時不変:時間に依存しない系でなければならないが、非線形であってもかまわなない
  • x にはスカラー値変数だけでなく、任意長のベクトルを指定できる。これにより、任意階数の微分方程式も扱えることにる。

ことに御注意ください。

解が であることが分かっていますが、時刻 0 での初期値が 1 である



の微分方程式系を kOde(..) を使って下のように解く事ができます。時刻 0 での初期値を 1 とし、10 秒までの数値解を求めます。





Python sf one-liner
方程式:x' = x が、 時刻 0 秒での初期値が 1 としたときの, 10 秒までのデフォルト 51 点の解を求める
kOde(λ x:x, 1, 10s`)
===============================
[  1.00000000e+00   1.22140276e+00   1.49182470e+00   1.82211881e+00
   2.22554095e+00   2.71828187e+00   3.32011698e+00   4.05520006e+00
   4.95303255e+00   6.04964764e+00   7.38905634e+00   9.02501383e+00
   1.10231768e+01   1.34637386e+01   1.64446476e+01   2.00855380e+01
   2.45325315e+01   2.99641018e+01   3.65982367e+01   4.47011874e+01
   5.45981538e+01   6.66863359e+01   8.14508749e+01   9.94843236e+01
   1.21510428e+02   1.48413172e+02   1.81272258e+02   2.21406437e+02
   2.70426434e+02   3.30299594e+02   4.03428836e+02   4.92749095e+02
   6.01845106e+02   7.35095275e+02   8.97847400e+02   1.09663329e+03
   1.33943093e+03   1.63598464e+03   1.99819616e+03   2.44060231e+03
   2.98095841e+03   3.64095084e+03   4.44706741e+03   5.43166042e+03
   6.63424504e+03   8.10308522e+03   9.89713068e+03   1.20883828e+04
   1.47647841e+04   1.80337481e+04   2.20264697e+04]
---- ClTensor ----

exp(10)
===============================
22026.4657948


kOde(..) 関数が求解する変数にはベクトルも許されるので より一般的な N 次微分の系も扱えますし、M 変数の系も扱えます。系に合わせた f(..) 関数の形と 初期値リストが必要になるだけです。見た目だけから一次の微分しか扱えないと看做すのは誤りです。



のような二次微分の系でも下のように 長さ 2 のベクトル変数で扱えば一次の微分系に落とせます。


方程式:x'' = -x, 時刻 0 秒での初期値が [1,0], 10 秒までのデフォルト 51 点の解を求める
kOde(~[`Y,-`X], [1,0], 10s`)
===============================
[[ 1.          0.        ]
 [ 0.98006658 -0.19866933]
 [ 0.92106099 -0.38941834]
 [ 0.82533561 -0.56464247]
 [ 0.6967067  -0.71735609]
 [ 0.54030229 -0.84147097]
 [ 0.36235774 -0.93203907]
 [ 0.16996713 -0.98544971]
 [-0.02919953 -0.99957358]
 [-0.2272021  -0.97384761]
 [-0.41614683 -0.9092974 ]
 [-0.58850111 -0.80849637]
 [-0.7373937  -0.67546315]
 [-0.85688872 -0.51550134]
 [-0.9422223  -0.33498812]
 [-0.98999245 -0.14111999]
 [-0.99829473  0.05837415]
 [-0.96679814  0.2555411 ]
 [-0.89675836  0.44252043]
 [-0.79096766  0.61185787]
 [-0.65364357  0.75680246]
 [-0.49026077  0.87157572]
 [-0.30733283  0.95160201]
 [-0.1121525   0.99369093]
 [ 0.087499    0.99616453]
 [ 0.28366218  0.9589242 ]
 [ 0.46851665  0.88345458]
 [ 0.63469284  0.77276441]
 [ 0.77556582  0.63126656]
 [ 0.88551945  0.46460211]
 [ 0.9601702   0.27941545]
 [ 0.996542    0.08308937]
 [ 0.99318482 -0.11654922]
 [ 0.95023249 -0.31154136]
 [ 0.86939738 -0.49411332]
 [ 0.75390215 -0.65698655]
 [ 0.60835122 -0.79366779]
 [ 0.43854725 -0.898708  ]
 [ 0.25125978 -0.96791956]
 [ 0.05395538 -0.99854322]
 [-0.14550005 -0.98935812]
 [-0.33915485 -0.94073042]
 [-0.51928862 -0.85459878]
 [-0.67871998 -0.73439697]
 [-0.81109292 -0.58491708]
 [-0.91113015 -0.41211839]
 [-0.97484349 -0.22288984]
 [-0.99969289 -0.02477538]
 [-0.9846877   0.17432679]
 [-0.93042611  0.36647911]
 [-0.83907137  0.54402106]]
---- ClTensor ----

cos(10s`)
===============================
-0.839071529076

`print([kOde(λ *x:(x[1],-x[0]), [1,0], 10s`,K)[K][0] for K in range(10,100, 10)])
`print([kOde(λ *x:[x[1],-x[0]], [1,0], 10s`,K)[K][0] for K in range(10,100, 10)])
`print([kOde(~[`Y,-`X], [1,0], 10s`,K)[K][0] for K in range(10,100, 10)])
[-0.8390710328088764,
 -0.83907120930045886,
 -0.83907139125088892,
 -0.83907131434765236,
 -0.83907137312275759,
 -0.83907142801738699,
 -0.83907145718651821,
 -0.83907146653632281,
 -0.83907149133404857]
===============================
None



Python sf one-liner


kOde:常微分方程式

Python sf の kOde(..) 関数を使えば、非線形な常微分方程式でも、その数値解を簡単に計算させられます。といっても、scipy.integrate package の odeint(..) 関数を簡単に利用するためのインターフェースを被せているだけであり、数値解を求める作業自体は scipy.integrate.odeint(..) 関数が行っています。でも、インターフェースを被せることで、より単純に常微分方程式の数値解を得られるようにしています。多くの常微分方程式問題を one-liner で扱えるようにしています。

scipy package は Matlab に近い機能を実装している Python package です。多くのユーザーにとって Matlab に近い機能が既に実装されています。多くの機能実装が Python でありながら Matlab に似せてあります。scipy.integrate.odeint(..) も そのような機能の一つです。でも これは時間に依存して変化する常微分方程式系でも扱えるようにしているぶん、記述が複雑になります。One-liner で扱うには、荷が重くなりがちです。

kOde(..) 関数は、常微分方程式系を時不変なものに限定することで、その扱いを単純にしました。二変数系の三次系ぐらいまでは one-liner で扱えるようにできました。

kOde(..) 関数のソース・コードはユーザー・カスタマイズの好例にもなっています。kOde(..) もソースを公開しています。40 行程度にすぎません。ユーザーの日常計算で、同じような計算が何度も続いて面倒だと感じるようになったら、ユーザー・カスタマイズに挑戦してみてください。小さな Python プログラムを書くだけです。

以下、「大振幅での振り子」「時不変な van der Poll系とリミット・サイクル」「水星の軌道」「特殊相対論での軌道」「時間に依存して変化する van der Poll系」の順序で Python sf の kOde(..) を適用していきます。

これらの常微分方程式系への適用例で Python sf one-liner の便利さが分かってもらえると思います。

調和振動子での波動関数の変化

こんどは調和振動子の時間発展の様子を見てみましょう。expm(..) 関数により行列自体の exp を計算させられます。expm(i t H) により、調和振動子の時間発展を計算させられます。



Python sf ワンライナー
N=64; Px,X=:Px128,X128;H=Px^2+X^2;v=eigh(H)[1][:,0];ψm :=~[expm(`i H t) v for t in arSqnc(0,N,300/N)]

=:ψm;=:ψm; renderMtCplx(ψm)


上の Python sf 式で arSqnc(0,N,300/N) は「0 から始まる、長さ N で刻み幅が 300/N の数列」の意味です。300 の値は位相の時間変化が一周期となるように選んだ結果です。位相の時間変化が全ての場所で一様に変化していくことが分かります。エネルギーの固有ベクトルの時間変化ですから、波動関数の絶対値の分布は一定です。

今度は一次の固有ベクトルの波動関数の時間変化を見てみましょう。上の Python sf 式で eigh(H)[1][:,0] の部分を eigh(H)[1][:,1] に書き換えているだけです。



Python sf ワンライナー
N=64; Px,X=:Px128,X128;H=Px^2+X^2;v=eigh(H)[1][:,1];ψm :=~[expm(`i H t) v for t in arSqnc(0,N,300/N)]

=:ψm;=:ψm; renderMtCplx(ψm)


位相変化の速度が、エネルギー・レベルに比例して早くなっています。


最後に時刻 0 での波動関数の分布がδ関数のときの、(波動関数が非固有関数状態の時となります)時間変化を見てみましょう。このときは X の変移に比例するばね定数を 100 倍しておきます。H=P^2+X^2 のままでは、自由粒子のときと殆ど同じ時間変化になってしまうからです。なお、ここでは Px256,X256 のより詳細な位置演算子、運動量演算子を使っています。



Python sf ワンライナー
N=256;Px,X=:Px256,X256;H=Px^2+100 X^2;v=kzrs(N);v[100]=1;ψm :=~[expm(`i H t) v for t in arSqnc(0,N,30/N)]

=:ψm;=:ψm; renderMtCplx(ψm)


波動関数で、δ関数が広がりながら波が伝播しますが、同時に調和振動子での振動も見られます。

ちなみに自由粒子のとき Hamiltonian は「H=p^2 」となり、初期値がδ関数の波動関数は、下のような時間変化となります。



Python sf ワンライナー
N=256;Px,X=:Px256,X256;H=Px^2;v=kzrs(N);v[100]=1;ψm :=~[expm(`i H t) v for t in arSqnc(0,N,30/N)]

=:ψm; renderMtCplx(ψm)


時間変化により波が端に到達したとき、反対側から、その波が再出現しています。有限領域での fft により微分演算子を定めているからです。逆に言って、波動関数が端に到達しないあいだは、実際の動きに近いと推測されます。

最後に
ここで行ったような Hamitonian 演算子の作成や、その系の時間変化は任意の系について適用できます。御自分でも是非とも試してください。こんな計算は Mathematica や Matlab ではやる気にならないと思います。Python sf ならば、ここで見てきたように簡単です。Python sf の評価版でも機能制限はありません。最初に 5 秒のディレーが入るだけです。ぜひとも御試しください。

調和振動子演算子:行列と固有値

前に作った X128, Px128 行列を使って調和振動子を扱ってみましょう。調和振動子の Hamiltonian は「P,X=:Px128,X128; H=P^2+X^2」の Python sf 式となります。その固有値を並べたベクトルは下の Pythnon sf 式で計算できます。eigvalsh(..) は Hermite 行列の固有値を計算します。eigval(..) でも計算できますが、eigvalsh(..) だと、固有値を昇順に並べてくれるので便利です。(eig(..), eigh(..) だと固有値と固有ベクトルの両方を計算してくれます。)



Python sf ワンライナー
P,X=:Px128,X128; H=P^2+X^2;eigvalsh(H)
===============================
[  0.01562500+0.j   0.04687500+0.j   0.07812500+0.j   0.10937500+0.j
   0.14062500+0.j   0.17187500+0.j   0.20312500+0.j   0.23437500+0.j
   0.26562500+0.j   0.29687500+0.j   0.32812500+0.j   0.35937500+0.j
   0.39062500+0.j   0.42187500+0.j   0.45312500+0.j   0.48437500+0.j
   0.51562500+0.j   0.54687499+0.j   0.57812502+0.j   0.60937488+0.j
   0.64062545+0.j   0.67187313+0.j   0.70313101+0.j   0.73435301+0.j
   0.76568259+0.j   0.79668628+0.j   0.82851392+0.j   0.85821144+0.j
   0.89244686+0.j   0.91719197+0.j   0.95912003+0.j   0.97451577+0.j
   1.03028588+0.j   1.03790528+0.j   1.10691858+0.j   1.11075316+0.j
   1.18924082+0.j   1.19141961+0.j   1.27716008+0.j   1.27854550+0.j
   1.37052430+0.j   1.37148116+0.j   1.46919984+0.j   1.46990011+0.j
   1.57308300+0.j   1.57361672+0.j   1.68209544+0.j   1.68251404+0.j
   1.79617784+0.j   1.79651281+0.j   1.91528474+0.j   1.91555646+0.j
   2.03938084+0.j   2.03960309+0.j   2.16843835+0.j   2.16862074+0.j
   2.30243508+0.j   2.30258456+0.j   2.44135315+0.j   2.44147479+0.j
   2.58517798+0.j   2.58527556+0.j   2.73389758+0.j   2.73397392+0.j
   2.88750198+0.j   2.88755922+0.j   3.04598287+0.j   3.04602258+0.j
   3.20933326+0.j   3.20935657+0.j   3.37754720+0.j   3.37755490+0.j
   3.55061220+0.j   3.55061965+0.j   3.72852388+0.j   3.72854631+0.j
   3.91128599+0.j   3.91132345+0.j   4.09889507+0.j   4.09894790+0.j
   4.29134811+0.j   4.29141689+0.j   4.48864245+0.j   4.48872806+0.j
   4.69077575+0.j   4.69087934+0.j   4.89774588+0.j   4.89786898+0.j
   5.10955092+0.j   5.10969546+0.j   5.32618911+0.j   5.32635751+0.j
   5.54765881+0.j   5.54785410+0.j   5.77395842+0.j   5.77418440+0.j
   6.00508638+0.j   6.00534782+0.j   6.24104111+0.j   6.24134404+0.j
   6.48182088+0.j   6.48217303+0.j   6.72742377+0.j   6.72783514+0.j
   6.97784751+0.j   6.97833122+0.j   7.23308921+0.j   7.23366284+0.j
   7.49314503+0.j   7.49383263+0.j   7.75800952+0.j   7.75884485+0.j
   8.02767451+0.j   8.02870650+0.j   8.30212697+0.j   8.30342925+0.j
   8.58134470+0.j   8.58303359+0.j   8.86528667+0.j   8.86755753+0.j
   9.15385692+0.j   9.15709481+0.j   9.44559565+0.j   9.45283663+0.j
   9.71253877+0.j   9.78445533+0.j   9.99535429+0.j  10.39654480+0.j]
---- ClTensor ----


皆様は、上の固有値の羅列に何らかの規則性を見出せますでしょうか。…………….. そんな時はグラフを描かせてみるのが定石です。



Python sf ワンライナー
P,X=:Px128,X128; H=P^2+X^2;plotGr(eigvalsh(H))


全体としては自乗関数の形です。でも最初の方は直線になっているようです。大きな固有値になってくると階段が目立ってきます。下のように差分値を見てみると、その規則性がより具体的に見えてきます。



Python sf ワンライナー
P,X=:Px128,X128; H=P^2+X^2; v=eigvalsh(H);v[1:]-v[:-1]
===============================
[  3.12500000e-02+0.j   3.12500000e-02+0.j   3.12500000e-02+0.j
   3.12500000e-02+0.j   3.12500000e-02+0.j   3.12500000e-02+0.j
   3.12500000e-02+0.j   3.12500000e-02+0.j   3.12500000e-02+0.j
   3.12500000e-02+0.j   3.12500000e-02+0.j   3.12500000e-02+0.j
   3.12500000e-02+0.j   3.12500000e-02+0.j   3.12499998e-02+0.j
   3.12500011e-02+0.j   3.12499939e-02+0.j   3.12500292e-02+0.j
   3.12498602e-02+0.j   3.12505637e-02+0.j   3.12476825e-02+0.j
   3.12578747e-02+0.j   3.12220069e-02+0.j   3.13295742e-02+0.j
   3.10036962e-02+0.j   3.18276415e-02+0.j   2.96975170e-02+0.j
   3.42354189e-02+0.j   2.47451147e-02+0.j   4.19280561e-02+0.j
   1.53957402e-02+0.j   5.57701087e-02+0.j   7.61940175e-03+0.j
   6.90132948e-02+0.j   3.83458841e-03+0.j   7.84876550e-02+0.j
   2.17878586e-03+0.j   8.57404757e-02+0.j   1.38542194e-03+0.j
   9.19787964e-02+0.j   9.56859884e-04+0.j   9.77186788e-02+0.j
   7.00276236e-04+0.j   1.03182883e-01+0.j   5.33726559e-04+0.j
   1.08478715e-01+0.j   4.18604185e-04+0.j   1.13663794e-01+0.j
   3.34974183e-04+0.j   1.18771930e-01+0.j   2.71721175e-04+0.j
   1.23824384e-01+0.j   2.22240775e-04+0.j   1.28835263e-01+0.j
   1.82395719e-04+0.j   1.33814336e-01+0.j   1.49475715e-04+0.j
   1.38768596e-01+0.j   1.21635290e-04+0.j   1.43703196e-01+0.j
   9.75738692e-05+0.j   1.48622020e-01+0.j   7.63453740e-05+0.j
   1.53528059e-01+0.j   5.72403224e-05+0.j   1.58423654e-01+0.j
   3.97101390e-05+0.j   1.63310672e-01+0.j   2.33168657e-05+0.j
   1.68190624e-01+0.j   7.69855867e-06+0.j   1.73057301e-01+0.j
   7.45544919e-06+0.j   1.77904229e-01+0.j   2.24230519e-05+0.j
   1.82739681e-01+0.j   3.74640032e-05+0.j   1.87571615e-01+0.j
   5.28323084e-05+0.j   1.92400208e-01+0.j   6.87875931e-05+0.j
   1.97225559e-01+0.j   8.56064600e-05+0.j   2.02047688e-01+0.j
   1.03594798e-04+0.j   2.06866533e-01+0.j   1.23102095e-04+0.j
   2.11681941e-01+0.j   1.44539030e-04+0.j   2.16493656e-01+0.j
   1.68400066e-04+0.j   2.21301294e-01+0.j   1.95293426e-04+0.j
   2.26104317e-01+0.j   2.25982002e-04+0.j   2.30901983e-01+0.j
   2.61440504e-04+0.j   2.35693282e-01+0.j   3.02937157e-04+0.j
   2.40476833e-01+0.j   3.52153213e-04+0.j   2.45250744e-01+0.j
   4.11362057e-04+0.j   2.50012376e-01+0.j   4.83704899e-04+0.j
   2.54757996e-01+0.j   5.73627998e-04+0.j   2.59482189e-01+0.j
   6.87600030e-04+0.j   2.64176889e-01+0.j   8.35336485e-04+0.j
   2.68829653e-01+0.j   1.03198877e-03+0.j   2.73420475e-01+0.j
   1.30228225e-03+0.j   2.77915443e-01+0.j   1.68889407e-03+0.j
   2.82253078e-01+0.j   2.27086185e-03+0.j   2.86299384e-01+0.j
   3.23789708e-03+0.j   2.88500841e-01+0.j   7.24097415e-03+0.j
   2.59702145e-01+0.j   7.19165523e-02+0.j   2.10898965e-01+0.j
   4.01190511e-01+0.j]
---- ClTensor ----


最初は 3.125e-2 で均等に並んでいたのが、最後のほうでは「0.273, 0.001, 0.278, 0.002, 0.282, 0.002, 0.288, 0.007, 0.260, 0.072, 0.211, 0.4]と一つ置きに変化するようになってきます。私はなんでこんな傾向が出てくるのか分かっていません。この説明ができる方がいらしたら、是非とも教えてください。

ただ最初のほうで 3.12500000e-2 と八桁精度で一定間隔に並ぶのも凄いことです。調和振動子らしく感じます。eigvalsh(H)[0] が 0.01562500 ちょうど 3.12500000e-2/2 になっているのも調和振動子の性質を的確に再現しています。「P,X=:Px128,X128; H=P^2+X^2」のモデルでも原点付近では的確なモデルになっていると思われます。

こんどは固有行列を見ておきましょう

eigh(H)がタプルで表現された「固有値の集まりのベクトル」と「固有ベクトルの集まりの行列の」ペアを返します。そのタプルの後ろ側の行列から縦ベクトルを取り出すと、固有ベクトルになります。複素数値ベクトルなどので絶対値のベクトルにしてグラフを描かせます



Python sf ワンライナー
P,X=:Px128,X128; H=P^2+X^2;plotGr(abs(eigh(H)[1][:,0]))


まずは 0 番目の固有ベクトルです。それらしい形をしています。



Python sf ワンライナー
P,X=:Px128,X128; H=P^2+X^2;plotGr(abs(eigh(H)[1][:,1]))


次に 1 番目の固有ベクトルです。これも それなりに納得できる形です。

位置・運動量 演算子:行列

X 座標 [-1,1] の位置に対応する演算子行列を Python sf を使って下のように計算できます。



Python sf ワンライナー
N=64; X=kzrs(2N,2N);for k in range(2N):X[k,k] = (-N+0.5+k)/N;X128:=X
===============================
[[-0.9921875  0.         0.        ...,  0.         0.         0.       ]
 [ 0.        -0.9765625  0.        ...,  0.         0.         0.       ]
 [ 0.         0.        -0.9609375 ...,  0.         0.         0.       ]
 ..., 
 [ 0.         0.         0.        ...,  0.9609375  0.         0.       ]
 [ 0.         0.         0.        ...,  0.         0.9765625  0.       ]
 [ 0.         0.         0.        ...,  0.         0.         0.9921875]]
---- ClTensor ----


上の Python sf 式を、少し詳しく解説します。

この one-liner は四つ Python sf 式・文より成り立っています。「;」デリミタで区切られています。Python では「;」で区切ることで、複数の式文を一行で書けます。ループ文であっても、実行させるのが一つのの式・文だけならば一行で書けます。通常のプログラミング言語でのように書くと下のようになります。



Python sf ブロック
//@@
N=64
X=kzrs(2N,2N)

for k in range(2N): X[k,k] = (-N+0.5+k)/N

putPv(X,"X128")
//@@@


kzrs(,) は Python sf が用意しているゼロ行列生成関数です。putPv(,) は picklable な変数をファイル変数としてカレント・ディレクトリに書き出す関数です。上の one-liner では「X128:=X」と表現しています。Python sf のファイル変数は「=:」で読み戻すことができ、何度でも再利用できます。

このような if 文を含まない、単純に並べていくだけのプログラムならば,「;」で区切って一行の one-liner 書いても可読性の悪化は少しだけで済みます。「;」や「:」を強調表示してやれば可読性は殆ど同じです。One-liner にすることで copy & paste 操作が簡単になります。Python sf は command line program であり、上の one-liner は「Python -u -m sfPP &quoteN=64; X=kzrs(2N,2N);for k in range(2N):X[k,k] = (-N+0.5+k)/N;X128:=X&quote」だけで計算できます。「Python -u -m sfPP &quote..&quote」を追加して shell に実行させる editor scipt なぞ簡単に書けます。(このような vim script をここに示してあります。)

エディタ上で「copy and paste ができる」「script でコンピュータに計算処理させられる」「一部のパラメータだけを修正編集して再実行させられる」 Python sf 式 one-liner は数式を使って試行錯誤しながら計算を繰り返す理系のツールとして強力なものになります。
One-liner であることで、IDE とは異なり、前回使った変数でも再度宣言しなければ使えません。でも、これは欠点ではなく利点です。Copy and paste では再宣言の手間はありません。再宣言することで、以前のコンテキストの影響を受けなくなります。考慮すべき範囲をを現在扱っている one-liner だけに限定できます。一ヶ月後に、その one-liner を実行させても、同じ計算結果を得られます。理工系分野で、様々の計算式を駆使しながら、また試行錯誤しながら思考の階梯を上っていく過程で Python sf の one-liner が強力な助けになることが分かってもらえると思います。このような視点からも以下の記事を読み続けてもらうことを希望します。

ここで作った X128 ファイル変数行列は [-1,1] の X 軸上の位置を 128 このδ関数の集まりで近似したモデルとも見なせます。次のようなグラフを描かせてみれば、それが見えやすいと思います。



Python sf ワンライナー
X=:X128; plotGr(X[100,:])


plotGr(..) は Python sf のグラフ表示関数であり、X[100,:] は Pyhon numpy での「行列の 100 行目をベクトルとして取り出す」操作です。

X 位置に対応する運動量演算子行列を Python sf を使って下のように計算できます。運動量演算子は微分演算子であり、それは p 空間では「f(p)=p」のような直線関数です。ですから、位置δ関数を Fourier 変換して、 p を掛けて、Fourier 逆変換してやれば X 空間での微分演算子行列を得られます



Python sf ワンライナー
N=64;md=sc.fft;p=pi sc.arange((-N+0.5)/N, 1, 1/N);P=kzrs(2N,2N,complex);mtE=sc.eye(2N);for k in range(2N):P[:,k]=md.ifft(p md.fft(mtE[k,:])); Px128:=P
===============================
[[  0.00000000e+00 +0.00000000e+00j  -2.45436926e-02 +9.99799194e-01j
   -2.45436926e-02 +4.99598340e-01j ...,  -2.45436926e-02 -3.32730723e-01j
   -2.45436926e-02 -4.99598340e-01j  -2.45436926e-02 -9.99799194e-01j]
 [ -2.45436926e-02 -9.99799194e-01j  -1.56125113e-17 +3.20997027e-16j
   -2.45436926e-02 +9.99799194e-01j ...,  -2.45436926e-02 -2.49196293e-01j
   -2.45436926e-02 -3.32730723e-01j  -2.45436926e-02 -4.99598340e-01j]
 [ -2.45436926e-02 -4.99598340e-01j  -2.45436926e-02 -9.99799194e-01j
    0.00000000e+00 -1.55892011e-16j ...,  -2.45436926e-02 -1.98995002e-01j
   -2.45436926e-02 -2.49196293e-01j  -2.45436926e-02 -3.32730723e-01j]
 ..., 
 [ -2.45436926e-02 +3.32730723e-01j  -2.45436926e-02 +2.49196293e-01j
   -2.45436926e-02 +1.98995002e-01j ...,  -3.29597460e-17 +3.60673409e-16j
   -2.45436926e-02 +9.99799194e-01j  -2.45436926e-02 +4.99598340e-01j]
 [ -2.45436926e-02 +4.99598340e-01j  -2.45436926e-02 +3.32730723e-01j
   -2.45436926e-02 +2.49196293e-01j ...,  -2.45436926e-02 -9.99799194e-01j
   -3.46944695e-18 -8.71535931e-17j  -2.45436926e-02 +9.99799194e-01j]
 [ -2.45436926e-02 +9.99799194e-01j  -2.45436926e-02 +4.99598340e-01j
   -2.45436926e-02 +3.32730723e-01j ...,  -2.45436926e-02 -4.99598340e-01j
   -2.45436926e-02 -9.99799194e-01j  -3.46944695e-18 +1.74223158e-16j]]
---- ClTensor ----


Python sf で多く使う ClTensor の行列・ベクトルクラスではベクトル同士の積演算は内積 scalar 値になります。でも運動量演算子を作るときに p を掛けるのはベクトルの要素同士を掛けることを意味します。一方で numpy の array ベクトルでの積演算は、ベクトルの要素同士を掛け合わせたベクトルを意味します。また sc は Python sf での numpy packege を意味します。上の式は numpy package での演算処理を行わせているので、少し複雑な Python sf 式になっています。

運動量演算子行列は複素行列です。100 行目の横一行は下のようになります。この数値の羅列から意味を読み取れる方は殆どいないと思います。



Python sf ワンライナー
P=:Px128; P[100,:]
===============================
[ -2.45436926e-02 +2.99065760e-02j  -2.45436926e-02 +2.84529606e-02j
  -2.45436926e-02 +2.70797918e-02j  -2.45436926e-02 +2.57790465e-02j
  -2.45436926e-02 +2.45436926e-02j  -2.45436926e-02 +2.33675379e-02j
  -2.45436926e-02 +2.22451063e-02j  -2.45436926e-02 +2.11715348e-02j
  -2.45436926e-02 +2.01424880e-02j  -2.45436926e-02 +1.91540857e-02j
  -2.45436926e-02 +1.82028430e-02j  -2.45436926e-02 +1.72856186e-02j
  -2.45436926e-02 +1.63995711e-02j  -2.45436926e-02 +1.55421219e-02j
  -2.45436926e-02 +1.47109232e-02j  -2.45436926e-02 +1.39038301e-02j
  -2.45436926e-02 +1.31188770e-02j  -2.45436926e-02 +1.23542567e-02j
  -2.45436926e-02 +1.16083021e-02j  -2.45436926e-02 +1.08794707e-02j
  -2.45436926e-02 +1.01663303e-02j  -2.45436926e-02 +9.46754697e-03j
  -2.45436926e-02 +8.78187364e-03j  -2.45436926e-02 +8.10814083e-03j
  -2.45436926e-02 +7.44524776e-03j  -2.45436926e-02 +6.79215449e-03j
  -2.45436926e-02 +6.14787495e-03j  -2.45436926e-02 +5.51147048e-03j
  -2.45436926e-02 +4.88204400e-03j  -2.45436926e-02 +4.25873466e-03j
  -2.45436926e-02 +3.64071288e-03j  -2.45436926e-02 +3.02717575e-03j
  -2.45436926e-02 +2.41734273e-03j  -2.45436926e-02 +1.81045153e-03j
  -2.45436926e-02 +1.20575430e-03j  -2.45436926e-02 +6.02513835e-04j
  -2.45436926e-02 -1.34224200e-17j  -2.45436926e-02 -6.02513835e-04j
  -2.45436926e-02 -1.20575430e-03j  -2.45436926e-02 -1.81045153e-03j
  -2.45436926e-02 -2.41734273e-03j  -2.45436926e-02 -3.02717575e-03j
  -2.45436926e-02 -3.64071288e-03j  -2.45436926e-02 -4.25873466e-03j
  -2.45436926e-02 -4.88204400e-03j  -2.45436926e-02 -5.51147048e-03j
  -2.45436926e-02 -6.14787495e-03j  -2.45436926e-02 -6.79215449e-03j
  -2.45436926e-02 -7.44524776e-03j  -2.45436926e-02 -8.10814083e-03j
  -2.45436926e-02 -8.78187364e-03j  -2.45436926e-02 -9.46754697e-03j
  -2.45436926e-02 -1.01663303e-02j  -2.45436926e-02 -1.08794707e-02j
  -2.45436926e-02 -1.16083021e-02j  -2.45436926e-02 -1.23542567e-02j
  -2.45436926e-02 -1.31188770e-02j  -2.45436926e-02 -1.39038301e-02j
  -2.45436926e-02 -1.47109232e-02j  -2.45436926e-02 -1.55421219e-02j
  -2.45436926e-02 -1.63995711e-02j  -2.45436926e-02 -1.72856186e-02j
  -2.45436926e-02 -1.82028430e-02j  -2.45436926e-02 -1.91540857e-02j
  -2.45436926e-02 -2.01424880e-02j  -2.45436926e-02 -2.11715348e-02j
  -2.45436926e-02 -2.22451063e-02j  -2.45436926e-02 -2.33675379e-02j
  -2.45436926e-02 -2.45436926e-02j  -2.45436926e-02 -2.57790465e-02j
  -2.45436926e-02 -2.70797918e-02j  -2.45436926e-02 -2.84529606e-02j
  -2.45436926e-02 -2.99065760e-02j  -2.45436926e-02 -3.14498356e-02j
  -2.45436926e-02 -3.30933385e-02j  -2.45436926e-02 -3.48493659e-02j
  -2.45436926e-02 -3.67322318e-02j  -2.45436926e-02 -3.87587261e-02j
  -2.45436926e-02 -4.09486772e-02j  -2.45436926e-02 -4.33256766e-02j
  -2.45436926e-02 -4.59180192e-02j  -2.45436926e-02 -4.87599427e-02j
  -2.45436926e-02 -5.18932780e-02j  -2.45436926e-02 -5.53696833e-02j
  -2.45436926e-02 -5.92537156e-02j  -2.45436926e-02 -6.36271305e-02j
  -2.45436926e-02 -6.85950256e-02j  -2.45436926e-02 -7.42948179e-02j
  -2.45436926e-02 -8.09097113e-02j  -2.45436926e-02 -8.86895090e-02j
  -2.45436926e-02 -9.79839134e-02j  -2.45436926e-02 -1.09298027e-01j
  -2.45436926e-02 -1.23389475e-01j  -2.45436926e-02 -1.41448786e-01j
  -2.45436926e-02 -1.65460136e-01j  -2.45436926e-02 -1.98995002e-01j
  -2.45436926e-02 -2.49196293e-01j  -2.45436926e-02 -3.32730723e-01j
  -2.45436926e-02 -4.99598340e-01j  -2.45436926e-02 -9.99799194e-01j
  -1.56125113e-17 +6.27969927e-17j  -2.45436926e-02 +9.99799194e-01j
  -2.45436926e-02 +4.99598340e-01j  -2.45436926e-02 +3.32730723e-01j
  -2.45436926e-02 +2.49196293e-01j  -2.45436926e-02 +1.98995002e-01j
  -2.45436926e-02 +1.65460136e-01j  -2.45436926e-02 +1.41448786e-01j
  -2.45436926e-02 +1.23389475e-01j  -2.45436926e-02 +1.09298027e-01j
  -2.45436926e-02 +9.79839134e-02j  -2.45436926e-02 +8.86895090e-02j
  -2.45436926e-02 +8.09097113e-02j  -2.45436926e-02 +7.42948179e-02j
  -2.45436926e-02 +6.85950256e-02j  -2.45436926e-02 +6.36271305e-02j
  -2.45436926e-02 +5.92537156e-02j  -2.45436926e-02 +5.53696833e-02j
  -2.45436926e-02 +5.18932780e-02j  -2.45436926e-02 +4.87599427e-02j
  -2.45436926e-02 +4.59180192e-02j  -2.45436926e-02 +4.33256766e-02j
  -2.45436926e-02 +4.09486772e-02j  -2.45436926e-02 +3.87587261e-02j
  -2.45436926e-02 +3.67322318e-02j  -2.45436926e-02 +3.48493659e-02j
  -2.45436926e-02 +3.30933385e-02j  -2.45436926e-02 +3.14498356e-02j]
---- ClTensor ----


でもグラフにしてやると、その意味が具体的に見えてきます。



Python sf ワンライナー
P=:Px128; plotGr(abs(P[100,:]))


δ関数の微分の近似らしくなってきます。

Python sf の renderMtCplx(..) 関数を使うと複素数値行列全体の複素数分布の様子も見ることが可能になる。高さで絶対値を表示すると同時に、Red(正の実数軸に対応します),Green,Blue の色の混ざり具合で位相回転を表すからです。



Python sf ワンライナー
P=:Px128; renderMtCplx(P)


ここで作った X128, Px128 行列を使うと、量子力学の様々の問題を具体的に計算可能な Python sf 式として表現できます。

パラメータ N を買えて X64,Px64, X256,P256 の行列も作っておきましょう。ラフでも早い計算をさせたいとき、時間をかけても精度のよい計算をさせたいときなどで便利に使い分けられます。



N=32;md=sc.fft;p=pi sc.arange((-N+0.5)/N, 1, 1/N);P=kzrs(2N,2N,complex);mtE=sc.eye(2N);for k in range(2N):P[:,k]=md.ifft(p md.fft(mtE[k,:])); Px64:=P

N=32; X=kzrs(2N,2N);for k in range(2N):X[k,k] = (-N+0.5+k)/N;X64:=X

N=128;md=sc.fft;p=pi sc.arange((-N+0.5)/N, 1, 1/N);P=kzrs(2N,2N,complex);mtE=sc.eye(2N);for k in range(2N):P[:,k]=md.ifft(p md.fft(mtE[k,:])); Px256:=P

N=128; X=kzrs(2N,2N);for k in range(2N):X[k,k] = (-N+0.5+k)/N;X256:=X


量子論での行列による力学

Python sf を使えば、量子力学での行列による力学を実際に計算可能です。教科書に文字で書いてあるだけの式を、数値実験可能にものにできます。以下、その様子を見ていきましょう

量子力学

量子力学では「波でもあり粒子でもある」とよく言われます。ふざけた説明です。こんなのは説明とは言えません。波の概念と粒子の概念は両立しません。