재료수치해석 강릉대학교 금속공학과 정 효 태 Numerical Analysis

차 례 F 비선형 방정식 F 연립 1 차 방정식 F 행렬의 고유치문제 F 수치 보간법 F 회 귀 분 석F 회 귀 분 석 F 수 치 적 분F 수 치 적 분 F 상미분 방정식 F 편미분 방정식 F Monte Carlo 법 F 최적화 수법 F 선 형 계 획F 선 형 계 획

1 장 비선형방정식 (Non-Linear Equations) F f(x) = 0 * e -x - cos(  x/2) = 0 * 3x 4 - x 3 - 5x 2 - x + 4 = 0 4 Incremental Search 법 4 Bisection 법 ( 반분법 ) 4 False Position 법 (Regula Falsi 법 ) 4 Newton-Raphson 법 4 Newton’s Second-Order 법 4 Bairstow 법 4 Bisection 법에 의한 연립비선형방정식의 해석 4 Newton-Raphson 법에 의한 연립비선형방정식의 해석

1.1 Incremental Search 법 F 변수값을 적정구간내에서 미세하게 변화시켜 해를 구함 F 예 ] f(x) = a 3 x 3 + a 2 x 2 + a 1 x 1 + a 0 = 0 0 y x y=f(x) xx xixi xi+xxi+x f(x i ) f(x i +  x) F f(x i )*f(x i +  x) < 0 이면 * x i 와 x i +  x 사이에 근이 존재 *  x 를 감소시키고, * x i 를 초기치로 설정하여 * 증분  x 가 허용오차보다 작을 때까지 반복계산

F 계산순서 1. 초기근사해 x i 를 정한다. 2. 초기증분  x 를 정한다. 3. f(x i ) 및 f(x i +  x) 를 계산한다. 4. f(x i )*f(x i +  x) 의 부호를 계산한다. 5.  f(x i )*f(x i +  x) > 0 이면 x i  x i +  x 입력후 순서 3 으로  f(x i )*f(x i +  x) < 0 이면  x <  이면 x i +  x 를 근으로 하고 계산 끝  x >  이면  x   x*0.1 입력후 순서 3 으로  f(x i )*f(x i +  x) = 0 이면 x i +  x 를 근으로 하고 계산 끝

/* Program 1-1 */ /* THE INCREMENTAL-SEARCH METHOD */ /* */ #include #define f(x) a[3]*pow((x),3)+a[2]*pow((x),2)+a[1]*(x)+a[0] main() { double a[4]; double xi, xd, y1, y2, eps=1.0e-4, del_x=0.1; int it_max, it=0; printf("%c Input a[3 - 0]:"); scanf("%lf%lf%lf%lf", &a[3], &a[2], &a[1], &a[0]); printf("%cInput [Max Iteration] & [initial x]:"); scanf("%d%lf", &it_max, &xi); printf("\n Initial dx = %f Tolerance = %f\n", del_x, eps);

while(1){ it++; xd=xi+del_x; y1=f(xi); y2=f(xd); if( y1*y2>0 )xi=xd; if( y1*y2<0 ) if( del_x<eps )break; else del_x=del_x/10.0; if( y1*y2==0 )break; if( it>it_max ){ printf("iteration= %d Check Error!!!!",it); exit(1); } printf("\nf(x) = %10.3lfx^3 + %10.3lfx^2 + %10.3lfx + %10.3lf\n", a[3], a[2], a[1], a[0]); printf("First= %10.3lf \n iteration= %d \n",xd,it); }

F(X) = A3*X**3 + A2*X**2 + A1*X + A0 READ(5,*)A3,A2,A1,A0,ITM,XI EPS=1.0E-4 DELX=0.1 IT=0 WRITE(6,200)DELX,EPS 200 FORMAT(1X,14HINITIAL DELX=,F10.3,5X,'TOLERANCE=',E10.3) 10 IT=IT+1 FX=F(XI) XD=XI+DELX FXD=F(XD) IF(IT.GT.ITM) GOTO 1000 IF(FX*FXD) 11,14,12 11 IF(DELX.LE.EPS) GOTO 14 DELX=DELX/10.0 GOTO XI=XD GOTO WRITE(6,210)A3,A2,A1,A0,XD,IT 210 FORMAT(1H,F10.3,'X**3',F10.3,'X**2',F10.3,'X',F10.3/ + 1H,6HFIRST=,F10.3/1X,'ITERATION='I3) STOP 1000 WRITE(6,220) IT 220 FORMAT(1H,5X,'ITERATION=',I2,5X,'Check Error') STOP END

1.2 Bisection 법 ( 반분법 ) F 변수구간을 반씩 나누어 근의 존재를 찾음 F [x i, x i+1 ] 을 [x i, x i+1/2 ] [x i+1/2, x i+1 ] 로 구분, x i+1/2 =(x i +x i+1 )/2 0 y x y=f(x) xixi x i+1 f(x i ) f(x i+1 ) F f(x i )*f(x i+1/2 ) < 0 이면 * x i 와 x i+1/2 사이에 근이 존재 * x i+1/2  x i+1 * 또는 f(x i )*f(x i+1/2 ) > 0 이면 x i+1/2  x i * |x i+1 -x i | <  또는 |f(x i+1 )| <  의 조건까지 반복계산

F 계산순서 1. 구간 x i, x i+1 을 정한다. 2. f(x i ) 및 f(x i+1 ) 을 계산한다. 3. x i+1/2 =(x i + x i+1 )/ f(x i )*f(x i+1/2 ) 의 부호를 계산한다. 5.  f(x i )*f(x i+1/2 ) > 0 이면 x i  x i+1/2 입력후 순서 2 로  f(x i )*f(x i+1/2 ) < 0 이면 |x i+1 -x i | <  이면 x i+1/2 를 근으로 계산 끝 |x i+1 -x i |>  이면 x i+1  x i+1/2 입력후 순서 2 로  f(x i )*f(x i+1/2 ) = 0 이면 x i+1/2 를 근으로 계산 끝

/* Program 1-2 */ /* The Bisection Method */ #include #define f(x) (a[3]*pow((x),3)+a[2]*pow((x),2)+a[1]*(x)+a[0]) main() { double a[4]; double low_x, high_x, mid_x, low_y, high_y, mid_y, eps=1.0e-4; int it_max, it=0; printf("Input a[3-0] : "); scanf("%lf%lf%lf%lf", &a[3], &a[2], &a[1], &a[0]); do{ printf("Lower, Upper x Guess : "); scanf("%lf%lf", &low_x, &high_x); }while( f(low_x) * f(high_x) >= 0 ); printf("Input Max Iteration : "); scanf("%d", &it_max);

while(1){ it++; mid_x=(low_x+high_x)/2.0; low_y=f(low_x); high_y=f(high_x); mid_y=f(mid_x); if( it>it_max ){ printf("Iteration= %d Check Error!!!!",it); exit(1); } if( low_y*mid_y>0 )low_x=mid_x; if( low_y*mid_y<0 ) if( fabs(high_x-low_x)<eps )break; else high_x=mid_x; if( low_y*mid_y==0. )break; }; printf("nf(x) = %lfx^3 + %lfx^2 + %lfx + %lf\n",a[3],a[2],a[1],a[0]); printf("Root= %lf \n Iteration= %d \n",mid_x,it); }

F(X)=A3*X**3+A2*X**2+A1*X+A0 READ(5,*) A3,A2,A1,A0,X1,X2,ITM EPS=1.0E-4 IT=0 10 IT=IT+1 FX1=F(X1) FX2=F(X2) X3=(X1+X2)/2.0 FX3=F(X3) IF(IT.GT.ITM) GOTO 1000 IF(FX1*FX3) 11,13,12 11 IF(ABS(X2-X1).LE.EPS) GOTO 13 X2=X3 GOTO X1=X3 GOTO WRITE(6,200) A3,A2,A1,A0,X3,IT 200 FORMAT(1X,F10.3,'*X**3',F10.3,'*X**2',F10.3,'*X',F10.3,/ + 1X,'ROOT=',F10.3/1X,'ITERATION=',I3) STOP 1000 WRITE(6,210) IT 210 FORMAT(1H,'ITERATION =',I2,5X,'Check Error'/) STOP END

1.3 False Position 법 (Regula Falsi 법 ) F 변수구간을 직선식으로 나누어 근의 존재를 찾음 F [x i, x i+1 ] 을 [x i, x i3 ] [x i3, x i+1 ] 로 구분, F x i3 = { x i f(x i+1 )-x i+1 f(x i ) } / { f(x i+1 )-f(x i ) } 0 y x y=f(x) xixi x i+1 f(x i ) f(x i+1 ) F 반분법과 동일한 계산과정 F 근의 수렴속도를 개선한 방법임 F 이 방법은 근의 수렴이 한쪽 방향으로만 이루어지므로 수치해는 근에 무한히 접근하나 일치하지는 않는다. x i3 f(x i3 )

/* Program 1-3 *//* The Method of False Position *//* - Linear Interpolation - */ #include #define f(x) (a[3]*pow((x),3)+a[2]*pow((x),2)+a[1]*x+a[0]) void main(void) { double a[4]; double x1, x2, x3, y1, y2, y3, eps=1.0e-4; int it_max, it=0; printf("Input a3, a2, a1, a0 : "); scanf("%lf%lf%lf%lf", &a[3], &a[2], &a[1], &a[0]); do{ printf("Lower, Upper x Guess : "); scanf("%lf%lf", &x1, &x2); }while( f(x1) * f(x2) >= 0 ); printf("Input Max Iteration : "); scanf("%d", &it_max);

while(1){ it++; y1=f(x1); y2=f(x2); x3=(x1*y2-x2*y1)/(y2-y1); y3=f(x3); if( it>it_max ){ printf("Iteration= %d Check Error!!!!",it); exit(1); } if( y1*y3>0 )x1=x3; if( y1*y3<0 ) if( fabs(x2-x1)<eps )break; else x2=x3; if( y1*y3==0. )break; }; printf("\nf(x) = %5.2lfx^3 + %5.2lfx^2 + %5.2lfx + %5.2lf\n", a[3],a[2],a[1],a[0]); printf("Root= %lf\nIteration= %d \n",x3,it); }

F(X)=2*X**3-11*X**2+2*X+15 READ(5,*)X1,X2,ITM EPS=1.0E-4 IT=0 10 IT=IT+1 FX1=F(X1) FX2=F(X2) X3=(X1*FX2-X2*FX1)/(FX2-FX1) FX3=F(X3) IF(IT.GT.ITM)GOTO 1000 IF(FX1*FX3)11,13,12 11 IF(ABS(X2-X3).LE.EPS)GOTO 13 X2=X3 GOTO IF(ABS(X3-X1).LE.EPS)GOTO 13 X1=X3 GOTO WRITE(6,200)X3,IT 200 FORMAT(1X,'ROOT=',F10.3/1X,'ITERATION=',I3) STOP 1000 WRITE(6,210) IT 210 FORMAT(1X,5X,'ITERATION=',I2,5X,'Check Error') STOP END

1.4 Newton-Raphson 법 F 미분가능한 함수에 대해 적용 F 수렴속도가 매우 빠르다 0 y x y=f(x) x n+1 xnxn (x n, f(x n )) F 곡선위의 점 [x n, f(x n )] 에서의 접선과 x 축과의 교점 [x n+1, 0] F f’(x n )=-f(x n )/{ x n+1 - x n }  x n+1 = x n - f(x n )/f’(x n ) F x n+1  x n 으로 설정 F | x n+1 - x n | <  | (x n+1 - x n )/ x n | <  | f(x n+1 )- f(x n ) | <  의 조건만족때까지 반복계산 

F 계산순서 1. 초기근사해 x n 을 정한다. 2. f(x n ) 및 f’(x n ) 을 계산한다. 3. 다음의 근사해 x n+1 을 구한다. x n+1 = x n - f(x n )/f’(x n ) 4. | x n+1 - x n | <  이면, x n+1 을 근으로 계산 끝 아니면, x n  x n+1 입력후 순서 2 로

/* Program 1-4 *//* Newton-Raphson Method */ #include #define f(x) (pow((x),2)-2*sin(x)) #define df(x) (2*(x)-2*cos(x)) void main(void) {double x, xn, eps=1.0e-4; int it=0, it_max; printf("Input [initial x], [max iteration]:"); scanf("%lf%d", &xn, &it_max); do{it++; x = xn; xn = x - f(x) / df(x); if( it > it_max ){ printf("Iteration= %d Check Error!!!!", it); exit(1); } }while( fabs(xn - x) > eps ); printf("Root= %f \n Iteration= %d \n", xn, it); }

F(X)=X*X-2*SIN(X) DF(X)=2*X-2*COS(X) READ(5,*)XN,ITM EPS=1.0E-4 IT=0 10 IT=IT+1 FX=F(XN) DFX=DF(XN) XNI=XN-FX/DFX IF(IT.GT.ITM) GOTO 1000 IF(ABS(XN-XNI).LE.EPS) GOTO 20 XN=XNI GOTO WRITE(6,200)XNI,IT 200 FORMAT(1X,'ROOT=',F10.5/1X,'ITERATION=',I3) STOP 1000 WRITE(6,210) IT 210 FORMAT(1X,5X,'ITERATION=',I2,5X,'Check Error') STOP END

1.5 Newton’s Second-Order 법 F 2 차 미분 가능한 함수에 대해 적용 F 수렴속도가 매우 빠르다 F 함수 f(x) 를 Taylor 급수전개하면, F 증분이 0 이되면 정확한 해가 얻어지므로, 급수합이 0 이 되도록 한다. (f(x)=0 인 x 값이 근이다.) F 위의 식에서 3 째 전개항까지 고려하고 그 이상을 0 이라고 하면,

F 계산순서 1. 초기근사해 x n 을 정한다. 2. f(x n ) 및 f’(x n ), f’’(x n ) 을 계산한다. 3. 다음의 근사해 x n+1 을 구한다. 4.  | f(x n+1 ) | <  이면, x n+1 을 근으로 계산 끝  아니면, x n  x n+1 입력후 순서 2 로

/* Program 1-5 *//* Newton's Second-Order Method */ #include #define f(x) ( 2*cos(x)-exp(x)) #define df(x) (-2*sin(x)-exp(x)) #define ddf(x) (-2*cos(x)-exp(x)) void main(void) {double x, xn, eps=1.0e-4; int it=0, it_max; printf("\nInput [initial x], [maximum iteration]: "); scanf("%lf%d",&xn,&it_max); do{it++; x = xn; xn = x - f(x)/(df(x) - ddf(x) * f(x)/(2.0 * df(x))); if( it > it_max ){ printf("Iteration= %d Check Error!!!!", it); exit(1); } }while( fabs(f(xn)) > eps ); printf("Root= %lf \nIteration= %d \n", xn, it); }

F(X)=2*COS(X)-EXP(X) D1(X)=-2*SIN(X)-EXP(X) D2(X)=-2*COS(X)-EXP(X) READ(5,*)XN,ITM EPS=1.0E-4 IT=0 10 IT=IT+1 FX=F(XN) D1X=D1(XN) D2X=D2(XN) XNI=XN-FX/(D1X-D2X*FX/(2.0*D1X)) FX=F(XNI) IF(IT.GT.ITM) GOTO 1000 IF(ABS(FX).LE.EPS) GOTO 20 XN=XNI GOTO WRITE(6,200)XNI,IT 200 FORMAT(1X,'ROOT=',F10.3/1X,'ITERATION=',I3) STOP 1000 WRITE(6,210) IT 210 FORMAT(1X,5X,'ITERATION=',I2,5X,'Check Error') STOP END

1.6 Bairstow 법 F 고차방정식을 2 차식의 곱으로 변형하여 근의 공식을 이용 F 실근, 허근을 포함 모든 근을 구할 수 있다. F 위의 방정식을 x 2 +px+q 로 나눈 식이 아래와 같으면, F 나머지를 g 1 x+g 2 라고 하면 다음의 항등식이 성립한다. F 위의 항등식을 원식과 비교하면 오른쪽의 관계식이 얻어진다.

F g 1, g 2 는 p, q 의 함수로서 g 1 (p, q)=0, g 2 (p, q)=0 이 성립하는 p, q 값에 대해서 f(x) 는 (x 2 +px+q) 의 함수를 갖는다.

F 계산순서 1. p 및 q 의 초기값을 정한다. 2. b -1 =0, b 0 =1, a 1,a 2,…,a n 에서 모든 b k 을 계산한다. 3. c -1 =0, c 0 =1, b 1,b 2,…,b n 에서 모든 c k 을 계산한다. 4. b n-1, b n, c n-3, c n-2, c n-1 에서 보정값  p,  q 을 계산한다. 5. p j+1 =p j +  p, q j+1 =q j +  q 4. |  p | <  & |  q | <  이면, x 2 +px+q=0 의 근을 구함 아니면, p j  p j+1, q j  q j+1 입력후 순서 2 로

/* program 1-6 */ /* Bairstow Method */ #include void root(double p, double q); void main(void){ double a[20], b[20], c[20]; double p, q, EPS, det; double xr, xi, delta_p, delta_q; int i, j, k, n, m, it, it_max=10000; scanf("%lf%lf%lf", &p, &q, &EPS); scanf("%d", &n); for(k=1;k<=n+1;k++) scanf("%lf", &a[k]); printf("\n Real part Imaginary part\n"); j = 0; lable1: m = n - 2*j; j++; if(m-2<0){xr = -1*a[2]/a[1]; xi = 0; printf("%12.5lg %12.5lg\n", xr, xi); exit(0);} if(m-2==0){p = a[2]/a[1]; q = a[3]/a[1]; root(p,q); exit(0); } it = 0;

lable2: b[1] = a[1]; b[2] = a[2] - p*b[1]; for(k=3;k<=m+1;k++) b[k] = a[k]-p*b[k-1]-q*b[k-2]; c[1] = b[1]; c[2] = b[2] - p*c[1]; for(k=3;k<=m;k++) c[k] = b[k] - p*c[k-1] - q*c[k-2]; det = pow(c[m-1],2) - c[m-2]*(c[m]-b[m]); if(det==0){ printf("Not converged, Check initial values of p and q\n"); for(k=1;k<=m;k++) printf("a[%d] = %lg\n", k, a[k]); exit(0);} delta_p = ( b[m]*c[m-1] - b[m+1]*c[m-2])/det; delta_q = ( b[m+1]*c[m-1] - b[m]*(c[m]-b[m]))/det; p += delta_p; q += delta_q; it++; if( (fabs(delta_p)<EPS) && (fabs(delta_q)<EPS) ){ root(p,q); for(k=1;k<=m+1;k++) a[k] = b[k]; goto lable1;} if( it<= it_max )goto lable2; printf("Not converged, Check initial values of p and q\n"); for(k=1;k<=m;k++) printf("a[%d] = %lg\n", k, a[k]); }

void root(double p, double q) { double d, xr1, xr2, xi1, xi2; d = pow(p,2) - 4.*q; if(d<0){ xr1 = -1*p/2.; xi1 = sqrt(-1*d)/2.; xr2 = xr1; xi2 = -1*xi1; } else{ xr1 = (-1*p + sqrt(d))/2.; xi1 = 0; xr2 = (-1*p - sqrt(d))/2.; xi2 = 0; } printf("%12.5lg %12.5lg\n", xr1, xi1); printf("%12.5lg %12.5lg\n", xr2, xi2); }

1.7 Bisection 법에 의한 연립비선형방정식의 해석 0 y x f 1 (x,y)=0 f 2 (x,y)=0 x L x M x U yUyMyLyUyMyL f 1 (x L,y L )=0 f 1 (x M,y M )=0 f 1 (x U,y U )=0 f 2 (x,y)=0 f 2 (x L,y L )=? f 2 (x M,y M )=? f 2 (x U,y U )=? 근의 존재조건 f 2 (x L,y L )*f 2 (x M,y M )<0

F 계산순서 1. 근의 존재구역의 상한과 하한 (x U, x L ) 을 정한다. 2. x M = (x U +x L )/2 3. F 1 =0 을 만족하는 근 x, y 를 구한다. 4. f 2 (x L,y L )*f 2 (x M,y M ) 을 계산 부호판정 f 2 (x L,y L )*f 2 (x M,y M ) < 0 ; x U  x M f 2 (x L,y L )*f 2 (x M,y M ) > 0 ; x L  x M f 2 (x L,y L )*f 2 (x M,y M ) = 0 ; (x M, y M ) 근을 구함 5. 수렴판정 | x U - x L | <  이면, (x M, y M ) 근을 구함 아니면, 순서 2 로

/* Program 1-7 */ /* Bisection Method for Two Variables */ #include #define g(x) (sqrt(log((x)*(x)-4.*(x)+10.)+log(2.))) #define f(x,y) ((x)*(x)-(y)*(y)-2.*sin(((x)+(y))/sqrt(2.))) void main(void) {int it_max, it=0; double x_low, x_up, x_mid,y_low, y_mid, f_low, f_mid, eps=1.0E-4; printf("\nInput [Max iteration], [lower x], [upper x]:"); scanf("%d%lf%lf",&it_max,&x_low,&x_up); do{y_low=g(x_low); f_low=f(x_low,y_low); x_mid=(x_low+x_up)/2.0; y_mid=g(x_mid); f_mid=f(x_mid,y_mid); it++; if(it>it_max){ printf("Iteration=%d Check Error!!!",it); exit(1);} if( f_low * f_mid >0 ) x_low=x_mid; if( f_low * f_mid <0 ) x_up=x_mid; }while( fabs(x_up - x_low) > eps && f_low * f_mid != 0 ); printf("The Answer is... \nx= %lf \ny= %lf \nIteration=%d",x_mid,y_mid,it); }

G(X)=SQRT(ALOG(X**2-4.*X+10)+ALOG(2.)) F(X,Y)=X**2-Y**2-2.*SIN((X+Y)/SQRT(2.)) EPS=1.0E-4 READ(5,*)XL,XU YL=G(XL) FL=F(XL,YL) IT=0 10 XM=(XL+XU)/2 YM=G(XM) FM=F(XM,YM) IT=IT+1 IF(IT.GT.20) GOTO 1000 IF(FL*FM.LT.0) THEN XU=XM ELSE XL=XM END IF IF((XU-XL).GT.EPS) GOTO 10 WRITE(6,200) XM,YM STOP 1000 WRITE(6,210) STOP 200 FORMAT(1X,5X,'**SOLUTION**'/ + 1X,5X,'X=',F6.3,5X,'Y=',F6.3) 210 FORMAT(1X,'PROGRAM ERROR') END

1.8 Newton-Raphson 법에 의한 연립비선형방정식의 해석 f(x) 를 x=x (0) 부근에서 Taylor 전개하여 (x-x (0) ) 의 2 차 이상의 항을 무시하면 f(x)=f(x (0) )+G(x (0) )(x-x (0) ) 단 G(x (0) ) 는 n*n 의 정방행렬로서 요소 i,j 의 G ij (x (0) ) 는 에 의해 주어진다. 여기서 f(x)=0, x=x (1) 으로 두면 x (1) = x (0) – G(x (0) ) -1 f(x (0) ) 가 얻어지며, 따라서 반복식은 x (i+1) = x (i) – G(x (i) ) -1 f(x (i) ) 가 된다. 여기서 G(x (i) ) -1 는 G(x (i) ) 의 역행렬이다.

F 계산순서 1. 초기 근사해 벡터 x (0) 를 정한다. 2. f(x (0) ), G(x (0) ) 를 계산한다. 3. 연립방정식의 계수행렬을 만든다. 4. 연립방정식을 풀어 -G(x (0) ) -1 f(x (0) ) 을 구한다. 5. x (1) = x (0) - G(x (0) ) -1 f(x (0) ) 을 구한다. 6. x (i) 의 각 성분값이 허용오차보다 작도록 반복한다.

/* Program 1-8 */ /* Newton-Raphson Method for Two Variables */ #include double x[10], f[10], y[10], g[10][11]; int p, n=2, it_max=20, it=0; double EPSE=1.0E-4, EPSZ; void SRS1(void); void SRS2(void); void SRS3(void); main( ){ int i; printf("\nx[1], x[2] =?"); scanf("%lf%lf",&x[1],&x[2]); printf("\nNUMBER OF EQUATION = %d",n); printf("\nINITIAL VALUES NO XO"); for(i=1;i<=n;i++) printf("\n %d %lf",i,x[i]); printf("\nALLOWABLE = %lf",EPSE); printf("\nMAXIMUN NUMBER OD ITERATION = %d",it_max);

p=0; SRS1(); if(p==1){printf("\nCOUNT OVER!!"); exit(1);} printf("\nNUMBER OF ITERATION = %d",it); printf("\n SOLUTION"); for(i=1;i<=n;i++)printf("\n %d %lf",i,y[i]); } void SRS1(void){ int j; for(it=1;it<=it_max;it++){ SRS2(); SRS3(); for(j=1;j<=n;j++) y[j]= x[j] + g[j][n+1]; for(j=1;j<=n;j++){ EPSZ= g[j][n+1]/y[j]; printf("\nJ=%d EPSZ=%lf EPSE=%lf", j, EPSZ, EPSE); if( fabs(EPSZ)>EPSE )break; } if(j==n)return; for(j=1;j<=n;j++)x[j]=y[j]; } p=1;}

void SRS2(void){ int j; f[1]=pow((x[1]-2),2)+pow((x[2]-2),2)-4.0; f[2]=pow((x[1]-5),2)+pow((x[2]-2),2)-4.0; g[1][1]=2.*(x[1]-2); g[1][2]=2.*(x[2]-2); g[2][1]=2.*(x[1]-5); g[2][2]=2.*(x[2]-2); for(j=1;j<=n;j++)g[j][n+1]=-1.0*f[j];} void SRS3(void){ int l,j,i; double w; for(l=1;l<=n;l++){ w = g[l][l]; for(j=l+1;j<=n+1;j++) g[l][j] = g[l][j]/w; for(i=1;i<=n;i++){ if(i==l)continue; w = g[i][l]; for(j=1;j<=n+1;j++) g[i][j]=g[i][j]-w*g[l][j]; }}

COMMON /C1/ SX(2),SF(2) COMMON /C2/ SY(2) COMMON /C3/ SG(2,3) INTEGER P DATA N/2/ DATA EPSE,M/ ,20/ READ(5,*) (SX(I),I=1,2) WRITE(6,200) N WRITE(6,210) DO 10 I = 1,N WRITE(6,220) I,SX(I) 10 CONTINUE WRITE(6,230) EPSE WRITE(6,240) M C P = 0 CALL SRS1(K,M,N,EPSE,P) IF (P.EQ.1) THEN WRITE(6,280) 280 FORMAT(/10X,'COUNT OVER') STOP ENDIF WRITE(6,250) K WRITE(6,260) DO 20 I=1,N WRITE(6,270) I,SY(I) 20 CONTINUE 200 FORMAT(/10X,'NUMBER OF EQUATION =',I2) 210 FORMAT(/10X,'INITIAL VALUES',6X,'NO',8X,'XO') 220 FORMAT(/30X,I2,8X,3F3.1) 230 FORMAT(/10X,'ALLOWABLE =',E12.5) 240 FORMAT(/10X,'MAXIMUM NUMBER OF ITERATIONS =',I3) 250 FORMAT(/10X,'NUMBER OF ITERATION =',I2) 260 FORMAT(/25X,' SOLUTION') 270 FORMAT(/22X,I2,3X,2F8.1/) STOP END

SUBROUTINE SRS1(K,M,N,EPSE,P) COMMON /C1/ SX(2),SF(2) COMMON /C2/ SY(2) COMMON /C3/ SG(2,3) DO 10 K=1,M CALL SRS2(N) CALL SRS3(N) DO 20 J=1,N SY(J)=SX(J) + SG(J,N+1) 20 CONTINUE DO 30 J = 1,N EPSZ = SG(J,N+1)/SY(J) IF (ABS(EPSZ).GT.EPSE) GO TO CONTINUE GO TO DO 40 J=1,N SX(J) = SY(J) 40 CONTINUE 10 CONTINUE P = RETURN END C SUBROUTINE SRS2(N) COMMON /C1/ SX(2),SF(2) COMMON /C3/ SG(2,3) R1 = 2.0 R2 = 2.0 X1 = 2.0 Y1 = 2.0 X2 = 5.0 Y2 = 2.0 A1 = R1**2 A2 = R2**2 SF(1) =(SX(1) - X1)**2 + (SX(2)-Y1)**2 - A1 SF(2) =(SX(1) - X2)**2 + (SX(2)-Y2)**2 - A2 C SG(1,1) =2*(SX(1)-X1) SG(1,2) =2*(SX(2)-Y1) SG(2,1) =2*(SX(1)-X2) SG(2,2) =2*(SX(2)-Y2) DO 10 J = 1,N SG(J,N+1) = -SF(J) 10 CONTINUE RETURN END C

SUBROUTINE SRS3(N) COMMON /C3/ SG(2,3) DO 10 L=1,N SW = SG(L,L) DO 20 J=(L+1),N+1 SG(L,J) = SG(L,J)/SW 20 CONTINUE DO 30 I=1,N IF(I.EQ.L) GO TO 30 SW = SG(I,L) DO 40 J=1,N+1 SG(I,J) = SG(I,J) - SW*SG(L,J) 40 CONTINUE 30 CONTINUE 10 CONTINUE RETURN END

2 장 연립 1 차 방정식 ( System of Linear Simulatenous Equations ) 1.Gauss 소거법 2.Gauss-Jordan 법 3.LU 분해법 4.Thomas 법 5.Gauss-Seidel 법 6.SOR 법 7. 역행렬 (Inverse Matrix) 직접법 간접법