LU 分解 解线性性方程组

Code:

1. #include<iostream>   
2. using namespace std;  
3. #define M 50   
4. int  main(){   
5.   
6.    double A[M][M],L[M][M],U[M][M];   
7.    cout<<"请输入矩阵的行列数："<<endl;   
8.    int n;                                   //n矩阵的行列数   
9.    cin>>n;   
10.    int num=n*n;                             //num矩阵的总数据个数   
11.    cout<<'/n'<<"请输入矩阵："<<endl;   
12.   
13.    double temp;   
14.    for(int i=1;i<n+1;i++)   
15.     for(int j=1;j<n+1;j++)   
16.     {   
17.       cin>>A[i][j];   
18.     }   
19.   cout<<"你输入的矩阵如下："<<endl;   
20.    for(int i=1;i<n+1;i++){   
21.     for(int j=1;j<n+1;j++)   
22.     {   
23.      cout<<A[i][j]<<' ';   
24.   
25.     }   
26.     cout<<endl;   
27.    }   
28.   
29.   for(int i=1;i<n+1;i++)   
30.      L[i][i]=1; //对角线 已经 知道全部为 1；   
31.   
32.   for(int i=1;i<n+1;i++)              //求U矩阵的第一行元素   
33.        U[1][i]=A[1][i]; //U 的第一行 等于 a 的第一行   
34.   for(int i=2;i<n+1;i++)   
35.   
36.    L[i][1]=A[i][1]/U[1][1];        //求L矩阵的第一列元素   
37.   for(int r=2;r<n+1;r++)   
38.   {   
39.        for(int i=r;i<n+1;i++)   
40.        {   
41.         temp=0;                     //用来做临时变量 来存取L[r][k]*U[k][i]的和 其中k=1...........r-1   
42.         for(int k=1;k<r;k++)   
43.          temp+=L[r][k]*U[k][i];   
44.         U[r][i]=A[r][i]-temp;   
45.   
46.        }   
47.   
48.        for(int i=r+1;i<n+1;i++)   
49.        {   
50.          temp=0;   
51.          for(int k=1;k<r;k++)   
52.          temp+=L[i][k]*U[k][r];   
53.          L[i][r]=(A[i][r]-temp)/U[r][r];   
54.        }   
55.   }   
56.   cout<<"LU分解的结果如下："<<endl;   
57.   
58.   cout<<"L矩阵："<<endl;   
59.   for(int i=1;i<n+1;i++){   
60.    for(int j=1;j<n+1;j++){   
61.     if(i>j)   
62.      cout<<L[i][j]<<' ';  //L全部在 i<j 部分   
63.     else  
64.      if(i==j)   
65.       cout<<'1'<<' ';    //对角线 为1   
66.      else  
67.       cout<<'0'<<' ';   // 其他部分 为 0   
68.   
69.    }   
70.    cout<<endl;   
71.   }   
72.   cout<<"U矩阵："<<endl;   
73.   for(int i=1;i<n+1;i++)   
74.   {   
75.    for(int j=1;j<n+1;j++)   
76.     if(i>j)   
77.      cout<<'0'<<' ';   
78.     else  
79.      cout<<U[i][j]<<' ';  //U矩阵 在 右上角 这时的 i<j;   
80.   
81.    cout<<endl;   
82.   }   
83. }   
84.