============================================================= Source code ============================================================= #include "stdafx.h" #include "math.h" int _tmain(int argc, _TCHAR* argv[]) { int i, j, k, m; double a[8][8] ={35.836, 20.0433, 26.3189, 25.1571, 15.4733, 20.6763, 24.0837, 23.301, 20.0433, 14.6742, 15.0619, 16.9476, 6.8295, 11.9028, 13.7083, 14.3384, 26.3189, 15.0619, 32.5065, 26.9765, 21.3285, 19.3257, 26.2256, 24.4308, 25.1571, 16.9476, 26.9765, 30.1095, 17.3402, 14.3453, 22.386, 19.7201, 15.4733, 6.8295, 21.3285, 17.3402, 19.1458, 13.5818, 16.8893, 16.952, 20.6763, 11.9028, 19.3257, 14.3453, 13.5818, 21.7242, 14.2585, 23.6775, 24.0837, 13.7083, 26.2256, 22.386, 16.8893, 14.2585, 25.8063, 19.4059, 23.301, 14.3384, 24.4308, 19.7201, 16.952, 23.6775, 19.4059, 27.9643}; double b[8] = {9.5974, 3.4039, 5.8527, 2.2381, 7.5127, 2.5510, 5.0596, 6.9908}; double l[8][8]; double lt[8][8]; double A[8][8]; double y[8]; double x[8]; double deta, detu; double AC[8][8]; double ax[8]; double xx[8]; double u1[8][8]; double l1[8][8]; for (i = 0; i <8; i++){ for (j = 0; j <8; j++){ l[i][j]=0; lt[i][j]=0; } } for (i = 0; i <8; i++){ y[i]=b[i]; } for(i=0; i<8; i++){ for(j=0; j<8; j++){ AC[i][j]=0; AC[i][j]=a[i][j]; } } /* 1(1.1) */ for (k = 0; k < 8-1; k++){ l[k][k]=sqrt(a[k][k]); for(m=k+1; m<8; m++){ l[m][k]=a[m][k]/l[k][k]; for(j=k+1; j<8; j++){ for(i=m; i<8; i++){ a[i][j]-=l[i][k]*l[j][k]; } } } } l[7][7]=sqrt(a[7][7]); printf("1(1.1) \n"); printf("Matrix L: \n"); for(i=0; i<8; i++){ for(j=0; j<8; j++){ printf(" %7.4f", l[i][j]); } printf("\n"); } for(i=0; i<8; i++){ for(j=0; j<8; j++){ lt[i][j]=l[j][i]; } } printf("Matrix Lt: \n"); for(i=0; i<8; i++){ for(j=0; j<8; j++){ printf(" %7.4f", lt[i][j]); } printf("\n"); } for(i=0; i<8; i++){ for(j=0; j<8; j++){ A[i][j]=0; for (k = 0; k < 8; k++) { A[i][j]+=l[i][k]*lt[k][j]; } } } printf("Matrix L*Lt=: \n"); for(i=0; i<8; i++){ for(j=0; j<8; j++){ printf(" %7.4f", A[i][j]); } printf("\n"); } /* 1(1.2) */ for (i = 0; i <8; i++) { y[i] /= l[i][i]; for (j = 0; j =0; i--) { x[i] /= l[i][i]; for (j = i+1; j < 8; j++){ x[i] -= l[j][i] * x[j]/l[i][i]; } } printf("1(1.2) \n"); printf("x= \n"); for(i=0; i<8; i++){ printf(" %7.4f", x[i]); printf("\n"); } for(i=0; i<8; i++){ ax[i]=0; for (k = 0; k < 8; k++) { ax[i]+=AC[i][k]*x[k]; } } for(i=0; i<8; i++){ xx[i]=ax[i]-b[i]; } printf("Ax-b= \n"); for(i=0; i<8; i++){ printf(" %7.4f", xx[i]); printf("\n"); } /* 1(1.3) */ deta=1; for(i=0; i<8; i++){ deta *= l[i][i]; } deta=deta*deta; printf("1(1.3) \n"); printf("detA = %13.8f\n", deta); /* LU */ for( i = 0; i < 8; i++ ){ for( j = 0; j < 8; j++ ){ if( i == j ){ l1[i][j] = 1; } else{ l1[i][j] = 0; } u1[i][j] = AC[i][j]; } } for( i = 0; i < 8; i++ ){ for( j = i+1; j < 8; j++ ){ l1[j][i] = u1[j][i] / u1[i][i]; for( k = i; k < 8; k++ ){ u1[j][k] = u1[j][k] - l1[j][i]*u1[i][k]; } } } detu=1; for(i=0; i<8; i++){ detu *= u1[i][i]; } printf("detU = %13.8f\n", detu); return 0; } ============================================================= Results ============================================================= 1(1.1) ============================================================= Matrix L: 5.9863 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 3.3482 1.8611 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 4.3965 0.1835 3.6254 0.0000 0.0000 0.0000 0.0000 0.0000 4.2024 1.5459 2.2664 2.2187 0.0000 0.0000 0.0000 0.0000 2.5848 -0.9805 2.7982 0.7444 1.7662 0.0000 0.0000 0.0000 3.4539 0.1818 1.1329 -1.3604 1.5147 2.0817 0.0000 0.0000 4.0231 0.1280 2.3485 -0.0187 0.0330 -1.1512 1.6620 0.0000 3.8924 0.7017 1.9830 -0.9990 1.5706 1.9799 0.7270 0.6897 Matrix Lt: 5.9863 3.3482 4.3965 4.2024 2.5848 3.4539 4.0231 3.8924 0.0000 1.8611 0.1835 1.5459 -0.9805 0.1818 0.1280 0.7017 0.0000 0.0000 3.6254 2.2664 2.7982 1.1329 2.3485 1.9830 0.0000 0.0000 0.0000 2.2187 0.7444 -1.3604 -0.0187 -0.9990 0.0000 0.0000 0.0000 0.0000 1.7662 1.5147 0.0330 1.5706 0.0000 0.0000 0.0000 0.0000 0.0000 2.0817 -1.1512 1.9799 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 1.6620 0.7270 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.6897 Matrix L*Lt=: 35.8360 20.0433 26.3189 25.1571 15.4733 20.6763 24.0837 23.3010 20.0433 14.6742 15.0619 16.9476 6.8295 11.9028 13.7083 14.3384 26.3189 15.0619 32.5065 26.9765 21.3285 19.3257 26.2256 24.4308 25.1571 16.9476 26.9765 30.1095 17.3402 14.3453 22.3860 19.7201 15.4733 6.8295 21.3285 17.3402 19.1458 13.5818 16.8893 16.9520 20.6763 11.9028 19.3257 14.3453 13.5818 21.7242 14.2585 23.6775 24.0837 13.7083 26.2256 22.3860 16.8893 14.2585 25.8063 19.4059 23.3010 14.3384 24.4308 19.7201 16.9520 23.6775 19.4059 27.9643 ============================================================= 1(1.2) ============================================================= x= 4.4683 5.5626 3.9636 -7.7353 5.4125 -18.7650 -7.6187 13.5611 Ax-b= -0.0000 -0.0000 0.0000 -0.0000 0.0000 -0.0000 -0.0000 -0.0000 ============================================================= 1(1.3) ============================================================= L = 1.0000 0 0 0 0 0 0 0 0.5593 1.0000 0 0 0 0 0 0 0.7344 0.0986 1.0000 0 0 0 0 0 0.7020 0.8306 0.6252 1.0000 0 0 0 0 0.4318 -0.5268 0.7718 0.3355 1.0000 0 0 0 0.5770 0.0977 0.3125 -0.6131 0.8576 1.0000 0 0 0.6721 0.0688 0.6478 -0.0085 0.0187 -0.5530 1.0000 0 0.6502 0.3770 0.5470 -0.4503 0.8892 0.9511 0.4374 1.0000 U = 35.8360 20.0433 26.3189 25.1571 15.4733 20.6763 24.0837 23.3010 0 3.4639 0.3416 2.8771 -1.8248 0.3385 0.2381 1.3060 0 0 13.1434 8.2167 10.1445 4.1071 8.5144 7.1891 0 0 0 4.9226 1.6517 -3.0183 -0.0416 -2.2165 0 0 0 0 3.1196 2.6752 0.0582 2.7740 0 0 0 0 0 4.3334 -2.3964 4.1215 0 0 0 0 0 0 2.7620 1.2082 0 0 0 0 0 0 0 0.4758 detA = 142653.24173371 detU = 142653.24173371 ============================================================= 2. A = QR, where Q is orthogonal matrix such that Q^tQ = I, R is a upper triangular matrix => Q^tA = Q^tQR = R. Thus, Ax = b => Q^tAx = Q^tb => Rx = Q^tb which is easy solve using backward substitution because R is upper triangular. ============================================================= 3. Q^tQ = I => det(Q^tQ) = det(I) = 1 => det(Q^t)det(Q) = 1 => (det(Q))^2 = 1 => det(Q) = 1 or -1. ============================================================= 4. A = QR => det(A) = det(QR) = det(Q)det(R) = det(R) or -det(R) = r11*r22*....*rnn or -r11*r22*....*rnn.