Rとjuliaとpythonとrubyで

• November 21, 2022

様々な言語で角運動量演算子
角運動量演算子の行列を求めるプログラムを4つの言語で書いてみたので，メモがてら公開しておこうと思う．言語名と角運動量演算子で検索してもあまりヒットしないので，誰かの役に立つかもしれないし． Rについては，以前もブログに書いたことがあるが，j=0の場合にも対応できるように一部を書き直した．

mjz<-function(j) diag(j:-j,2*j+1)
mjp<-function(j) matrix(diag(sqrt(diffinv(j:-j)*2))[-1,-2*j-2],2*j+1)
mjm<-function(j) t(mjp(j))
mjx<-function(j) (mjp(j)+mjm(j))/2
mjy<-function(j) (mjp(j)-mjm(j))/2i
me<-function(j) diag(2*j+1)

juliaは，比較的新しい言語だが，数値計算には便利かも知れないので，書いてみた．いつかこのルーチンを使ってみたい．

using LinearAlgebra 
mjz(j)=diagm(j:-1:-j)
mjp(j)=diagm(1=>2j:-2:2-2j|>cumsum.|>sqrt)
mjm(j)=mjp(j)'
mjx(j)=(mjp(j)+mjm(j))/2
mjy(j)=(mjp(j)-mjm(j))/2im
me(j)=I(2j+1)

pythonはなぜか人気が衰えない言語で，書きにくいけど需要があるかも知れないので，書いてみた．

import numpy as np 
def mjz(j):return np.diag(j-np.arange(2*j+1))
def mjp(j):return np.diag(np.sqrt(np.cumsum((j-np.arange(2*j))*2)),1)
def mjm(j):return mjp(j).T
def mjx(j):return (mjp(j)+mjm(j))/2
def mjy(j):return (mjp(j)-mjm(j))/2j
def me(j):return np.identity(int(2*j+1))

rubyは数値計算にはあまり向かない言語だと思っていたが，numoを使えば数値計算もできそうなので，書いてみた．

require "numo/linalg"
def mjz(j) (j-Numo::DFloat[0..2*j]).diag end
def mjp(j) Numo::NMath.sqrt((j-Numo::DFloat[0...2*j]).cumsum*2).diag(1) end
def mjm(j) mjp(j).transpose end
def mjx(j) (mjp(j)+mjm(j))/2 end
def mjy(j) (mjp(j)-mjm(j))/2i end
def me(j) Numo::DFloat.eye(2*j+1) end

Rはそのままで実行可能だが，それ以外の3つの言語では，何らかのライブラリを組み込まないといけないのが，少し面倒だ．Rやjuliaは半整数をあまり意識しなくても動くが，pythonとrubyは，rangeが半整数だと動かないので，2*jが整数であることを生かして少し工夫しないといけなかった．R以外の3つの言語では，対角要素の一つ上に成分を並んだ行列を簡単に作れるので，上昇演算子を簡単に書けるが，Rだとそこに工夫が必要で，もう少し良い手法もあるかも知れない．