2015年12月25日金曜日

熱伝導方程式を陰解法で解いてみるで

熱伝導方程式の数値解法は割と現代(近代)科学の真骨頂やと思ったのでまとめるで。 こういう風に数式からプログラムまでなんとなくまとめてある資料はなかったので大学生とかは重宝すると思ってるで。 けどコピーはいかんで。

熱伝導方程式

一次元の熱伝導方程式は次式で表される。
\begin{equation} \frac{\partial T}{\partial t} = \alpha \frac{\partial^2T}{\partial^2x} \end{equation} $T$は温度、$\alpha$は温度伝導率であり、次式より求まる。 \begin{equation} \alpha=\frac{\lambda}{\rho c_{p}} \end{equation} $\lambda$は熱伝導率、$\rho$は密度、$c_{p}$は定圧比熱である。

偏導関数の近似

差分法では偏導関数を前進差分、後退差分で次のように近似する。 \begin{equation} \frac{\partial f}{\partial x}=\frac{f_{i+1}-f_{i}}{\Delta x} \end{equation} \begin{equation} \frac{\partial f}{\partial x}=\frac{f_{i}-f_{i-1}}{\Delta x} \end{equation} 2階の偏導関数は1階の偏導関数であるから前進差分と後退差分を用いれば次式を得られる。 \begin{equation} \frac{\partial^2 f}{\partial^2 x} = \frac{f_{i+1}-2f_{i}+f_{i-1}}{\Delta x^2} \end{equation}

陽解法

前進差分を用いて数値解を求める方法は陽解法と呼ばれる。 $T_{x,t}$として陽解法による差分近似式は次式になる。 \begin{equation} \frac{T_{i,j+1} - T_{i,j}}{\Delta t} = \alpha\frac{T_{i+1,j}-2T_{i,j}+T_{i-1,j}}{\Delta x^2} \end{equation} $T_{i,j+1}$について解くと次式になる。ここで$\gamma$は$\frac{\Delta t}{\Delta x^2}$となる。 \begin{equation} T_{i,j+1}=\alpha \gamma T_{i+1,j} + (1-2 \alpha \gamma)T_{i,j}+\alpha \gamma T_{i-1,j} \end{equation}

陰解法

後退差分を用いて数値解を求める方法は陰解法と呼ばれる。 差分近似式は次式となる。 \begin{equation} \frac{T_{i,j+1}-T_{i,j}}{\Delta t} = \alpha \frac{T_{i+1,j+1}-2T_{i,j+1}+T_{i-1,j+1}}{\Delta x^2} \end{equation} これを整理すると次式になる。 \begin{equation} -\alpha \gamma T_{i+1,j+1}+(1+2 \alpha \gamma)T_{i,j+1} - \alpha \gamma T_{i-1,j+1} = T_{i,j} \end{equation} 陽解法とは異なり$T_{i,j}$から、次の時刻での$T_{i-1,j+1},T_{i,j+1},T_{i+1,j+1}$の3点を求める。 前式を行列で表すと次の式のようになる。 \[ \left( \begin{array}{ccccc} 1+2 \alpha \gamma & -\alpha \gamma & 0 & \ldots & 0 \\ -\alpha \gamma & 1+2 \alpha \gamma & -\alpha \gamma & \ldots & 0 \\ & \vdots & & \ddots & \vdots \\ 0 & \ldots & -\alpha \gamma & 1+2 \alpha \gamma & -\alpha \gamma \\ 0 & \ldots & 0 & -\alpha \gamma & 1+2 \alpha \gamma \end{array} \right) \left( \begin{array}{c} T_{1,j+1} \\ T_{2,j+1} \\ \vdots \\ T_{M-1,j+1} \\ T_{M,j+1} \end{array} \right) = \left( \begin{array}{c} T_{1,j} \\ T_{2,j} \\ \vdots \\ T_{M-1,j} \\ T_{M,j} \end{array} \right) \] 前式の連立一次方程式を数値解法によって解く。このとき行列の$T_{i,j+1}$と$T_{M,j+1}$には境界条件を設定する。

連立一次方程式の数値解法

熱伝導方程式を陰解法で解くためには連立一次方程式を解く必要がある。 連立一次方程式を解く方法は数多くあるが、ここでは反復法のヤコビ反復法とガウス・サイデル法について記載する。 以下の三元連立一次方程式を考える。 \begin{equation} \left( \begin{array}{ccc} a_{11} & a_{12} & a_{13} \\ a_{21} & a_{22} & a_{23} \\ a_{31} & a_{32} & a_{33} \end{array} \right) \left( \begin{array}{c} x_{1} \\ x_{2} \\ x_{3} \end{array} \right) = \left( \begin{array}{c} b_{1} \\ b_{2} \\ b_{3} \end{array} \right) \end{equation} 前式の行列は次式に書き下せる。 \begin{eqnarray} \begin{cases} a_{11} x_{1} + a_{12} x_{2} + a_{13} x_{3} = b_{1} & \\ a_{21} x_{1} + a_{22} x_{2} + a_{23} x_{3} = b_{2} & \\ a_{31} x_{1} + a_{32} x_{2} + a_{33} x_{3} = b_{3} & \\ \end{cases} \end{eqnarray} この時、$x_{i}$について次式が導かれる。 \begin{eqnarray} \begin{cases} x_{1}^{(k+1)} = \cfrac{b_{1} - a_{12}x_{2}^{(k)} - a_{12}x_{3}^{(k)}}{a_{11}} & \\ x_{2}^{(k+1)} = \cfrac{b_{2} - a_{21}x_{1}^{(k)} - a_{22}x_{3}^{(k)}}{a_{22}} & \\ x_{3}^{(k+1)} = \cfrac{b_{3} - a_{31}x_{1}^{(k)} - a_{31}x_{3}^{(k)}}{a_{33}} \end{cases} \end{eqnarray} $x_{i}^{(0)}$を初期値として与え、上式に基づき$x_{i}^{(0)},x_{i}^{(1)},x_{i}^{(2)},...$と計算する。 反復回数を$k$とし、$k$を十分大きくすると,$x_{i}^{(k)}$は解に収束し、n元連立一次方程式の場合は次式により計算する。 \begin{equation} x_{i}^{(k+1)} = \frac{b_{i}}{a_{ii}} - \sum_{j=1,j \ne i}^{n} \frac{a_{ij}}{a_{ii}}x_{j} \end{equation} このように解を求める方法がヤコビ反復法である。 解が収束したかどうかは$x_{i}^{(k)}$と$x_{i}^{(k+1)}$の値を比べてその差が許容誤差内に収まっているかどうかで判断する。 ヤコビ反復法では$k$ステップの$x$だけを使って$k+1$の値を求める。 ガウス・サイデル法では$x_{1},x_{2}...x_{n}$と順番に計算した場合、$x_{2}^{k+1}$を求める時には$x_{1}^{k+1}$がすでに計算されており、これらの最新の値を$k+1$の値の計算で使う。

プログラム例

陽解法の場合、CFL条件によりタイムステップ幅に厳しい条件を付ける必要がある。 陰解法の場合はタイムステップ幅により解が不安定になることはない。 ここでは陰解法を用いてガウスサイデル法により一次元の熱伝導方程式を解いたPythonのプログラム例とその結果を以下に示す。
最後の絵はプログラムの出力をプロットしたもので右軸が長さ、上軸が温度を表し、時間ごとに左の線から右の線に変わっている。 図を見ると最初は左端だけ300度であるのが徐々に全体に温度が拡がっていく様子がわかる。

実践的な解法

後退差分による差分近似式を$T_{i,j+1}$について整理すると次式になる。 \begin{equation} T_{i,j+1}=\frac{T_{i,j}+\alpha \gamma (T_{i+1,j+1} + T_{i-1,j+1})}{1+2 \alpha \gamma} \end{equation} 上式を線形反復式として解くこともできる。この場合は行列を意識することがなく簡潔になる。 プログラム例は以下になり、同様の結果を得ることができる。
正直をいうと教科書的なガウスサイデル法を使っている例はほとんどなくて実践的な解法に記載したようなサンプルが多かったので絶賛混乱中や。

2015年10月13日火曜日

クロスプラットフォームの幻想

今現在クロスプラットフォームに対応しました!という場合、Mac,Windows,Unixとかだろう。 モバイルだったらiOS,Androidもあるかもしれない。

この場合、バージョンは一年に一回上がるとして対応するバージョンが5個増えることになる。 その間にも新しいOSが増えるかもしれない。 アーキテクチャが変わるかもしれない。 新しいアーキテクチャが出てくるかもしれない。 呼び出していたAPIがなくなるかもしれない。 ブラウザとかが関わるのであればさらに掛け算して評価が必要だ。 Javaや.NETだったらRuntimeが対応していないかもしれない。 3DTouchのような操作やGUIの見た目が微妙に違うかもしれない。 ライセンス形態が変わるかもしれない。 Windows3.1とか古いOSはもう対応しなくていいかもしれない。 同じものでもOSごとに性能が違うかもしれない。 仮想環境で動かさないといけないかもしれない。 今だとコンテナ技術なんてのもある。 TVとか家電用のOSもあるかもしれない。 使いたい言語をライセンス的に使ったらいけないかもしれない。

要はJavaが言っていたような Write Once Run Anywhere というのは幻想で、 こんなに色んな技術が出てくるとハードウェアやOSに合う方法や言語を使うのが 結局早くなってしまうんでないかなあと思った次第です。

2015年10月9日金曜日

コピペプログラミング

ほとんどこのブログをqiita的に使っていたので素直にqiitaに書くことにした。 qiitaとかstackoverflowとか見て思うのが最近はほとんどググれば答えが載っていること。 訳のわからないバグもエラーコードで大体理由がわかってしまう。 コピペプログラミングはダメだとかいう人もいるけどそもそも意味がわかっていなければコピペで組むこともできないわけで。 ソフトウエアはオブジェクト指向に代表されるように組み合わせるが基本だと思うので、コピペプログラミングも技術の一つと言ってもいいと思う。 実はコピペプログラミングは生産性を10倍ぐらい上げているはずだと思うんだけどそういうことを研究している人はいるんだろうか?

2015年9月24日木曜日

古き良きteapot


3D的なことをしたくなったのでデバッグ用にいろいろOpenGLを調べていたのですが今はOpenGL3.1ごろやらで必須になったCore Profileとかいうのが主流らしいです。 ただデバッグ用に確認したいだけで高速化を求めなければGLUTで十分だなあと思った次第です。 teapotを回転するだけですが時代を感じてよいですね。

2015年9月19日土曜日

C++からCを呼んでみる

かなり使い古されていると思うけどC++からCを呼んでみました。 数値計算や簡単なロジックだけであればC言語だけで十分なのでC++の闇に触れないでC言語で作って他の言語から呼ぶのが、 今から性能を求めるシステムを作る場合はよいのではないかと思います。

2015年9月15日火曜日

iPhoneでOpenGLを使って四角を書く

iPhoneで四角だけを書くものを作りました。 OpenGL ESがわからない上にObjective-Cはしんどいです。 これ正方形を書いたつもりなんだけど長くなっているということは 解像度によって調整しないといけないということなんやろうか。。 ただこの前作った流体シミュレーションと組み合わせると面白いことができそうです。

2015年9月13日日曜日

もくもくする何か

もくもくする何かを作ってみました。 ナビエストークス方程式的な流体をスタガード格子でコンパクトに書いているサンプルがないので作ってみました。 もともと同じようなものがあったのでほとんど佐野ってきています。 時間ができたらiPhoneアプリみたいにしたいですね。

2015年8月16日日曜日

いろんな言語でOOPする

そろそろObjective-Cの勉強もしたくなったのいろんな言語でOOPしてみました。 Javaで書いた1+2だけをするプログラムををC++/Objective-C/Swiftで書いてみます。 これで少しiPhoneアプリのソースも読めそうです。 あとObjective-Cの硬い感じに慣れるとSwiftのチャラっとした感じは受け付けないのも少し理解できました。 ただ結局iPhoneもAndroidもロジックは同じC++が使えるのでC++が結局生き残るだろうなとも思います。 リソース管理はスマートポインタで少し改善できているし、最終的にC++でないとできないことも多いだろうし。

Java

C++

Objective-C

Swift

2015年6月6日土曜日

ポインタの威力

C言語の復習にポインタの威力を改めて勉強してみました。 やはり 3つ以上のポインタは訳がわからなくなりそうですね。

2015年3月10日火曜日

Singleton Collection

いろんな言語でSingletonを実装してみると言語ごとの違いがよく出て面白い。 ただ、C言語はかなりノリで作ったものなので、こうこともできるという参考ということで。

2015年2月11日水曜日

glMatrix2.x系の互換性

WebGLのサンプルを見ているとよくglMatrix0.9を使っていることがあるけど新しいglMatrix2.2を使うと若干引数の位置とかが違ったのでメモ。

2015年1月27日火曜日

明けましておめでとうございます。

明けましておめでとうございます。
WebGLとGithub Pagesの勉強がてら日章旗を描いてみました。

  http://sanofc.github.io/webgl/08/ 

色々めんどくさいことも多いですが、来年ははためかせたいです。。

2015年1月23日金曜日

一人バイキング

楽しかったです。