數值分析作業 Lu-decomposition

數值分析 第五次作業 LU D eco mpo sition 組別： 第六組 學生： 陳器奇、謝鎮宇 學號：9301062a、9301035a 程式碼(使用 C 語言)： #include #include #include #define double

<stdio.h> <stdlib.h> <math.h> dim 10

//

說明(1)

//

說明(2)

//

說明(3)

//

說明(4)

printf( "\nInput the dimension of Matrix: " ); scanf( "%d",&n );

//

說明(5)

if(n>=dim) { printf( "Dimension over preset value %d\n",dim ); exit(0); } printf( "\n Input the a matrix\n\n" );

//

說明(6)

//

說明(7)

//

說明(8)

a[dim][dim], L[dim][dim], U[dim][dim], x[dim], b[dim], bp[dim];

mprint( int , int , double [][dim] , double [] ); main() { int n,i,j,k; double temp;

for(i=1; i<=n; i++) { for(j=1;j<=n;j++) { printf( "Input the element of a[%d,%d]: ",i,j ); scanf( "%lf",&a[i][j] ); } } printf( "\n Input the b vector\n\n" ); for(i=1; i<=n ; i++) { printf( "input the value of b[%d]: ",i ); scanf( "%lf",&b[i] );

} mprint(n,n,a,b);

//

說明(9)

//

說明(10)

for(i=1; i<=n; i++) { L[i][1]=a[i][1]; } for(j=1; j<=n; j++) { U[1][j]=a[1][j]/L[1][1]; } for(j=2; j<=n; j++) { for(i=j; i<=n;i++) { temp=0.0; for(k=1;k<=j-1;k++) { temp=temp+L[i][k]*U[k][j]; } L[i][j]=a[i][j]-temp; } U[j][j]=1.0;

//

說明(11)

//

說明(12)

for(i=j+1;i<=n;i++) { temp=0.0; for(k=1; k<=j-1;k++) { temp=temp+L[j][k]*U[k][i]; } U[j][i]=(a[j][i]-temp)/L[j][j]; }

//

說明(13)

//

說明(14)

/*

開始 LU 分解程序

*/

for(i=0; i
} /*

Backward substitution (反代法)

bp[1]=b[1]/L[1][1]; for(i=2; i<=n; i++) { temp=0.0; for(k=1; k<=i-1;k++) { temp=temp+L[i][k]*bp[k]; } bp[i]=(b[i]-temp)/L[i][i];

*/

} x[n]=bp[n]; for(j=n-1; j>=1; j--) { temp=0.0; for(k=j+1; k<=n; k++) { temp=temp+U[j][k]*x[k]; } x[j]=bp[j]-temp; }

}

mprint(n,n,L,b); mprint(n,n,U,bp); for(j=1; j<=n; j++) { printf(" x[%d]=%10.5lf\n",j,x[j]); } /* 結束主程式 */

/* 副程式 */ mprint(int m,int n,double aa[][dim],double bb[]) { int i,j; for(i=1; i<=m; i++){ for(j=1; j<=n;j++){ printf(" %10.5e ",aa[i][j]); } printf(" %10.5e \n",bb[i]); } }

//

說明(15)

//

說明(16)

//

說明(17)

//

說明(18)

程式說明： LU Decomposition 定義：將一個矩陣 A 分解成上三角矩陣 U 與下三角矩陣 L。 解法方程式為： Ax = b =>LUx = b 設定 Ux=b'，求 x，由下列兩個方程式求解: Lb' = b Ux = b' 先算出 b'，再代回算出 x(反代法)。 矩陣的 LU 分解(以四乘四為例)

將 L 矩陣的一，二，三，四列與 U 矩陣的第一行相乘得到: L11 = a11 , L21 = a21 , L31 = a31 , L41 = a41 將 L 矩陣的第一列和 U 矩陣的二，三，四行相乘則得到：

同理 L 矩陣的二，三，四列乘以 U 的矩陣第二行可得：

同理可以推出：

所以我們可以得出通 式：

當 j=1 時, Li1 = ai1 。而當 i=1 時,

最後使用 Backward Substitution(反代法)求出 x：

LU 分解的優點：不需儲存 LU 中的 0 與 1，原本的矩陣 A 即能表示 LU 兩個矩陣

說明(1)：定義 dim 為 10，為我們設定的矩陣維度最大值。 說明(2)：宣告 a, L, U 為二維矩陣，x, b, bp 為一維矩陣。 說明(3)：使用 mprint 副程式。 說明(4)：宣告 n,i,j,k 為整數，temp 為雙精度浮點數。 說明(5)：螢幕顯示輸入矩陣的維度。 說明(6)：如果掃描輸入的 n 大於或等於 dim，螢幕顯示輸入維度超過預定值。 說明(7)：螢幕顯示逐次輸入 A 矩陣的元素，掃描。 說明(8)：螢幕顯示逐次輸入 b 向量矩陣的元素值，掃描。 說明(9)：使用 mprint 副程式，螢幕顯示 a,b 矩陣。 說明(10)：設 L,U 矩陣內元素為 0。 說明(11)：計算出下三角矩陣 L 的方法。 說明(12)：設上三角矩陣 U 對角線元素值為 1。 說明(13)：計算出上三角矩陣 U 的方法。 說明(14)：計算出 b' 矩陣的方法。 說明(15)：j 遞減反代去求出 x。 說明(16)：螢幕顯示 L,U,b,b' 的值。 說明(17)：螢幕顯示最後算出來的 x 矩陣值。

說明(18)：副程式 mprint，逐次顯示出傳進來的 aa 及 bb 值。 程式執行結果及討論： 題目希望我們求以下聯立方程式之 x , y ,

z 的解：

x + y + z = -2 2x - y + 3z = 5 3x + 2y - 2z = 1

[ 情況一 ：使用 double 型 別，在說明(8)及說明(9)中，a[i][j]及 b[i]掃描為雙精度浮點數 ] Input the dimension of Matrix: 3 Input the a matrix Input Input Input Input Input Input Input Input Input

the the the the the the the the the

element element element element element element element element element

of of of of of of of of of

a[1,1]: 1 a[1,2]: 1 a[1,3]: 1 a[2,1]: 2 a[2,2]: -1 a[2,3]: 3 a[3,1]: 3 a[3,2]: 2 a[3,3]: -2

Input the b vector input the value of b[1]: -2 input the value of b[2]: 5 input the value of b[3]: 1 1.00000e+00 1.00000e+00 1.00000e+00 -2.00000e+00 2.00000e+00 -1.00000e+00 3.00000e+00 5.00000e+00 3.00000e+00 2.00000e+00 -2.00000e+00 1.00000e+00 1.00000e+00 2.00000e+00 3.00000e+00

0.00000e+00 0.00000e+00 -2.00000e+00 -3.00000e+00 0.00000e+00 5.00000e+00 -1.00000e+00 -5.33333e+00 1.00000e+00

1.00000e+00 1.00000e+00 1.00000e+00 -2.00000e+00 0.00000e+00 1.00000e+00 -3.33333e-01 -3.00000e+00

0.00000e+00 x[1]= x[2]= x[3]=

0.00000e+00

1.00000e+00

-7.50000e-01

2.00000 -3.25000 -0.75000

我們可以算出下三角矩陣 L： [ 1.00000e+00 [ 2.00000e+00 [ 3.00000e+00

0.00000e+00 0.00000e+00 ] -3.00000e+00 0.00000e+00 ] -1.00000e+00 -5.33333e+00 ]

上三角矩陣 U： [ 1.00000e+00 1.00000e+00 1.00000e+00 ] [ 0.00000e+00 1.00000e+00 -3.33333e-01 ] [ 0.00000e+00 0.00000e+00 1.00000e+00 ] 可以算出 x = 2，y = -3.25，z = -0.75 [ 情況二 ：使用 double 型別，在說明(8)及說明(9)中， a[i][j]及 b[i]掃描為單精度浮點數 ] Input the dimension of Matrix: 3 Input the a matrix Input Input Input Input Input Input Input Input Input

the the the the the the the the the

element element element element element element element element element

of of of of of of of of of

a[1,1]: 1 a[1,2]: 1 a[1,3]: 1 a[2,1]: 2 a[2,2]: -1 a[2,3]: 3 a[3,1]: 3 a[3,2]: 2 a[3,3]: -2

Input the b vector input the value of b[1]: -2 input the value of b[2]: 5 input the value of b[3]: 1 5.26354e-315 5.26354e-315 5.26354e-315 5.30499e-315 1.58735e-314 5.32571e-315 5.32571e-315 5.30499e-315 1.59150e-314 5.26354e-315 0.00000e+00 5.30499e-315 1.05685e-314 5.32571e-315 -2.07226e-317

1.59150e-314 5.35680e-315 5.26354e-315

0.00000e+00 1.59150e-314 0.00000e+00 5.35680e-315 1.05893e-314 5.26354e-315

1.00000e+00 1.00000e+00 1.00000e+00 3.02362e+00 0.00000e+00 1.00000e+00 1.96078e-03 -1.01088e+00 0.00000e+00 0.00000e+00 1.00000e+00 -1.02560e+00 x[1]= x[2]= x[3]=

5.05808 -1.00887 -1.02560

我們可以算出下三角矩陣 L： [ 5.26354e-315 0.00000e+00 [ 5.30499e-315 1.05685e-314 [ 5.32571e-315 -2.07226e-317

0.00000e+00 ] 0.00000e+00 ] 1.05893e-314 ]

上三角矩陣 U： [ 1.00000e+00 1.00000e+00 1.00000e+00 ] [ 0.00000e+00 1.00000e+00 1.96078e-03 ] [ 0.00000e+00 0.00000e+00 1.00000e+00 ] 最後可以算出 x = 5.05808，y = -1.00887，z = -1.02560 [ 情況三 ：將 double 改成 float，在說明(8)說明(9)中，a[i][j]及 b[i]掃描為單精度浮點數 ] Input the dimension of Matrix: 3 Input the a matrix Input Input Input Input Input Input Input Input Input

the the the the the the the the the

element element element element element element element element element

of of of of of of of of of

a[1,1]: 1 a[1,2]: 1 a[1,3]: 1 a[2,1]: 2 a[2,2]: -1 a[2,3]: 3 a[3,1]: 3 a[3,2]: 2 a[3,3]: -2

Input the b vector input the value of b[1]: -2 input the value of b[2]: 5 input the value of b[3]: 1 1.00000e+00 2.00000e+00 3.00000e+00

1.00000e+00 1.00000e+00 -2.00000e+00 -1.00000e+00 3.00000e+00 5.00000e+00 2.00000e+00 -2.00000e+00 1.00000e+00

1.00000e+00 2.00000e+00

0.00000e+00 0.00000e+00 -2.00000e+00 -3.00000e+00 0.00000e+00 5.00000e+00

3.00000e+00

-1.00000e+00

-5.33333e+00

1.00000e+00

1.00000e+00 1.00000e+00 1.00000e+00 -2.00000e+00 0.00000e+00 1.00000e+00 -3.33333e-01 -3.00000e+00 0.00000e+00 0.00000e+00 1.00000e+00 -7.50000e-01 x[1]= x[2]= x[3]=

2.00000 -3.25000 -0.75000

結果與情況一相同 結論 ：從情況一及三中，我們可以解出 x = 2，y = -3.25，z = -0.75 (與課本後面解答相 同)。但如果我們使用 double 當作變數宣告型別，a[i][j]及 b[i]掃描為單精度浮點數，則會 使最後的 x，y，z 值產生不小的誤差：x = 5.05808，y = -1.00887，z = -1.02560，與正確 答案相差不少，所以我們建議使用變數宣告型別與掃描型別最好是同樣的精度浮點數(單 精度浮點數或是雙精度浮點數)，以提高計算正確度。

～完～ ～感謝閱讀～