F# で熱伝導の数値解析を頑張ってみた

F# Advent Calendarの18日目です。

F# Advent Calendar 2012 : ATND

さて、F#で何かしてみよう！

ということで、

今回は熱伝導の数値解析（のようなも）をしてみました！

＜はじめに＞

数値解析をするにあたり、できるだけ簡素化して考えてみました。

ざっくりとした仮定はこんな感じ。

１）対象の範囲は20 * 20の二次元配列とする

２）対象の範囲の３辺は壁、残る１辺を熱源とする

３）壁と熱源以外の要素の値はその要素の上下左右の要素の平均とする

＜プログラムについて＞

文字にすると仮定がわかりにくいですが、プログラムの流れとしては下記の通り

１）初期値の配列を生成する

２）関数に配列を渡し、仮定３）に基づき計算の上新しい配列を返す

３）２）の繰り返し

実際のソースの一部を以下に示します。

■初期値の作成（※ 熱源は９、壁、その他の要素は０とします。）

//----------------------------------------------------------
// 定数の定義
//
let width  = 20
let height = 20

//----------------------------------------------------------
// 熱源・壁を定義するメソッド
//
let heatSource (tempHeatLevel : int) (tempWidth : int) = Array.create tempWidth tempHeatLevel

let heatSource_heatLevelIs_0 (tempWidth : int) = heatSource 0 tempWidth
let heatSource_heatLevelIs_9 (tempWidth : int) = heatSource 9 tempWidth

//----------------------------------------------------------
// 配列の初期化メソッド
//
let rec initializeArray (tempWidth : int) (tempHeight : int) =
  if tempHeight = 1
    then
      heatSource_heatLevelIs_9 tempWidth
  else
      Array.append (initializeArray tempWidth (tempHeight - 1)) (heatSource_heatLevelIs_0 tempWidth)

■引数の配列から新しい配列を生成して返す（メインの処理）

//----------------------------------------------------------
// メインの演算メソッド（tempIndexには(array.Length - 1) - width * 2 を指定）
//
let rec getResultArray (tempWidth : int) (tempIndex : int) (tempArray : int[]) =
  if tempIndex = 0 
    then
      [|0|]
  elif ((tempIndex % tempWidth) = 0) || ((tempIndex % tempWidth) = (tempWidth - 1))
    then
      Array.append
        (getResultArray tempWidth (tempIndex - 1) tempArray)
        [|0|]
  else
    if ((double tempArray.[tempIndex + tempWidth + 1] +
         double tempArray.[tempIndex + tempWidth - 1] +
         double tempArray.[tempIndex + tempWidth * 2] +
         double tempArray.[tempIndex]                 +
         2.0) / 4.0 * rateOfHeatCondactor) < 9.0
      then
        Array.append
          (getResultArray tempWidth (tempIndex - 1) tempArray)
          [|int (Operators.round (double tempArray.[tempIndex + tempWidth + 1] +
                                  double tempArray.[tempIndex + tempWidth - 1] +
                                  double tempArray.[tempIndex + tempWidth * 2] +
                                  double tempArray.[tempIndex]                 +
                                  2.0) / 4.0 * rateOfHeatCondactor)|]
    else
        Array.append
          (getResultArray tempWidth (tempIndex - 1) tempArray)
          [|9|]

//----------------------------------------------------------
// メインの演算結果に境界条件を追加するメソッド
//
let getResultArrayWithBoundaryCondition (tempFiestRow : int[]) (tempLastRow : int[]) (tempArray : int[]) =
  Array.append
    tempFiestRow
    (Array.append
       tempArray
       tempLastRow)

let getResultArrayWithBoundaryCondition_firstRowIs_9_lastRowIs_0 (tempArray : int[]) =
  getResultArrayWithBoundaryCondition (heatSource_heatLevelIs_9 width) (heatSource_heatLevelIs_0 width) tempArray

■メインの処理の繰り返し

let rec getResult (tempLoopCount : int) (tempArray : int[]) =
  if tempLoopCount = 0 
    then
      getResultArrayWithBoundaryCondition_firstRowIs_9_lastRowIs_0 (getResultArray width (tempArray.Length - 1 - width * 2) tempArray)
  else
      getResultArrayWithBoundaryCondition_firstRowIs_9_lastRowIs_0 (getResultArray width (tempArray.Length - 1 - width * 2) (getResult (tempLoopCount - 1) tempArray))

■このほかにも熱伝導率的なパラメータなんかも付け加えています。

＜結果＞

処理の結果をまとめるとこんな感じになります。

■５回処理後