« メルセンヌツイスタの高速版(SFMT)がExcel VBAで簡単に使える! | トップページ | 三体問題の8の字解をGeoGebra4.0で表示。 »

2012年8月31日 (金)

Numerical Recipes in CのCholesky分解ルーチンをExcel VBAに移植

ニューメリカルレシピをExcel VBAに移植するシリーズ。今回はコレスキー分解。

これは非常に簡単に、

Sub choldc(a() As Double, n As Long, p() As Double)
     Dim i As Long, j As Long, k As Long
     Dim sum As Double
    
     For i = 1 To n
        For j = i To n
            sum = a(i, j)
            For k = i - 1 To 1 Step -1
                sum = sum - a(i, k) * a(j, k)
            Next k
            If (i = j) Then
                If (sum <= 0#) Then
                    MsgBox "choldc failed"
                End If
                p(i) = Sqr(sum)
            Else
                a(j, i) = sum / p(i)
            End If
        Next j
    Next i

End Sub

と書ける。

使い方は、aを2次元配列で行列を表し、正の値を持つ対称行列とする。で上三角行列だけを入力しておけばOKで、下三角行列にコレスキー分解した値が入る。ただし、対角成分は一次元配列のpに入る。

例えば5x5のHilbert行列 ( aij = 1/(i+j-1) )をコレスキー分解するには、

Option Explicit

Private Sub CommandButton1_Click()
    Dim n As Long
    Dim a() As Double, p() As Double, b() As Double
    Dim i As Long, j As Long
   
    n = 5
    ReDim a(n, n)
    ReDim b(n, n)
    ReDim p(n)
   

   
    For i = 1 To n
        For j = i To n
            a(i, j) = 1# / (i + j - 1)
            Worksheets("sheet1").Cells(i + 1, j + 1) = a(i, j)
        Next j
    Next i

    Call choldc(a, n, p)
   
   
   
    For j = 1 To n
        For i = j To n
            Worksheets("sheet1").Cells(i + 1, j + 7) = a(i, j)
        Next i
    Next j
    For i = 1 To n

            Worksheets("sheet1").Cells(i + 1, n + 9) = p(i)

    Next i
End Sub

で呼び出す。Maximaで計算した結果がこちら。

Hilbert

でばっちり合ってます。

さて、昨日メルセンヌツイスタで正規乱数を計算できるようになって、本日コレスキー分解ができるようになったということは?(続く)。

----

Excel VBAにNumerical Recipesのルーチンを移植するシリーズ:

三重対角行列

http://sci.tea-nifty.com/blog/2012/07/numerical-recip.html

LU分解

http://sci.tea-nifty.com/blog/2010/07/excel-vbanumeri.html

« メルセンヌツイスタの高速版(SFMT)がExcel VBAで簡単に使える! | トップページ | 三体問題の8の字解をGeoGebra4.0で表示。 »

学問・資格」カテゴリの記事

コメント

コメントを書く

(ウェブ上には掲載しません)

トラックバック

この記事のトラックバックURL:
http://app.cocolog-nifty.com/t/trackback/512682/55542902

この記事へのトラックバック一覧です: Numerical Recipes in CのCholesky分解ルーチンをExcel VBAに移植:

» PCを自宅で楽々マスター [楽ぱそDVD]
マイクロソフトのワード・エクセル・パワーポイント・アクセスの操作方法を自宅で楽々マスターできる 動画パソコン教室です。 [続きを読む]

« メルセンヌツイスタの高速版(SFMT)がExcel VBAで簡単に使える! | トップページ | 三体問題の8の字解をGeoGebra4.0で表示。 »

最近の記事

最近のコメント

2017年12月
          1 2
3 4 5 6 7 8 9
10 11 12 13 14 15 16
17 18 19 20 21 22 23
24 25 26 27 28 29 30
31            
フォト
無料ブログはココログ