c 矩陣求逆的lu分解實現

2021-10-01 19:17:26 字數 3809 閱讀 9984

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...