首先... SyntaxHighlighter沒法高亮Mathematica的語法,所以下面只能以這種形式表示了
沒有row interchange的LU分解函數(LU Decomposition without row interchanging):
LUFactor[mat0_] :=
Module[{m, n, i, j, k = 1, L, U = mat0},
{m, n} = Dimensions[U];
L = IdentityMatrix[m];
For[j = 1, j < n + 1, j++,
If[U[[k, j]] != 0,
For[i = k + 1, i < m + 1, i++,
L[[i, k]] = U[[i, j]]/U[[k, j]];
U[[i]] = U[[i]] - U[[k]]*U[[i, j]]/U[[k, j]];
];
k++
]
];
Print["L=", L // MatrixForm];
Print["U=", U // MatrixForm];
Return[{L, U}];
]
將這個放在Mathematica的檔案中執行一次後就能呼叫這個函數,在裡面放入要分解的矩陣就可。Print那兩行是顯示L和U,不需要的話可刪除。如果想將這個函數的L和U指定到另外的變數中,可以寫成{a,b}=LUFactor[matrix],這樣L和U就會指定到a和b這兩個變量上。
事源是這樣的,Mathematica本身內置了一個叫LUDecomposition的函數,但它沒法選擇做不做row interchange,而筆者在溫習線性代數的課題時,在LU分解的練習中有些題目是要在不做row interchange的情況下解的,只好放棄。
然後想說其實筆者也會一點編程,應該查一下Mathematica的編程語法就會了吧,於是就邊查邊找別人寫的函數的例子,就編出上面的function了。
一開始只是寫了4x4的矩陣作試驗,然後擴建成nxn,再把參數改一改,令它能處理mxn,最後加一個If和改改參數變成能處理pivot的問題。
不過Mathematica的資料真難找,搜尋常常夾了一堆不相干的東西,像要找function的定義方法就找了很多數學function的相關東西出來,找programming就變mathematical programming,好麻煩。
2015年12月10日 星期四
訂閱:
張貼留言
(
Atom
)
沒有留言 :
張貼留言