理系は大学入ったらプログラムやったらいいよ、という話

理系学生のみなさんは大学初年度に関数電卓を買うと思います。
私も当時生協で売っていたカシオの関数電卓を買い、定期試験などでよく使用したのを覚えています。

でも普段の学習では関数電卓より、プログラミング（言語は何でもいいのですが）で計算したほうが色々はかどるぞ！ってのが今回の話です。
今もし大学1年の自分に何か伝えられるとしたら、「プログラミングと線形代数はしっかりやっとけ！」って言いたい・・(´・ω・｀)
僕がプログラミングに興味をもったのは、四年生の冬だったので遅かった・・・
もっと早く知っておけば～だったのにと思うことが、よくあるので、記事にしてみました。
（線形代数については機会があればまた話します）
以下、大学に入って早い時期にプログラミングを学ぶ利点を列挙します。言語はpython

利点①単純計算のミスが減る

具体例から示しましょう
たとえば水素原子の原子模型で、ボーア模型というものがあります。
これによるとエネルギー最安定状態（基底状態）の第1半径は
${ a_{n=1}=\frac{4\pi\epsilon_0\hbar^2}{me^2} }$
です。（要は水素原子核からだいたいこれくらいの位置に電子が存在する、という値）
SI単位系でこれがどれほどの値か計算してみましょう

手計算

f:id:Misaki_yuyyuyu:20170826200534j:plain
ボーア半径計算一発勝負！
・・途中プランク定数を書き間違えそうになってますが、かろうじて正解です。
ただこれ以上複雑な計算だと、保証はできない。。。といった印象
時間は物理定数調べる時間＋検算を合わせて、2，3分

プログラムで計算

#!/usr/bin/env python
# -*- coding: utf-8 -*-
#used packages
import numpy as np
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt
from pandas import *
from pylab import *
import scipy

#physical constants
e=1.602*10**(-19)
pi=np.pi
h=6.626*10**(-34)
hbar=h/2/pi
e0=8.854*10**(-12)
me=9.109*10**-31


print 4*pi*e0*hbar**2/(me*e**2)

この短いコードを走らせれば、5.29e-11m、つまり0.5Å（0.05nm）と正しい値がアウトプットされます。
実行にかかる時間はせいぜい1秒以内、式を入力する時間を合わせても10数秒

なんといっても良いところは、コードの再利用ができること。基底状態のエネルギーE_0
${ E_{n=1}=E_0=-\frac{me^4}{32\pi^2\epsilon_0^2\hbar^2} }$
を知りたければ、最後のprint以下を変えればいいだけ
この計算過程で唯一計算ミスを誘発しそうなのは、物理定数の入力のところですが、こういうコードを再利用することで入力を避け、計算ミスを減らすことができます。関数電卓にはない機能かな
あとこういう手計算って、だいたい本質的な問題ではないのでラクしちゃったほうがいいと思います（分野にもよりますが）
理系学生が計算で苦労するのは、三次元ラプラシアンの極座標表示の導出くらいでええんや・・・

利点②　グラフによる可視化ができる

たとえば学科の授業で苦労して水素原子の3dz^2波動関数を導出したとします。結果
${ \psi=\frac{1}{81\sqrt{2\pi}}(\frac{1}{a_0})^\frac{3}{2}(\frac{r}{a_0})^2e^{-r/3a_0}(\sqrt{3}cos^2\theta-\frac{1}{\sqrt{3}}) }$
という式を導出したとします。ではこの軌道はどのような空間分布をしているのだろうか？
関数の形からΦ依存性がないので、z軸に無限回転軸があるとか、θ=90に節みたいなところがあるとかはなんとなく分かるんですが、
それ以上は良くわからない。ということでマッピングしてみます。
波動関数は三次元の(x,y,z)各点に値を持つので視覚的な表示が難しいですが、ここではx=0としてyz平面による断面を表示してみましょう。

#!/usr/bin/env python
# -*- coding: utf-8 -*-
#used packages
import numpy as np
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt
import pylab
import scipy
#proportional to electron density
def P(r,theta):
    return r**2*np.exp(-r)*(3**0.5*(np.cos(theta)**2-3**(-0.5)))

#physical constants
e=1.602*10**(-19)
pi=np.pi
h=6.626*10**(-34)
hbar=h/2/pi
e0=8.854*10**(-12)
me=9.109*10**-31
j=1j


X=np.arange(-5,5,0.03)
Y=np.arange(-5,5,0.03)

#xy to rtheta
x,y=pylab.meshgrid(X,Y)
r=(x**2+y**2)**0.5
theta=np.arccos(x/r)
p=P(r,theta)**2

fig = plt.figure(figsize=(6,6))
ax = fig.add_subplot(111)

#plotting
im=plt.imshow(p,extent=(-5,5,-5,5))
plt.colorbar(im)
plt.axis('off')
plt.savefig('3dz2.png')
plt.show()

結果は以下です。規格化していないので、値自体に意味はないのですが、電子密度に比例した量となっています。
f:id:Misaki_yuyyuyu:20170826204837p:plain
wikipediaによると3dz軌道の等電子密度面は以下左上の形なので、計算はだいたいあっている模様

特に量子力学とかは計算が大変になりがちで、始めたころはそっちに意識がいっちゃうけど、プログラム使えれば大局的に考えられるようになる気がする。

利点③　代数計算・行列計算

計算にx,yなどの変数定数が含まれていても代数演算ができるのがsympyというパッケージ

代数計算

${ f(x)=\cos(x^2)\sin(x^2) }$
として二階微分f''(x)を求める。関数はテキトーだけど考えただけでもめんどくさそう・・・

#!/usr/bin/env python
# -*- coding: utf-8 -*-
#used packages
import sympy as sym

x =sym.Symbol('x')
print sym.diff(sym.sin(x**2)*sym.cos(x**2), x, 2)
# output value is 2*(-8*x**2*sin(x**2)*cos(x**2) - sin(x**2)**2 + cos(x**2)**2)

べんり

行列計算

#!/usr/bin/env python
# -*- coding: utf-8 -*-
#used packages
import sympy as sym

x =sym.Symbol('x')
y=sym.Symbol('y')
A =sym.Matrix([[1,x], [y,1]])
print A**2
# output is Matrix([[x*y + 1, 2*x], [2*y, x*y + 1]])