Numerical Methods for Material Scientists 20080930 3rd presentation by Lee, jisan

Assignment n x n linear system Gauss-Jordan Elimination: Using Augmented Matrix

Pivoting Pivoting: Pivot 이 0이면 Gauss elimination 이 불가능하고, Computer가 처리하는 수의 체계는 정밀도에 한계가 있기 때문에 0에 가까운 경우에는 반올림오차가 발생할 수 있다. 따라서 각 행간의 순서를 바꾸는 것을 Pivoting이라고 한다. Partial Pivoting: Pivot 요소 아래 열에서 가장 큰 계수를 찾은 후, 그 항이 Pivot 이 되도록 행의 위치를 교환하는 방법. Scaled Partial Pivoting: 단순히 각 행의 첫 계수만을 비교하는 것이 아니라 행을 Normarlizing 한 후 Partial pivoting 을 하는 방법.

C++로 n x n matrix 받기 # include<stdlib.h> Void main() { int r, i; double **a; //2d 배열 double *x; //1d 배열 scanf(“%d”, &r); a = (double **)malloc(r * sizeof(double)); x=(double*)malloc(r*sizeof(double)); for(i = 0 ; i < r ; i++) a[i] = (double *)malloc((r+1)*sizeof(double)); } -> … Value 입력받으면 됩니다.

1. Scaled Partial Pivoting 각 행의 계수 중 가장 큰 값으로 normalization 한 값을 비교하여 pivoting 한다. //scale factor 구하기: for (i=0; i<r; i++) { s[i] = abs(a[i][0]); for (j=1; j<r; j++) { if (abs(a[i][j])>s[i]) s[i] = abs(a[i][j]); }

1. Scaled Partial Pivoting 각 행의 계수 중 가장 큰 값으로 normalization 한 값을 비교하여 pivoting 한다. void scaled_pivot(double **a, double *s, int r, int n) { double dum=0; int i, maxi=n; for (i=n; i<r; i++) if (abs(a[i][n]/s[i])>dum) { dum = abs(a[i][n]/s[i]); maxi = i; } if (maxi!=n) for (i=n; i<(r+1); i++) { dum = a[maxi][i]; a[maxi][i] = a[n][i]; a[n][i] = dum; dum = s[maxi]; s[maxi] = s[n]; s[n] = dum; 각 행을 Normalizing 한 후 가장 큰 수를 저장. 만약 현재 Pivot 이 가장 큰 수가 아니면 아래 행과 교체 한다.

2. Gauss Elimination 위 행에 적당한 factor 를 곱하여 제일 첫 열을 소거시키면서 upper triangular matrix형태로 만든다. void elimination(double **a, int r, int n, int *err) { double factor1; int i, j; *err = -1; for (i=n+1; i<r; i++) { factor1 = a[i][n]/a[n][n]; for (j=n; j<=r; j++) { a[i][j] = a[i][j] - (factor1 * a[n][j]); } for (i=n+1; i<r; i++) if (a[i][n+1]!=0) *err = 0; Factor 계산하고 아래 행에 pivot항을 곱한 값을 뺀다. Pivot 이 0이 되면 error 출력 (err= -1 로 두고, 여기서 0으로 다시 바뀌지 않으면 error)

Result : Pivoting + Gauss Elimination Ex) a 4x4 System

3. Back Substitution Gauss Elimination 한 결과를 바탕으로 최종 해를 구한다. void substitution (double **a, int r, double *x) { int i, j; double sum=0; for (i=(r-1); i>=0; i--) { sum = 0; for (j=i+1; j<r; j++) { sum = sum + a[i][j]*x[j]; } x[i] = (a[i][r] - sum) / a[i][i];

4. Inverse Matrix A와 Identity의 Augmented Matrix에서 Back subsitution 까지는 같은 방법으로 하고, 그 후 A가 Identity가 될 때까지 역으로 elimination -> Inverse Mtx. 즉, A를 I로 바꾼 Operation B가 A의 Inverse가 된다. void inversion (double **a, double **b, int r) { int i, j, k; double dum1, dum2; for (i=0; i<r; i++) { dum1 = a[i][i]; for (j=0; j<r; j++) a[i][j] = a[i][j] / dum1; b[i][j] = b[i][j] / dum1; } for (i=r-1; i>0; i--) for (j=i-1; j>=0; j--) { dum2 = a[j][i]; for (k=r-1; k>=0; k--) { b[j][k] = b[j][k] - dum2*b[i][k]; a[j][k] = a[j][k] - dum2*a[i][k]; Matrix A의 Pivot 을 1로 맞춰준다. 다시 아래로부터 elimination

Result : Back Substitution & Inverse Matrix

5. Analysis: Pivoting to reduce error. Computer는 한정된 유효수자만 취급하기 때문에 항상 반올림오차가 발생할 수 있다. 대략 100개 이상의 방정식을 다루거나, 또는 계수 크기 자체가 반올림 오차에 영향을 끼친다. 다음과 같은 두 개의 간단한 Linear System 을 생각해 보자. 0.000000000000003X1 + 3X2 = 2.000000000000001 X1 + X2 = 1 Real solution: X1 =1/3, X2 =2/3 2X1 + 100000000000000000X2 = 100000000000000000 X1 + X2 = 2 Real solution: X1 =1.00000000000000002, X2 =0.99999999999999998

5. Analysis: Pivoting to reduce error. 0.000000000000003X1 + 3X2 = 2.000000000000001 X1 + X2 = 1 Real solution: X1 =1/3, X2 =2/3

5. Analysis: Pivoting to reduce error. 0.000000000000003X1 + 3X2 = 2.000000000000001 ……1) X1 + X2 = 1 ……2) 왜 이렇게 큰 오차가??? 식에 1/0,000000000000003곱하고 두번째 식으로부터 X1을 소거하면, -999999999999999 X2 = -666666666666666 X2 =2/3 이 결과를 1)식에 대입하여 X1을 구하면, X1 = ( 2.000000000000001 - 3(2/3) ) / 0.000000000000003 여기서 이 값은 뺄셈의 무효화로 인해 유효숫자의 개수에 매우 민감하다. 만약 Pivoting 을 수행한다면, X2 = 2/3 은 그대로고 이를 1)에 대입시키면 X1 = (1 – 2/3) / 1 을 얻을 수 있다. 이 때 위 결과보다 훨씬 정확한 해를 얻을 수 있다.

5. Analysis: Pivoting to reduce error. 2X1 + 100000000000000000X2 = 100000000000000000 X1 + X2 = 2 Real solution: X1 =1.00000000000000002, X2 =0.99999999999999998

5. Analysis: Pivoting to reduce error. 2X1 + 100000000000000000X2 = 100000000000000000 ……1) X1 + X2 = 2 ……2) 왜 이렇게 큰 오차가??? Partial Pivoting 만 한 경우, -49999999999999999 X2 = -49999999999999998 X2 =1 이 결과를 1)식에 대입하여 X1을 구하면, X1 = ( 100000000000000000 - 100000000000000000) / 2 X1 = 0 이런 오차가 나오는 이유는 order가 매우 큰 수 간의 연산에서 1은 무시되기 때문이다. 2) Scaled Partial Pivoting한 경우, X2 =1 를 식 2)에 대입하면, X1 = 2-X2 = 1 역시 위 결과보다 훨씬 정확한 해를 얻을 수 있다.

6. Summary C++을 이용하여 n x n Linear System 을 Gauss-Jordan Elimination 방법을 이용하여 해를 구하고 역행렬을 구하였다. 이 방법에서 최대한 오차를 줄이기 위해 Scaled partial Pivoting 을 수행하였다. Double Precision 에서도 severe한 오차가 날 수 있는 두 가지 간단한 Linear System 을 통해 Pivoting 을 하지 않았을 때의 해와 Partial Pivoting을 하였을 때, 그리고 Scaled Partial Pivoting 을 하였을 때의 해를 비교하였다. Pivoting 이 적절하지 않을 때 오차가 커지는 이유는, 비슷한 수 끼리의 뺄셈이나, order차이가 매우 큰 수 간의 덧셈, 뺄셈을 컴퓨터가 구분하지 못하기 때문이다. 결국 Computer 의 유효 숫자 처리에도 한계가 있기 때문에 Scaled Partial Pivoting 을 실행하여 반올림 오차를 줄여야 한다.