Misakichi’s ログblog

好きなものを紹介したり備忘録のため

pythonで特殊関数を使う(&カーブフィット)

pythonの科学技術計算に特化したモジュールscipyで使える特殊関数が豊富
Special functions (scipy.special) — SciPy v0.19.0 Reference Guide
これまで使ってなかったんだけど、フリーでこれだけ充実しているのはすごいです。その点いい時代ですね



みなさんは特殊関数と言われてなにを思い出しますか?
材料系の僕が最初に思いついたのは、エラーファンクション(誤差関数,erf)でした。
異種元素の拡散を扱うときによく出てくるやつです。
たとえば銅の表面に亜鉛メッキを行ってしばらく放置したとき、亜鉛がどこまで拡散するか?深さ方向の濃度分布は?
という疑問は拡散方程式に初期値を与えて、方程式を解くことで解決します。つまり拡散方程式
{\displaystyle
\frac{\partial C}{\partial t}=D\frac{\partial^2 C}{\partial x^2}
}
を束縛条件
{\displaystyle
t=0,x>0:C=0\\
t>0,x=0:C=C_0
}
で解きます。これは解析的にとけて、時刻t位置xの亜鉛濃度は
{\displaystyle
C(x,t)=C_0 erfc(\frac{x}{2\sqrt{Dt}})
}
となります。
たとえば拡散係数をD=1.0×10^{-9}m^2/sとするとこんな感じになりました。
tというのは亜鉛メッキを行った後の時間に相当します
f:id:Misaki_yuyyuyu:20170217011644p:plain


ちなみにここでのZn濃度Cは、時間と位置と拡散係数の関数です。
したがってある時間において、濃度と位置の関係を測定してやれば、逆にそこから拡散係数Dをもとめることもできちゃいます。
scipyにはカーブフィッティングが搭載されているみたいなので、それを利用して拡散係数を導出してみましょう。
ここでは
①拡散係数をD=1.0×10^{-9}m^2/sとして、true curveを作成
②測定による誤差を含んだデータを作成(noised data)
③誤差データをカーブフィッティングすることにより拡散係数Dの導出
という流れになっています。
測定から推測した拡散係数はD=1.02×10^{-9}m^2/sと真の値に近いものが得られました。

f:id:Misaki_yuyyuyu:20170218192757p:plain

固体中の拡散係数Dの導出は、だいたいこのようにモデルを立て、フィッティングを行うことで算出していまつ