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で計算した結果がこちら。
でばっちり合ってます。
さて、昨日メルセンヌツイスタで正規乱数を計算できるようになって、本日コレスキー分解ができるようになったということは?(続く)。
----
Excel VBAにNumerical Recipesのルーチンを移植するシリーズ:
三重対角行列
https://sci.tea-nifty.com/blog/2012/07/numerical-recip.html
LU分解
https://sci.tea-nifty.com/blog/2010/07/excel-vbanumeri.html
« メルセンヌツイスタの高速版(SFMT)がExcel VBAで簡単に使える! | トップページ | 三体問題の8の字解をGeoGebra4.0で表示。 »
「学問・資格」カテゴリの記事
- 浜村渚の計算ノート 10さつめ ラ・ラ・ラ・ラマヌジャンを読んだ。九章算術、ベクトル、四元数、電卓、そしてラマヌジャン!タクシー数も1+2+3+…=-1/12もいろんな公式も出てきます。カプレカー数も。高精度計算サイトkeisan.casio.jpにUPしているものとも関連していてよかった。(2023.09.26)
- 出遅れましたがFizzBuzzをExcelのMAP,LAMBDA,SEQUENCE,IFS関数を使って一行(というかただの1セル入力)で作る。(2023.09.24)
- (速報続報)iPhone15 Proが早くも分解。USモデルなのでミリ波アンテナ3つが見えてる。USモデルとその他でMLBの形から違う!(USはeSIM、他はSIMカード)、なのでスペースがなくて他モデルはミリ波アンテナ部分に部品乗せてる!(2023.09.23)
- 離散リアプノフ方程式 AXAᴴ - X + Q = 0がクロネッカー積とvecで計算できることを思い出した!(Vec Trickというそう)。Pythonのscipy.linalg.solve_discrete_lyapunovとnumpy.kronの両方で計算してちゃんと合うことを確認。(2023.09.22)
- 高周波(RF・マイクロ波・ミリ波・5G)関連ニュース2023年9月19日 Microwave Magazineの特集はRFIDや氷を検出する話、Microwave Journalで車載アンテナ評価でRanLOSというのを初めて知る、Gapwavesの多層導波管、PythonのRFライブラリScikit-RFに高木分解を使うTUG multiline TRLが。(2023.09.19)
トラックバック
この記事へのトラックバック一覧です: Numerical Recipes in CのCholesky分解ルーチンをExcel VBAに移植:
» PCを自宅で楽々マスター [楽ぱそDVD]
マイクロソフトのワード・エクセル・パワーポイント・アクセスの操作方法を自宅で楽々マスターできる 動画パソコン教室です。 [続きを読む]
« メルセンヌツイスタの高速版(SFMT)がExcel VBAで簡単に使える! | トップページ | 三体問題の8の字解をGeoGebra4.0で表示。 »
コメント