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

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

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

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

[code lang=”csharp”]
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.[j] - p * F.[i] for k in 0..n do A.[j,k] <- A.[j,k] - p * A.[i,k] // 後退代入 F.[n] <- F.[n]/A.[n,n] for i in n-1..-1..0 do for j in i+1..n do F.[i] <- F.[i] - A.[i,j]*F.[j] F.[i] <- F.[i]/A.[i,i] F [/code]

「有限要素法概説」にあった前進消去「上三角化」したあとに後退代入でひとつずつ計算する。

[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
[/code]

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

[code lang=”csharp”]
val it : Vector = vector [|3.0; 4.0|]
[/code]

こんな風に結果が出る。

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

[code lang=”csharp”]
let M = matrix[[2.;-1.];
[-1.;2.]]
let F = vector[2.0;5.0]

let sol = MathProvider.LinearAlgebra.solve M F
[/code]

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

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