[F#] 連立一次方程式を解く

逆行列、行列式とできたので、連立一次方程式を解かせる。有限要素法のいわゆる「solve」というやつで、解析の手順として、

• pre : 構造を設定、要素に分割
• solve : 連立一次方程式を解く（それだけじゃない？）
• post : ひずみの結果を表示する

という手順になっている。
Gauss の消去法を使って直接的に解くわけだが、これって、そのまんま LU 分解なわけで。なので、solve 自体は反復法なども含めた意味で使わないとダメっぽい。

let solve (A' : matrix) (F': vector) = 
    let A = A'.Copy()
    let F = F'.Copy()
    let n = A.NumCols-1
    // 前進消去
    for i in 0..n-1 do
        for j in i+1..n do
            let p = A.[j,i]/A.[i,i]
            F.[j] <- F.&#91;j&#93; - p * F.&#91;i&#93;
            for k in 0..n do
                A.&#91;j,k&#93; <- A.&#91;j,k&#93; - p * A.&#91;i,k&#93;
    // 後退代入
    F.&#91;n&#93; <- F.&#91;n&#93;/A.&#91;n,n&#93;
    for i in n-1..-1..0 do
        for j in i+1..n do
            F.&#91;i&#93; <- F.&#91;i&#93; - A.&#91;i,j&#93;*F.&#91;j&#93;
        F.&#91;i&#93; <- F.&#91;i&#93;/A.&#91;i,i&#93;
    F
&#91;/code&#93;
<p>
「有限要素法概説」にあった前進消去「上三角化」したあとに後退代入でひとつずつ計算する。
</p>
[code lang="csharp"]
// 2x -y = 2
// -x +2y= 5 を解く

let M = matrix[[2.;-1.];
               [-1.;2.]]
let F = vector[2.0;5.0]

solve M F 

こんな風にいれておいて、実行すると

val it : Vector<float> = vector [|3.0; 4.0|]

こんな風に結果が出る。

実は、F# MathProvider にも同じものがあって、

let M = matrix[[2.;-1.];
               [-1.;2.]]
let F = vector[2.0;5.0]

let sol = MathProvider.LinearAlgebra.solve M F

として解くと同じ事ができる。

有限要素法の場合、ひとつの節点の自由度は6になるので、節点の数の6倍の計算量になる。さらに、要素が持つ節点（頂点）は三角形二次要素場合は6点になるので、要素数の36倍の計算量なる。が、計算量が多くなるといっても、多くなっても2桁増えるぐらいだから大したことはない？ことはないか。前進消去で、n^3/2 ぐらいのオーダーになるので、n が 100倍になると、100万倍ぐらいの計算量が増える。となると、この素直な方法で解いてしまうと、メッシュ数を2倍細かくすると約8倍の計算量が必要になる。確かに、高速化が必要な分野であるかも。
このあたりは後で実測して試してみよう。