正向替換包裝數組

所以我試圖創建一個腳本，將採取一個數組，將其轉換成壓縮形式，然後對它運行Cholesky因式分解，最後使用正向和反向替換來解決它（Ax=b）。正向替換包裝數組

我現在遇到的麻煩是爲打包數組創建正向替換（列明智）算法。我已經正確地完成了後向替換，但是我似乎無法正確轉換算法。

這是當前的後向代入算法我有（這是完美的）

x = b; 
p = n*(n+1)/2; 

for j = n:-1:1 
    x(j) = x(j)/u(p); 
    p = p - 1; 
    for i = j-1:-1:1 
    x(i) = x(i) - u(p) * x(j); 
    p = p - 1; 
    end 
end

凡n是原始A矩陣（其是正方形和SPD）的大小。從A矩陣u的映射如下：

A(i,j) = u(i+j*(j-1)/2)

我現在正向替換算法的迭代如下：

x = b; 
p = 1; 

for j = 1:n 
    x(j) = x(j)/u(p); 
    p = p + 1; 
    for i = j+1:n 
     x(i) = x(i) - u(p) * x(j); 
     p = p + 1; 
    end 
end

我似乎無法找出我做錯了什麼。無論我嘗試結束，只是給我NaN或Inf作爲x的答案。比我聰明的人能算出算法應該是什麼樣子嗎？

來源

2016-10-19 Bradley Gray

從我所瞭解的您試圖解決的兩個問題導致Ax = b的解決方案。

哪裏A = LL'也就是說A是首先用着替代解決Ly = b等於乘以其移調（L'）

下三角矩陣（L）。

最後通過解決L'x = y與回代換。

凡你的背部substituion你u(p)被定義爲這樣的（對於n = 3）：

 | u(1) u(2) u(4) | 
L' = | 0 u(3) u(5) | 
    | 0 0 u(6) |

我認爲你正在運行到與前substition問題來形成什麼L將被定義爲：

| u(1) 0 0 | 
L = | u(2) u(3) 0 | 
    | u(4) u(5) u(6) |

你的價值p需要有以下進展來完成你想與前substition什麼：

{1, 2, 4, 3, 5, 6}

我已經修改了您當前的轉發代理代碼來實現此模式。

x = b; 
p = 1; 

for j = 1:n 
    d = j-1; 
    p = j + (d^2+d)/2; 
    x(j) = x(j)/u(p); 
    d = d + 1; 
    for i = j+1:n 
     p = j + (d^2+d)/2; 
     x(i) = x(i) - u(p) * x(j); 
     d = d + 1; 
    end 
end

我能夠測試這個，並且它產生了預期的結果。我希望這可以幫助你。

來源

2016-10-19 15:06:54

正向替換包裝數組

回答

相關問題