TIM Labs

2016年9月アーカイブ

全体の流れを図にしてみた。
図中の赤枠は、格子点1つに対する偏微分方程式を示している。

  1. まず、サイズNの正方行列Bの外周に対して温度を決める。
  2. Bを列ベクトルbに整形する。
  3. サイズN*Nの正方行列(係数行列)を作る。
  4. 係数行列の逆行列を両辺に左から掛けることで、温度列ベクトルtが求まる。
  5. tをサイズNの正方行列Tに整形すると、全格子点の温度が求まる。

poisson2d-8b.png
上図や今までに説明したことを元にプログラムを書くと、以下のようになる。
 
# 2元定常熱伝導の偏微分方程式を解こう

# 条件設定
TH = 1600;
T0 = 300;
N = 50; NN = N*N;
L = 0.1;
x = linspace(0,L,N);

# 境界条件の設定
B = zeros(N,N);
B(1,:)= B(N,:)= B(:,N) = T0;
B(1:N,1) = sin(pi*x/L)*(TH-T0) + T0;
b = reshape(B,NN,1);

# 係数行列の作成
C = full(gallery("poisson",N));
I = eye(NN);
C(1:N,:) = I(1:N,:);
C(1:N:NN,:) = I(1:N:NN,:);
C(N:N:NN,:) = I(N:N:NN,:);
C(NN-N+1:NN,:) = I(NN-N+1:NN,:);
C = sparse(C);

# 微分方程式を解く
t = C\b;
T = reshape(t,N,N);

# 温度分布を図示
[X,Y] = meshgrid(x,x); 
surf(X,Y,T);
view(50,30);
poisson2d-2.png
偏微分方程式の一番簡単な例を示してきたが、考え方が分かれば3次元化も可能であるが、格子点の数が劇的に多くなり、そのままでは殆んど計算できなくなってしまう。
対象領域が正方領域で、x、y両方をを同じ間隔で区切ったが、異なる間隔で区切ることも、対象領域が正方形ではなく長方形であっても、係数行列を工夫すれば計算可能になるのも分かるであろう。
さらに、対象領域が長方形でなく、任意な閉領域で、領域内の格子点に対して偏微分方程式が成り立つ場合も、面倒であるが同様にして計算できることが分かるであろう。
これらは、ここに書くのは面倒なので省略する。

次回は、別のことを書こうと思う。

偏微分方程式を解くために、まず、行列のreshapeを理解しないといけない。
評価したい対象領域をN*Nに区切ったので、対象領域はN*Nの行列(2次元の配列)になっている。実際には温度のN*N配列である。

しかし、温度が行列のままだと、C * T = B という感じで、偏微分を差分表現したのを行列演算に持ち込めない。
なお、Cは係数行列、Tは温度、Bは差分表現の右辺あるいは境界条件である。

そこで、行列を列ベクトルにすることを考える。

行列を整形する関数は非常にたくさんあるのだが、今回使うのはreshape()である。

>> A = 1:9
A =

   1   2   3   4   5   6   7   8   9

>> A = reshape(A,3,3)
A =

   1   4   7
   2   5   8
   3   6   9

>> A = reshape(A,9,1)
A =

   1
   2
   3
   4
   5
   6
   7
   8
   9

>>
したがって、 C * T = B において、元々2次元配列であったTとBを1次元列ベクトルにできるはずだ。

求めたいのは温度Tなので、C * T = B の両辺にCの逆行列をかけると次のようになる。

T = C \ B

差分表現の右辺あるいは境界条件をN*Nの行列として作る。 外周以外は、偏微分方程式が成り立ち右辺が0なので、zeros()で全部0のN*Nの配列を作る。 外周は温度を外気温T0=300を指定するのだが、過熱している辺についてだけ、sinカーブになるようにする。 そして、最後に、N*Nの列ベクトルbを作っている。
B = zeros(N,N);
B(1,:)= B(N,:)= B(:,N) = T0;
B(1:N,1) = sin(pi*x/L)*(TH-T0) + T0;
b = reshape(B,NN,1);

これに合わせて、前回 gallery("poisson",N) で作った行列うち、上記の行列Bで温度を指定した箇所に対応する行を、単位行列の行と同じにするために、単位行列の対応行をコピーするという強引な方法をやっている。

C = full(gallery("poisson",N));
I = eye(NN);
C(1:N,:) = I(1:N,:);
C(1:N:NN,:) = I(1:N:NN,:);
C(N:N:NN,:) = I(N:N:NN,:);
C(NN-N+1:NN,:) = I(NN-N+1:NN,:);
C = sparse(C);

N=4の場合の係数行列は次のようになる。
前回と比較すると、対角線上に1が並んでいる箇所が元の領域の外周に相当する行を示す。
C =

   1   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0
   0   1   0   0   0   0   0   0   0   0   0   0   0   0   0   0
   0   0   1   0   0   0   0   0   0   0   0   0   0   0   0   0
   0   0   0   1   0   0   0   0   0   0   0   0   0   0   0   0
   0   0   0   0   1   0   0   0   0   0   0   0   0   0   0   0
   0  -1   0   0  -1   4  -1   0   0  -1   0   0   0   0   0   0
   0   0  -1   0   0  -1   4  -1   0   0  -1   0   0   0   0   0
   0   0   0   0   0   0   0   1   0   0   0   0   0   0   0   0
   0   0   0   0   0   0   0   0   1   0   0   0   0   0   0   0
   0   0   0   0   0  -1   0   0  -1   4  -1   0   0  -1   0   0
   0   0   0   0   0   0  -1   0   0  -1   4  -1   0   0  -1   0
   0   0   0   0   0   0   0   0   0   0   0   1   0   0   0   0
   0   0   0   0   0   0   0   0   0   0   0   0   1   0   0   0
   0   0   0   0   0   0   0   0   0   0   0   0   0   1   0   0
   0   0   0   0   0   0   0   0   0   0   0   0   0   0   1   0
   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   1
 

次の図で見たほうが分かりやすいだろう。
poisson2d-6.png

以上で、一通りの説明を終えたので、次回はプログラムを示す。
前回、偏微分方程式の数値解の結果を示したが、これをどうやって計算したかについて説明していこうと思う。
まず、対象とする一辺Lの正方形領域をメッシュに切ろう。頂点がN*N個になるように切って計算したのが前回の図である。すると、格子の一辺の長さΔはL/(N-1)なのだが、そんなものは全体に共通なので係数を掛けたり割ったりだけに過ぎず本質的でないので、格子の一辺の長さは1としよう。これで計算が楽になる。

次に、2次のポアソンの偏微分はどうすれば計算できるかを考える。微分、偏微分では、隣の点との差分を計算すればよく、2階編微分になっているので、差分の差分をとれば大丈夫ということで、次図のイメージで計算すればよいと思われる。

poisson2d-4.pngところで、上図の各格子点に対して、上下左右の隣点の値(温度)の和から対象点(i,j)の値(温度)の4倍を引けば2階編微分に相当する値が出るはずである。今回は、この値が0で、方程式が格子点(i,j)に対して1つ出来上がる。
そして、こういう点がn*n個あるので、方程式の数がn*n個と膨大になる。全ての格子点の温度T(i,j)が変数なので、変数の数がn*n個。そして、方程式の数がn*n個だと、これはn*n元の連立一次方程式になり解けるはずなのだ。

係数行列をCとすると、Cはサイズn*nの正方行列になるわけで、格子点(i,j)に対する式は、i+(j-1)*n 行目の係数で表される。この点については、どういう風に計算するかの説明が必要なのだが、境界条件の説明するときに合わせて行う。


この係数行列は規則的なので作ろうと思えば作れるのだが、Octaveは非常に多数の便利な行列を自動生成してくれる機能があり、そのための関数がgallery()である。

>> full(gallery("poisson",4))
ans =

   4  -1   0   0  -1   0   0   0   0   0   0   0   0   0   0   0
  -1   4  -1   0   0  -1   0   0   0   0   0   0   0   0   0   0
   0  -1   4  -1   0   0  -1   0   0   0   0   0   0   0   0   0
   0   0  -1   4   0   0   0  -1   0   0   0   0   0   0   0   0
  -1   0   0   0   4  -1   0   0  -1   0   0   0   0   0   0   0
   0  -1   0   0  -1   4  -1   0   0  -1   0   0   0   0   0   0
   0   0  -1   0   0  -1   4  -1   0   0  -1   0   0   0   0   0
   0   0   0  -1   0   0  -1   4   0   0   0  -1   0   0   0   0
   0   0   0   0  -1   0   0   0   4  -1   0   0  -1   0   0   0
   0   0   0   0   0  -1   0   0  -1   4  -1   0   0  -1   0   0
   0   0   0   0   0   0  -1   0   0  -1   4  -1   0   0  -1   0
   0   0   0   0   0   0   0  -1   0   0  -1   4   0   0   0  -1
   0   0   0   0   0   0   0   0  -1   0   0   0   4  -1   0   0
   0   0   0   0   0   0   0   0   0  -1   0   0  -1   4  -1   0
   0   0   0   0   0   0   0   0   0   0  -1   0   0  -1   4  -1
   0   0   0   0   0   0   0   0   0   0   0  -1   0   0  -1   4


第1引数に、行列の種類を文字列で指定する。"poisson"のとき、第2引数に元の分割数n(4)を指定すると、サイズがn*n(16)の正方行列が作られる。
full()は、gallary()が出力する行列がスパースな行列なので、そのまま表示すると、0以外が入っている箇所しか示さないので、普通の行列に展開するために使った。よく使われる行列の多くは、ほとんどの要素が0というのが多く、値の入っているところだけデータを持つことで、メモリが非常に圧縮できるし、計算も早くなったりする。

最初の図に従うと、点(i,j)の係数は-4で、隣の点は1でなければいけないのだが、この行列は符号が逆転している。でも、どうせ偏微分の結果は0なのだから、符合が逆転しても結果は同じになるので、このまま使ってしまおう。


配列は大きくなると、一覧性がどんどん悪くなる。それで、配列の状態を一瞬で把握できるような関数も用意されている。

>> colormap("gray")
>> imagesc(full(gallery("poisson",6)))
poisson2d-5.png
6x6の格子点に対するポアソンの行列なので、36次の正方行列になった。
対角線の白マスが4を示し、グレイが0、黒が-1になっている。
上下左右の4点を指定しようとしても、外周上の格子点の場合、4つの隣点のうちの少なくとも1つは存在しない点になるため、横の行を見ると、黒が本来4箇所ないといけないのだが、3箇所や2箇所になっている場所がある。これは境界に対して発生しているわけである。

さて、境界の問題が出てきたが、そもそもポアソンの方程式、解こうと思ったら何らかの条件が必要になる。今回は、外周上の格子点の温度を指定することで解く。
このような境界条件を指定しなければいけないのだが、それに合わせて、得られた行列も調整する必要があるのだが、それは次に説明視よう。

ここでは、非常に足早に必要なことも省略しながら説明しているので、この分野の計算についてちゃんと知りたい場合は、「偏微分方程式の数値解法」などのキーワードでググってみよう。あちこちの大学の講義資料や書籍がみつかるので、それらできちんと勉強しよう。

せっかくOctaveについて書いているのに、Octaveが得意な科学技術計算についてちょっとしか書いていない。
ということで、今回は、とつぜん偏微分方程式の数値解に挑戦してみようと思う。

物理現象を偏微分方程式で表すことはよくあるのだが、偏微分方程式をきちんと解くのはほとんど不可能である。
なので、現実には、適当な境界条件のもとで数値解を求めることになる。

とりあえず、結果を先に示そう。
2次元の熱伝導について数値解を求めてみよう。

1辺0.1m(10cm)の板の1つの辺の中央に熱が加わり、周囲の温度が下図のようになる場合について計算してみよう。
温度は、外周で300度(絶対温度、室温)、一辺の中央が1600度で、その一辺の温度分布はsinカーブになっているものとする。
高さ方向は温度になっている。

poisson2d-1.png外周は上図の温度であるが、内部の温度は次の偏微分方程式にしたがって変化する。
右辺が0になっているのは定常状態を示している。
というより、0が計算が楽なのでこれで今回は説明をしてしまう。

poisson2d-3.png
これを、正方形の板をN*Nのメッシュに切って、各点の温度を求めて図示するとこんな感じになる。
高温が赤、低温が青になっているが、これはOctaveで勝手に色づけしてくれるので便利。

poisson2d-2.png
こんな図になるように計算するにはどうすれば良いのだろうか。
この図は、N=50の場合の計算結果である。

実際には、物体は3次元空間に存在するのだから、x,y だけでなく、zも加えて、N*N*Nのメッシュに切って計算するべきなのだが、こうなると点の数が非常に増えて困ることになる。
点の数は、2次元の場合で N*N、3次元の場合に N*N*N になるから大変なのだ。
いや、本当は次元が増えるというのはそんなに甘いものではない。もっと急激にデータ爆発が起きるのだが、そのあたりは次で説明しよう。

RAII の話題は C++ の例が多いが、それもそのはず、ここ最近は GC が面倒を見てくれることがほとんどなので、 RAII の需要そのものが減っている。とはいえネットワークコネクションやらファイルハンドルやらと、あまりリソース解放のタイミングを遅らせたくないものに関しては、 RAII パターンで寿命を明示しておきたいのは GC があったとしても同じこと。

Ruby もそういったリソースに関してはブロックを渡すことで解放もしてくれる RAII パターンを実装したインターフェースが提供されているので、積極的に使っていきたいところだ。そういったリソースがひとつだけなら、素直に書けばいいだけのことなのであまり考える必要はない。

open(...) do |fp|
  something(fp)
end # ここでファイルハンドルが解放される

ところが、配列になっているものに対し、すべての要素を開き、同時に使いたいとなると意外と面倒なことに気付く。配列 arr = [a1, a2, a3, ...] があったとき、意味的にはこういうことがしたい。

open(a1) do |b1|
  open(a2) do |b2|
    open(a3) do |b3|
      ...
      something([b1, b2, b3, ...])
    end
  end
end

が、配列の要素数がわからない以上、この方針でネストを続けていくわけにはいかない。 果て...。どうすればいいだろう?

よくあるモンテカルロ法による円周率を求めるというのをOctaveでやってみる。
図のように、x,y それぞれ-1〜1の範囲の乱数を生成し、一辺2の正方形と単位円を重ねて、単位円の中に入っている点の割合から円周率を求める。
pimontecarlo.png
まずは、ごく普通のプログラミング言語風の書き方で試す。
一応、時間を計測するために、time()を使っている。

function [p, t] = pislow( n )
    t1 = time();
    cnt = 0;
    for i=1:n
        x = rand()*2-1;
        y = rand()*2-1;
        if( x^2+y^2 < 1 )
            ++cnt;
        endif
    endfor
  
    p = cnt / n * 4;
    t = time() - t1;
endfunction

実行すると、こんな感じになる。 精度を上げるために100万回試したら、27秒もかかった。 プログラミング言語としてはかなり遅い方だろう。

>> [p,t] = pislow(1000)
p =  3.1560
t =  0.027001
>> [p,t] = pislow(1000000)
p =  3.1434
t =  27.356

では、Octave的なプログラムにするとどうなるか、見てみよう。
 
function [p,t] = pifast( n )
    t1 = time();
    x = rand(n,1)*2-1;
    y = rand(n,1)*2-1;
    cnt = length(find(x.^2+y.^2<1));
    
    p = cnt / n * 4;
    t = time() - t1;
endfunction
forループが消え、100万個の乱数の配列x,yを作って、配列同士の演算をしている。 findは、条件を満たす場合の添え字の配列を返すので、その個数をlength()で数えることで、条件を満たす点数をカウントしている。

これを実行すると、100万回の場合、340倍ほど高速になっている。 これだけの差が出るので、forループを捨てて、配列演算をどんどん使う重要性が分かるだろう。
つまり、forループなどを使わず、配列演算を利用して短く高速に、バグが入る余地も少なくして横着に高速なプログラムを書くのがOctave風なのだ。

>> [p,t] = pifast(1000)
p =  3.1160
t =  0.0010002
>> [p,t] = pifast(1000000)
p =  3.1421
t =  0.080005

でも、注意すべきことがある。配列演算はとても高速なのだが、配列を確保するので、メモリをどんどん消費する。 要素数100万個の配列で、1要素が8バイトだと、1つの配列で8MB消費する。 1000万だと80MB、1億だと800MB、10億だと8GBにもなってしまい、実際に走らすとシステムが反応しなくなる。

最後に、最初に示した図を描くプログラム(スクリプト)を示そう。
 
hold on;

n = 1000;
x = rand(n,1)*2-1;
y = rand(n,1)*2-1;
lt = find(x.^2+y.^2<1);
ge = find(x.^2+y.^2>=1);

plot(x(lt),y(lt),".r");
plot(x(ge),y(ge),".b");

th = 0:pi/20:2*pi;
plot(cos(th),sin(th),"linewidth",2,"color","red");

hold off
hold on; hold off; があるが、これはFigureの重ね描きを制御するためである。

findで、引数の条件を満たす添え字の並びが得られるが、これを元の配列の引数として与えると、条件を満たす要素(値)だけ抽出することができる。これを利用して、条件を満たす場合と満たさない場合をplotを2回使って色分けして表示している。
Octaveで行列などが扱いやすく数値計算に向いていることは分かったが、結果の数字だけを眺めても様子が分からない。
計算結果を図示することがとても重要だ。図示することで、一瞬でデータの意味が分かる。

とりあえず、x-y 平面の原点からの距離に応じて高さを決めることで、波を描いてみた。

>> a = linspace(-20,20,131);
>> [x,y] = meshgrid(a,a);
>> r = sqrt( x.^2 + y.^2) + eps;
>> z = sin(r) ./ r;
>> surface(x,y,z);
wave.png
aは、-20から20までを130等分した配列である。131になっているのは、両端を含む等分点全部となっているため。
meshgridにより、xにはx座標値が入った配列、yにはy座標値が入った配列が作られる。
rには、原点から座標(x,y)までの距離が入るが、epsを加えるというオマジナイをしている。
そして、波の高さをzとして計算している。
面の表示はsurfaceで、メッシュのx座標配列、y座標配列、z座標配列を引数に与えると、上図が現れる。
マウスを動かすことで、好きな方向から眺められる。

次に、もうちょっと複雑な形を作ってみよう。でも、同じくたった5行でできる形だ。

>> [u,v] = meshgrid(linspace(0,10*pi,200),linspace(0,2*pi,30));
>> x = 2*(1-exp(u/(6*pi))).*cos(u).*cos(v/2).^2;
>> y = 2*(-1+exp(u/(6*pi))).*sin(u).*cos(v/2).^2;
>> z = 1-exp(u/(3*pi))-sin(v)+exp(u/(6*pi)).*sin(v);
>> surface(x,y,z);
snail.png巻貝のような形になっているが、[u,v] の媒介変数で曲面を表しただけである。
uは0から10π、vは0から2πのパラメトリック曲面のx,y,z座標値を求めて表示しているだけなので、処理の原理は同じである。
x,y,zの式をいろいろ弄れば、色々な3次元曲面が作れる訳である。

このアーカイブについて

このページには、2016年9月に書かれたブログ記事が新しい順に公開されています。

前のアーカイブは2016年8月です。

次のアーカイブは2016年10月です。

最近のコンテンツはインデックスページで見られます。過去に書かれたものはアーカイブのページで見られます。