固有値と変換行列

最近はPCの性能が上がり、便利なフリーソフトが増えてきたので、PCを使って簡単にいろいろなことができるようになってきている。 様々な学習にもPCを利用すると、理解が深まることがある。 私は式変形にはmaximaを、数値計算にはRを使って、理解の助けにしている。 例として、角運動量をRを使っていじってみよう。

角運動量演算子は、Jzの固有関数を使うと行列表示できる。 これをRで表すと、以下のようになる。

mjz<-function(j){return( diag(seq(j,-j,-1)) )}
mjp<-function(j){
 jp<-matrix(0,nrow=2*j+1,ncol=2*j+1,byrow=TRUE)
 for(i in 1:(j*2)){jp[i,i+1]<-sqrt(i*(j*2+1-i))}
 return(jp)
}
mjm<-function(j){return( t(mjp(j)) )}
mjx<-function(j){return( (mjp(j)+mjm(j))/2 )}
mjy<-function(j){return( (mjp(j)-mjm(j))/2i )}
me<-function(j){return( diag(2*j+1) )}

ここで、mjpとmjmはそれぞれJ+とJ-を、meは2J+1次元の単位行列を表す。 Jxを対角化するには、例えばeigen(mjx(3/2))とすると良い。 valuesに対角成分に対応する固有値が、vectorsには変換係数が出てくる。 このような振る舞いを見てみると、いろいろと理解が深まる。

2020/1/29追記 2019/3/21に指摘したバグも含めて、単純化したものも公開しておきます。そのプログラムがすぐに見つからなくて、公開した方が自分が便利なので。

mjz<-function(j){ diag(j:-j,2*j+1) }
mjp<-function(j){ matrix(sqrt(cumsum(j:-j)*2)*!1:(2*j+1)^2%%(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) }

2020/3/17追記 さらに、シンプルに書くことができたので、更新しておきます。

mjp<-function(j){ diag(sqrt(diffinv(j:-j)*2))[-1,-2*j-2] }

androidとPCの間でのデータ転送
最近、タブレットを使う人が増えてきているように感じる。しかし、ipadとPCの間で、直接データの転送をするのは、なかなか面倒らしい。androidの場合には、どうなのか調べてみた。 linuxにUSBで接続すると、dmesgでは認識しているけど、マウントされない。調べてみると、androidはMTPプロトコルを使っているかららしい。そこで、linuxにjmtpfsをインストールして、android側でMTPを有効にすると、自動的にマウントされるようになり、データの転送も自由にできる。アンマウントが あまりうまく行かないけど。

データフレームの扱い

Rを使って結晶内の特定の距離の範囲内にある原子を数え上げるプログラムを作っていたのだが、実行が異様に遅い。いろいろと試してみると、データフレームを結合するために使っているrbindのところで、ほとんどの時間を使っていることが分かった。webで検索してみると、dplyrというlibraryのbind_rowsを使うと速くなるというので、試してみたら、なんとか許容範囲の実行時間になって来た。

しかし、もっと速くできるのでは無いかといろいろと試行錯誤していたら、かなり速くすることに成功した。expand.gridを使うと、指定したすべての組み合わせのデータフレームを作ってくれるので、それを利用した。極端な例でこれを説明しよう。

球の中にある格子点の数を数えるプログラムを作るとしよう。座標をふって球の中にある場合に、データフレームに加えていくということをする場合には、以下のようなプログラムになる。

d<-data.frame(x=c(),y=c(),z=c())
for(x in -10:10){
for(y in -10:10){
for(z in -10:10){
 if(x^2+y^2+z^2<100){
  d<-rbind(d,data.frame(x=x,y=y,z=z))
 }
}}}

このrbindが時間をかなり使っているようだ。これをexpand.gridで格子点を作って、それを一度に条件判定をするプログラムは、以下のようになる。

d<-expand.grid(x=-10:10,y=-10:10,z=-10:10)
d<-subset(d,x^2+y^2+z^2<100)

このときの実行時間は、数百倍も違う。こんなに違うのかと、少し驚いたが、今後Rを使うときには、このようなことを考慮しようと思う。

久々のlubuntuのインストール
SSDをつんだノートにubuntu18.04を入れたら、以外に苦労した。まず、USBメモリにisoを入れて、そこからブートしてインストールを開始した。しかし、swapパーティションを無くすのを忘れたと思ってと中で止めたのがいけなかったか、それ以降苦労することになった。 partitionerが途中で止まって動かないようになってしまったので、レスキューモードで立ち上げたら、最後までインストール出来た。そのためにユーザーの設定がされない状態になってしまった。そこで、同じくレスキューモードでユーザーを作成した。 後で調べたら、今のubuntuはswapパーティションが無くなって、swapファイルになっていたので、最初のインストールを止める必要が無かったことが分かった。 さらに、aptでupdateをすると、「ハッシュサムが適合しません」と出てしまう。そのようなバグも報告されていたようである。/var/lib/apt/lists/* を消去することによって、この問題も解決した。そして、upgradeをしたら、それなりに使えるようになってきた。

Rのrglを使った原子軌道の角度部分の絵

原子軌道の角度部分の立体的な絵を、ソフトを使って書いてみた。webを調べるとpovrayを使って書いている人が多かったが、特殊関数などの定義が面倒だし、Rのrglを使って書いた方が良い気がして、その方針でやってみた。ソースは以下のとおりです。

library(misc3d)
library(gsl)
theta<-function(x,y,z)acos(z/sqrt(x^2+y^2+z^2))
phi<-function(x,y,z)atan( (sqrt(x^2+y^2)-x)/y )*2
th<- 0:100/100*pi
ph<- 0:200/200*2*pi
fc <- function(x,y,z)ifelse(fr(theta(x,y,z),phi(x,y,z))<0,'blue','red')
fx <- function(t,p)abs(fr(t,p))*sin(t)*cos(p)
fy <- function(t,p)abs(fr(t,p))*sin(t)*sin(p)
fz <- function(t,p)abs(fr(t,p))*cos(t)
for(l in 0:5){
 for(m in 0:l){
  fr <- function(t,p)legendre_Plm(l,m,cos(t))*cos(m*p)
  parametric3d(fx, fy, fz, th, ph,color=fc)
  play3d(spin3d(c(0, 0, 1), 10),5)
}}

ライブラリとしては、特殊関数を使うためにrglを使って、プロットにはmisc3dのparametric3dを使った。極座標のphiは、tan(2phi)から計算するようにして、条件分岐を減らすことができた。とりあえず、h軌道ぐらいまで書いて、五秒間10rpmで回すようにしてみたが、分かりやすくなったかな。

ルジャンドル陪多項式

Rを使って、球面調和関数展開をしようと思ったら、ルジャンドル陪多項式が見当たらない。調べてみたらgslにあった。r-cran-gslをaptで入れたら、

library(gsl)
legendre_Plm(l,abs(m),cos(theta))*exp(-1i*m*phi)

という感じで使うことができた。これで、必要なら軌道の絵なども書けるようになった。しかし、普段使っているノートは、少し古いubuntuで、これにはパッケージが無くて使えないという罠があった。

raspberry pi3を通じて

windowsで制御している装置のデータを、LANを通じて別のPCから見れるようにする際に、そのwindows PCはインターネットに接続したくなかったので、間にraspberry pi3を入れてみた。

raspberry pi3とそのwindows PCを有線LANで繋ぐ計画だったのだが、装置の制御に有線LANを使っていたので、まずはIPを調べる必要があった。その結果、PCが1で、装置が10と20のようだった。重複しないように、raspberry piは53にしてみた。なんとなく素数にしたくなってしまう。routerは未設定にした。これで、制御PCと装置とraspberry piのみのLANが構築できる。

次は、windows10でwindows共有をするのだが、これは非常に苦労した。ネットの情報を元にして、いろいろなところをいじったり、何度も再起動していたら、なんとか共有をすることができた。特に最近のwindowsはwindows共有をするのが非常に難しい。

あとはraspberry piにその共有フォルダをマウントするように/etc/fstabを書き換えた。このとき、単にautoとしてもダメで、下記のようにしたらうまく行くようだ。

//ip.1/Data /home/pi/public cifs username=user,password=pass,noauto,x-systemd.automount 0 0

そして、raspberry piを無線LANで繋ぎ、共有したフォルダをftpでanonymousから読み取りのみでアクセスできるようにしておいた。本当はsftpにしたいのだが、windowsからアクセスするのに面倒らしいので、ftpで我慢している。

今の所、装置との干渉も気にならないようだし、自動mountもうまく働いているので、このまま問題が起きないと良いな。あとはUSBの穴を塞いで。

二系統のLAN

ラズベリーパイで二系統のLANに接続するために、内蔵のLANとUSB-LANを使おうと思っていたが、よくよく考えると、raspberry pi 3で内蔵LANとWifiを使えば良いことに気がついた。

まずは、sshを有効にするために、

sudo systemctl enable ssh

とする。そして、wifiのパスワードを

sudo sh -c 'wpa_passphrase SSID PASSPHRASE >> /etc/wpa_supplicant/wpa_supplicant.conf'

として設定する。そして、/etc/dhcpcd.confを編集する。

 interface eth0
 static ip_address=192.168.0.3/24
 interface wlan0
 static ip_address=172.16.0.3/24
 static routers=172.16.0.1
 static domain_name_servers=172.16.0.1

こんな感じです。これで、一応二系統のLANに接続できた。本当はさらにいろいろと設定して、やりたいことがあるけど、まだ接続するPCが無いのでできません。

micropythonでadt7310を読む

micropythonで温度測定するために、adt7310を使う方法を模索してみた。adt7310のpinとesp8266は以下のように接続した。

1 SCK  D5 14
2 DOUT D6 12
3 DIN  D7 13
4 CS   D8 15

まず、初期化は以下のように行う。

from machine import Pin,SPI
spi=SPI(baudrate=1000000,sck=Pin(14),mosi=Pin(13),miso=Pin(12),firstbit=SPI.MSB,polarity=1,phase=1)
cs=Pin(15, Pin.OUT)
cs(1)

ここで、polarityはclockが何もしていない時にはhighであることを示し、phaseはclockが元に戻るedgeで読み取りを行うことを示す。 adt7310のモードの設定は、以下のように行う。

cs(0)
#spi.write(bytearray([1<<3,0xa0])) #one shot
spi.write(bytearray([1<<3,0x80])) #continuous
cs(1)

そして、温度の読み取りは

cs(0)
spi.write(bytearray([0x40|2<<3]))
ret=spi.read(2)
cs(1)
(ret[0]*256+ret[1])/128

とすれば良い。本当は符号の処理も行わないといけないけど、零下になることはまず無いので、良しとしよう。bytearrayの扱いもよく分かっていないけど、まあ良いか。

optimizeとoptim

Rをつかって複雑な式でfittingをやることがある。線形の場合にはlmを、非線型の場合にはnlsを使えば良い。しかし、単純な関数で表せなかったり、媒介変数で表されるような場合には、これらは使えない。どうしたら良いかを考えていたら、fittingとの差を表すような関数を定義して、それを最小化すれば良いことに気がついた。例えば、xをax^2で近似するときの最適なaを求める場合には、差を表す関数を定義して、optimizeでその関数と変数の範囲を指定して、最小値を求めれば良い。

x<-0:10/10
f<-function(a){return(sum((x-a*x^2))^2)}
res<-optimize(f,c(0,2))
res$minimum

二変数の場合には、optimを使う。このとき、関数の引数はベクトルとして与え、optimには初期値をベクトルとして与える必要がある。

x<-0:10/10
f<-function(a){return(sum((x-a[1]*x^2-a[2]))^2)}
res<-optim(c(0,0),f)
res$par

ここでの例は、単純で通常の方法でもできるだろうが、複雑なfittingではこれが有効になった。