アップデート 2006.9.12

熱方程式
放物型偏微分方程式


u(x , t)を長さLの1次元の棒の熱の分布とする。
放熱の方程式は
ut=c2uxx (0<x<L)
境界条件は
u(0 , t)=u(L , t)=0
初期条件は
u(x , 0)=f(x)
であるとする。
空間軸(x軸)を間隔hで等分し、時間軸(t軸)を間隔kで等分する。これを差分方程式に変換して数値的に解く。その解法は
u(xi , tj+1)=u(xi , tj)+r×{u(xi−1 , tj)−2u(xi , tj)+u(xi+1 , tj)}
ただし
i=ih
tj=jk
r=c2k/h2 (r<0.5が解が安定する条件である。)
である。


 EXCELのシートは次のようになる。
 緑色で示したB2からB8までのセルに必要なデータを入力する。なおB2からB8までのセルには名前を定義しなければならない。それぞれの名前は左隣のA2からA8までのセルに与えている。C2からC8までのセルに入力したデータの説明を書いておく。このようにしないと後で何を与えたのか忘れてしまう。
 白色で表示されたセルB10にはrの値をEXCELのマクロで計算して与えている。この値が0.5を越えると数値解が収束することが保証されない。マクロではE1を先頭にして数値解を与えるようにしている。セルB11に空間の分割数を表示している。これによりE列からY列のセルに数値解が表示される。セルB12に時間の分割数をマクロで計算して表示している。これにより2行から102行までのセルに数値解が表示される。セルB13に数値解が表示されるセルの数を表示する。この値が大きいと計算に時間がかかる。
 水色で示したE1からY1までのセルの数値が空間座標Xの値である。E2からY2までのセルの数値が時刻t=0での数値解、すなわち初期関数f(x)の値である。初期関数f(x)はこのシートには表示されていない。この関数はマクロを見ないと分からない。マクロではf(x)=x2(2-x)がこの例題で与えた初期関数である
 黄色で表示されたD2からD102までのセルの数値が時刻である。
 灰色で表示されたE2からY102までの数値が数値解である。

A

B

C

D

E

F

G

H
1 熱方程式


0

0.1

0.2

0.3

2 coeff

1

伝導率

0

0

0.019

0.072

0.153

3 length

2

棒の長さ

0.005

0

0.036

0.086

0.164

4 step_x

0.1

空間刻み幅

0.01

0

0.043

0.1

0.175

5 step_t

0.005

時間刻み幅

0.015

0

0.05

0.109

0.186

6 u_l_value

0

棒の左端の温度

0.02

0

0.045

0.118

0.1945

7 u_r_value

0

棒の右端の温度

0.025

0

0.059

.01245

0.203

8 time_end

0.5

最終時間

0.03

0

0.06225

0.131

0.209625

9

0.035

0

0.0655

0.1359375

0.209625

10 ratio

0.5

安定条件

0.04

0

0.06797

0.140875

0.21625

11

20

空間分割数

0.045

0

0.07044

0.1447188

0.2214688

12

100

時間分割数

0.05

0

0.07236

0.1515859

0.2266875

13

2000

表示セル数

0.055

0

0.07428

0.1564094

0.2349375

14


0.060

0

0.07579

0.546094

0.2349375


 この例題の場合、E2からY102までの範囲のセルに計算されたu(x , t)が記録される。D1からY102までの範囲選び、グラフウイザードで等高線を選択すると下図が描かれる。メモリが0から2までと表示された左下の辺が空間軸である。メモリが0から0.5までと表示された右下の辺が時間軸である。上方に伸びたメモリのついた軸は温度分布である。時刻0で最も高い曲線になっている温度分布は時間と共に指数関数的に減少していることが分かる。


 空間軸を横とし、温度を縦軸として散布図を描くと、温度分布の変化が分かる。最も高い非対称の曲線はt=0のときの初期温度分布でf(x)=x2(2-x)である。温度は順次棒の両端から失われて下がっていく様子が観察される。xが小さい場所では、棒の左端から失われる熱よりも棒の右方向から供給される熱の方が大きく、当初の温度が高くなっていることが観察される。



放物型偏微分方程式を解くマクロ
問題1:時間と温度との関係を調べよ。
問題2:f(x)の形を変更して、解を図示せよ。ただしf(0)=f(L)=0でなくてはならない。
問題3:r<0.5でないと安定して収束しないと言うことは、どのようなことか説明せよ。
参考文献:足立紀彦・竹内敬治、メカトロのための数値計算、森北出版、1988