c++矩陣求逆的lu分解實現
初學c++,嘗試利用基礎知識編寫矩陣求逆。但發現求解伴隨陣的過程非常複雜且難以實現,而我正好看到一篇求三角陣伴隨矩陣的文章,故嘗試程式設計實現。在這種方法下,計算量明顯減小,實現方法,思路適合初學者。
參考文獻:三角形矩陣求伴隨矩陣的一種方法(曾月新)
求逆矩陣思路:
1.求矩陣的crout(lu)分解,其中l為下三角矩陣,u為上三角矩陣
2.求l,u矩陣的伴隨陣,見參考文獻
3.求l,u矩陣的逆(伴隨陣a* /det(a))
4.inv_a = inv_u * inv_l
附上**如下:
#include
#include
#include
#include
using namespace std;
const int n=3;
//逆矩陣階數,求不同階數矩陣請修改n的值
int main()
}void
out(double matrix[
][n]
,const char* s)
;void
product
(double matrix1[
][n]
, double matrix2[
][n]
, double matrix_result[n]
[n])
;//初始化
double l[n]
[n], u[n]
[n], c[n]
[n], ad_l[n]
[n], ad_u[n]
[n];
memset
(inv_a,
0, n * n *
sizeof
(double));
memset
(l,0
, n * n *
sizeof
(double));
memset
(u,0
, n * n *
sizeof
(double));
memset
(c,0
, n * n *
sizeof
(double));
memset
(ad_l,
0, n * n *
sizeof
(double));
memset
(ad_u,
0, n * n *
sizeof
(double));
//lu分解
for(int i =
0; i < n; i++
) l[i]
[i]=1;
for(int i =
0; i < n -
1; i++
)for
(int j = i +
1; j < n; j++)}
u[n -1]
[n -1]
= a[n -1]
[n -1]
;for
(int i =
0; i < n -
1; i++
) u[n -1]
[n -1]
-= l[n -1]
[i]* u[i]
[n -1]
;//矩陣l
//out(l, "矩陣l");
//矩陣u
//out(u, "矩陣u");
//l*u
//product(l,u, c);
//out(c, "矩陣l*u");
//u的逆矩陣
for(int i =
0; i < n; i++
)else
if(j - i ==1)
else
if(j - i >=2)
while
(next_permutation
(permutation, permutation + j - i));
ad_u[i]
[j]= sum * deltas_aii;}}
} double det_u =1;
for(int k =
0; k < n; k++
) det_u *= u[k]
[k];
if(det_u <
1e-16
)for
(int i =
0; i < n; i++
)for
(int j =
0; j < n; j++
) ad_u[i]
[j]/= det_u;
//inv_u
//out(ad_u, "inv_u");
//u*u-1
//memset(c, 0, n*n * sizeof(double));
//product(u, ad_u, c);
//out(c, "矩陣u*u-1");
//l的逆矩陣
for(int i =
0; i < n; i++
)else
if(i - j ==1)
else
if(i - j >=2)
while
(next_permutation
(permutation, permutation + i - j));
ad_l[i]
[j]= sum * deltas_aii;}}
} double det_l =1;
for(int k =
0; k < n; k++
) det_l *= l[k]
[k];
if(det_u <
1e-16
)for
(int i =
0; i < n; i++
)for
(int j =
0; j < n; j++
) ad_l[i]
[j]/= det_l;
//矩陣inv_l
//out(ad_l, "矩陣inv_l=");
//l*l-1
//memset(c, 0, n*n * sizeof(double));
//product(l, ad_l, c);
//out(c, "矩陣l*inv_l=");
//矩陣a
out(a,
"矩陣a=");
//inv_a
product
(ad_u, ad_l, inv_a)
;out
(inv_a,
"逆矩陣inv_a=");
//驗證a*inv_a
memset
(c,0
, n * n *
sizeof
(double));
product
(a, inv_a, c)
;out
(c,"矩陣a*inv_a=");
system
("pause");
return0;
}void
out(double matrix[
][n]
,const char* s)
cout << endl;
}void
product
(double matrix1[n]
[n], double matrix2[n]
[n], double matrix_result[n]
[n])
矩陣LU分解的MATLAB與C 實現
矩陣的lu分解目的是將乙個非奇異矩陣 a 分解成 a lu 的形式,其中 l 是乙個主對角線為 1 的下三角矩陣 u 是乙個上三角矩陣。比如 a begin 1 2 4 3 7 2 2 3 3 end 我們最終要分解成如下形式 a l cdot u begin 1 0 0 3 1 0 2 1 1 e...
矩陣sum 矩陣LU分解的MATLAB與C 實現
矩陣的lu分解目的是將乙個非奇異矩陣 比如 現在主要的問題是如何由矩陣 計算得到矩陣 和 呢?我們將在下面詳細討論。首先從矩陣 入手,因為它是乙個上三角矩陣,所以很容易想到高斯消元法,依次把矩陣 主對角線左下角的元素消為 就得到 了。然後計算矩陣 這裡有個技巧,可以這樣想,正是因為有了 所以 的左下...
PLU 分解以及求逆矩陣
plu 分解是對lu分解的一種改進,其增加了選主元的操作增加了計算的穩定性,及在第i次迴圈中將 j w here max a i n,i j where max a i n,i j wher e max a i n i 行和第i行進行交換來比避免對角元素出現0的情況,計算結果 p a lu pa l...