1 / 98

Numerical Analysis Spring Semester 2011 Homework (part 1)

Numerical Analysis Spring Semester 2011 Homework (part 1). 환경공학과 20071477 오호식. I.1 Taylor 급수를 유도하고 절단 오차를 설명하라. 절단오차.

leo-murphy
Télécharger la présentation

Numerical Analysis Spring Semester 2011 Homework (part 1)

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript

  1. Numerical Analysis Spring Semester 2011Homework (part 1) 환경공학과 20071477 오호식

  2. I.1 Taylor 급수를 유도하고 절단 오차를 설명하라.

  3. 절단오차 수학적으로 엄밀하게 주어지는 함수f 의 값을, 유한의 사칙 연산의 반복 계산 식fa 로 근사하는 경우, (fa-f)의 오차가 생긴다. 이것을, 절단 오차(truncation error)라 한다. 예를 들면, 테일러급수를 이용하여 삼각함수를 계산할 때, 무한급수의 계산을 유한 항까지의 계산으로 중단하기 위해서, 절단오차가 발생한다.

  4. I.2 다음의 역행렬을 구하는 문제에 대하여 1) 알고리즘을 설명하고 프로그램을 작성 실행하라. 2) 책의 프로그램을 Visual Fortran, Visual C 등의 Compiler를 이용하여 실행하고 알고리즘 및 계산 결과를 설명하라. 3) 본 프로그램을 Visual Basic 프로그램으로 변환하라.

  5. (1)알고리즘 • 알고리즘:알고리즘: 어떠한 문제를 해결하기 위한 여러 동작들의 유한한 모임이다. • 입력 : 외부에서 제공되는 자료가 0개 이상 존재한다. • 출력 : 적어도 1개 이상의 결과를 내어야 한다. • 명확성 : 각 명령어들은 명확하고 모호하지 않아야 한다. • 유한성 : 알고리듬의 명령어들은 유한번의 수행후에 종료되어야 한다. 이것 은 수행 시간의 현실적인 유한성을 의미한다. • 효과성 : 모든 명령어들은 원칙적으로 종이와 연필만으로 수행될 수 있는 기본적인 것이어야 한다.

  6. 2) 책의 프로그램을 Visual Fortran, Visual C 등의 Compiler를 이용하여 실행하고 알고리즘 및 계산 결과를 설명하라. Visual Fortran <Visual Fortran> c < 주 프로그램 > implicit double precision (a-h, o-z) parameter (nmax=8) dimension a(nmax,nmax), wv(nmax), lsw(nmax) c *** 단계 1 자료 입력 *** write(6,*) 'LU 분해에 따른 역행렬의 계산' write(6,*) 'LU 분해의 계산' write(6,*) '행렬의 크기 n =' read(5,*) n write(6,*) '행렬 요소의 입력' call s_inmg(a,n,n,'a')

  7. c *** 단계 2 X=A^-1의 계산 *** call s_invm(a,wv,n,lsw,ierr) c *** 단계 3 계산 결과의 표시 *** if(ierr.eq.1) then write(6,*) '행렬이 정칙이 아님' else write(6,*) '계산 결과' call s_outm(a,n,n,'x') end if stop end c < LU 분해에 따른 역행렬 계산의 subroutine > subroutine s_invm(a,wv,n,lsw,ierr) c a : n*n 입력 행렬 /계산 결과 (in/out) c wv : n 워크 벡터 (work) c n : 스칼라 행렬의 크기 (in) c lsw : n*1 행의 교환 정보 (work) c ierr : 스칼라 에러 코드 0:정상 1:비 정상 (out) implicit double precision (a-h, o-z) dimension a(n,n), wv(n), lsw(n)

  8. c *** 단계 1 Doolittle법에 따른 LU 분해 *** call s_dot(a,n,lsw,ierr) if(ierr.eq.1) then return end if c *** 단계 2 y의 반복 계산 *** do i=1,n do j=i+1,n sum=0.0 do k=i+1,j-1 sum=sum+a(j,k)*a(k,i) enddo a(j,i)=-a(j,i)-sum enddo enddo

  9. c *** 단계 3 x의 반복 계산 *** do j=n,1,-1 do k=j,n wv(k)=a(j,k) enddo do i=1,n sum=0.0 do k=j+1,n sum=sum+wv(k)*a(k,i) enddo if(j.lt.i) a(j,i)=-sum/wv(j) if(j.eq.i) a(j,i)=(1.0-sum)/wv(j) if(j.gt.i) a(j,i)=(a(j,i)-sum)/wv(j) enddo enddo c *** 단계 4 pivot 교환에 따른 열의 교환 *** do i=1,n do j=1,n wv(j)=a(i,j) enddo do k=1,n a(i,lsw(k))=wv(k) enddo enddo ierr=0 return end

  10. <Visual C> /* < 주 프로그램 > */ #include "smp.h" main() { int i,N,*LSW; double *a, *wv; /*** 단계 1 자료 입력 ***/ printf("LU 분해에 따른 역행열의 계산\n"; printf("행렬 A의 크기 N = "); scanf("%d", &N); a=calloc(N*N, sizeof (double)); wv=calloc(N, sizeof (double)); LSW=calloc(N, sizeof (int)); if(a==NULL || wv==NULL || LSW==NULL){ printf("\n 메모리가 확보되지 않았음"); exit(1); } printf("행렬 A의 요소의 입력\n"); S_INMG(a,N,N,"a");

  11. /*** 단계 2 X=A^-1의 계산 ***/ i=S_INVM(a,wv,N,LSW); /*** 단계 3 계산 결과의 표시 ***/ if(i){ printf("행렬이 정칙이 아님") }else{ printf("계산 결과\n"); S_OUTM(a,N,N,"x") } } /* < LU 분해에 따른 역행열 계산의 함수 > */ int S_INVM(double *a, double *wv, int N, int *LSW) /* 반환값 : 에러 코드 0:정상 1:이상 a : N*N 입력 행렬 /역행 열 (In/Out) wv : N*1 작업 용 (Work) N : 스칼라 행렬의 크기 (In) LSW : N*1 행의 교환 정보 (Work) */ { int i, j, k, ks; double sum;

  12. /*** 단계 1 Doolittle법에 따른 LU 분해 ***/ i=S_DOT(a,N,LSW) if(i==1) return 1; /*** 단계 2 Y의 반복 계산 ***/ for(i=0;i<N;i++){ for(j=i+1;j<N;j++){ sum=0.0; for(k=i+1;k<j;k++) sum=sum+a[N*j+k]*a[N*k+i]; a[N*j+i]=-a[N*j+i]-sum; } } /*** 단계 3 X의 반복 계산 ***/ for(j=N-1;j>=0;j--){ for(k=j;k<N;k++) wv[k]=a[N*j+k]; for(i=0;i<N;i++){ sum=0; for(k=j+1;k<N;k++) sum=sum+wv[k]*a[N*k+i]; if(j<i) a[N*j+i]=-sum/wv[j]; if(j==i) a[N*j+i]=(1.0-sum)/wv[j]; if(j>i) a[N*j+i]=(a[N*j+i]-sum)/wv[j]; } }

  13. /*** 단계 4 Pivot 교환에 따른 열의 교환 ***/ • for(i=0;i<N;i++){ • for(j=0;j<N;j++) wv[j]=a[N*i+j]; • for(k=0;k<N;k++) a[N*i+LSW[k]]=wv[k]; • } • return 0; • }

  14. III.1 3원 연립방정식에 대하여 행렬식(Determinant)을 이용하여 역행렬을 구하는 방법과 연립방정식을 푸는 방법을 설명하라 (수학책 참조).

  15. 행렬식 • 정방 행렬A의 LU분해가 A=LU인 경우, 그 행렬식|A| 는 |A|=|L| |U| • 로 표현할 수 있다. Doolittle법에 따라LU 분해에서는, L의 대각 요소는 모두 1이고, 삼각 행렬의 행렬식은 그 대각 요소의 적(積)이기 때문에 • 가 된다. 따라서, |A| 은 U의 모든 대각 요소의 적으로서 구해진다. 단,LU 분해에 있어서 Pivot의 선택이 있는 경우는, Pivot을 교환할 때마다 행렬식의 부호를 바꾸어주기 위해, 교환 회수가 m인 경우

  16. 가 된다. 따라서, |A| 은 U의 모든 대각 요소의 적으로서 구해진다. 단,LU분해에 있어서 Pivot의 선택이 있는 경우는, Pivot을 교환할 때마다 행렬식의 부호를 바꾸어주기 위해, 교환 회수가 m인 경우 • 가 된다. • Doolittle법에 따른LU분해의 결과, 구해진U와 Gauss의 소거법에 있어서 전진 소거가 끝난 단계에서의 식 (4.14)의 계수 행렬이 같은 것은 앞서 서술하였다. 이것으로부터, Gauss 소거법에 있어서도 모든 Pivot의 적

  17. 에 따라 행렬식을 구하는 것이 가능하다. 행렬식만 구하는 경우는, Gauss의 단순 소거법에 있어서 우변의 계산과 후진 대입의 계산을 생략하고, 전진 소거를 하여, 모든 Pivot의 적을 구하면 된다. 따라서, 행렬식의 계산 순서는 다음과 같다.

  18. 역행렬 • 행렬 A의 역행렬 이란, 다음의 관계를 만족시키는 행렬이다. • 여기에서,I는 단위 행렬이다. 역행렬의 계산은, 식(4.4)로부터 구하는 것보다, 소거법 등을 이용하는 편이 계산 횟수가 적게 든다. 여기에서는,LU분해에 따른 계산 방법을 서술한다.

  19. 은 미지이므로 X라고 치환하면, 식 (5.5)로부터 AX=I • (5.6)으로 나타내진다. Doolittle법의 LU 분해에 따라 A=LU를 구하여 LY=I • 의 관계로부터 행렬 를 계산하고, 계속해서 UX=Y

  20. 의 관계로부터 행렬X를 계산할 수 있다. 여기에서, 는 LY=I 의 관계로부터 대각 요소가 모두 1인 하삼각행렬이 되고, 의 순서로, 그리고 j를 1부터 n까지라고 하면

  21. 으로 구해진다. 또, 는 의 순서로, 그리고 j를 n부터 1까지라고 하면 • 으로 구해진다. 또한, Y는 하삼각행렬이기 때문에 L의 배열에 중복되며, X도 U의 영역에 중복되어, 이에 따라 배열 영역을 절약할 수 있다. 나아가, LU분해의 결과를 A에 중복시키면, X의 계산 결과도 A에 중복시키는 것이 가능하다. 단, LU분해에 있어서 Pivot을 교환하는 경우는, 그것에 대응하는 X의 열교환이 필요하게 된다. Pivot 교환이 있는 경우의 계산 순서를 서술하였다.

  22. III.2 3원 연립방정식을 이용하여 Gauss 소거법의 알고리즘을 유도하라.

  23. III.3 다음의 3원 연립방정식의 해를 Gauss의 단순 소거법 및 Pivot 선택법을 이용하여 구하라. (1) (2) (3)

  24. III.4 위의 3원 연립방정식의 해를 교과서에서 제시된 C 및 Fortran 프로그램을 운영하여 구하라.

  25. <Fortran프로그램> c < 주 프로그램 > implicit double precision (a-h, o-z) integer ierr parameter (max=8) dimension x(max), a(max, max), b(max) c *** 단계 1 자료 입력 *** write(6,*) 'Gauss의 소거법에 따른 ax=b (a:n*n, x, b:n*1)의 계산‘ write(6,*) 'n=' read(5,*) n call s_inmg(a,n,n,'a') call s_inmg(b,n,1,'b') c *** 단계 2 Gauss 소거법의 계산 *** call s_gaus(x,a,b,n,ierr) c *** 단계 3 계산 결과의 표시 *** if(ierr.eq.1) then write(6,*) '부정확 또는 불능' stop end if call s_outm(x,n,1,'x') stop end

  26. c < Gauss 소거법의 subroutine > subroutine s_gaus(x,a,b,n,ierr) c c x : n * 1 방정식의 해 (out) c a : n * n 입력 행렬 (in/work) c b : n * 1 입력 벡터 (in/work) c n : scalar 방정식의 원래 해 (in) c ierr : scalar 에러 코드 0 : 정상 1 : 이상 (out) c implicit double precision (a-h, o-z) dimension x(n), a(n,n), b(n) parameter (e_eps=1.0e-14)

  27. c *** 단계 1 전진 소거 *** • do k=1,n-1 • call s_pivo(a,lc,n,k,ierr) • if(ierr.eq.1) then • return • else if(lc.ne.k) then • tmp = b(k) • b(k) = b(lc) • b(lc) = tmp • end if • do i=k+1,n • p=a(i,k)/a(k,k) • do j=k+1,n • a(i,j)=a(i,j)-p*a(k,j) • enddo • b(i)=b(i)-p*b(k) • enddo • enddo • c *** 단계 2 후진 대입 *** • do k=n,1,-1 • if(abs(a(k,k)).le.e_eps) then • ierr=1 • return • end if • s=0.0 • do j=k+1,n • s=s+a(k,j)*x(j) • enddo • x(k)=(b(k)-s)/a(k,k) • enddo • ierr=0 • return • end

  28. c < pivot 선택(행의 교환)의 subroutine > • subroutine s_pivo(a,m,n,k,ierr) • c a : n * n pivot 교환을 하는 행렬 (out/in) • c m : 스칼라 교체된 번호 (out) • c n : 스칼라 행렬의 크기 (in) • c k : 스칼라 pivot의 선택의 열의 위치 (in) • c ierr : 스칼라 에러 코드 0 : 정상 1 : 이상 (out) • implicit double precision (a-h, o-z) • dimension a(n,n) • parameter (e_eps=1.0e-14) • c *** 단계 1 초기 설정 *** • m = k • d = abs(a(k,k)) • c *** 단계 2 최대 pivot의 탐색 *** • do i=k+1,n • if(abs(a(i,k)).gt.d) then • m = i • d = abs(a(i,k)) • end if • enddo

  29. c *** 단계 3 이상 처리 *** • if(abs(d).le.e_eps) then • ierr=1 • return • else if(m.eq.k) then • ierr=0 • return • end if • c *** 단계 4 pivot의 교환 *** • do i=k,n • tmp = a(k,i) • a(k,i) = a(m,i) • a(m,i) = tmp • enddo • ierr=0 • return • end

  30. c < 행렬의 key 입력 subroutine > • subroutine s_inmg(a,n,m,name) • c a : n * m key 입력을 하는 행렬 (out) • c n : 스칼라 행의 크기 (in) • c m : 스칼라 열의 크기 (in) • c name : 1문자 행렬의 이름 (in) • implicit double precision (a-h, o-z) • character name*1 • dimension a(n,m) • write(6,*) • do i=1,n • if(m.eq.1) then • write(6,30) name, i • 30 format(a,'(',i2,')의 입력') • else • write(6,10) name, i, name, i, m • 10 format(a,'(',i2,', 1)~',a,'(',i2,',',i2,')의 입력') • end if • read(5,*)(a(i,j),j=1,m) • enddo • return • end

  31. c < 행렬의 화면 표시 subroutine > • subroutine s_outm(a,n,m,name) • c a : n * m 화면 표시를 하는 행렬 (out) • c n : 스칼라 행의 크기 (in) • c m : 스칼라 열의 크기 (in) • c name : 1 문자 행렬의 이름 (in) • implicit double precision (a-h, o-z) • character*1 name • dimension a(n,m) • write(6,10) name • 10 format(1h,a,'=') • do i=1,n • write(6,20)(a(i,j),j=1,m) • 20 format(8f15.6) • enddo • return • end

  32. <C 프로그램> • #include<stdio.h> • #include<conio.h> • void main() • { • int i,j,k,h; • float l,m,n,o,p,q,r; • clrscr(); • float A[3][4]={0, }; • for(i=0;i<3;i++) • { • for(j=0;j<4;j++) • { • printf("Input number A[%d][%d]=",i,j); • scanf("%f",&A[i][j]); • } • }

  33. for(j=3;j>=0;j--) • { • A[0][j]=A[0][j]/A[0][0]; • } • l=A[1][0]/A[0][0]; • m=A[2][0]/A[0][0]; • for(k=1;k<2;k++) • { • for(j=0;j<4;j++) • { • A[k][j]=A[k][j]-A[0][j]*l; • A[k+1][j]=A[k+1][j]-A[0][j]*m; • } • } • for(j=3;j>0;j--) • { • A[1][j]=A[1][j]/A[1][1]; • } • o=A[0][1]/A[1][1]; • p=A[2][1]/A[1][1]; • for(k=0;k<1;k++)

  34. { • for(j=1;j<4;j++) • { • A[k][j]=A[k][j]-(A[1][j]*o); • A[k+2][j]=A[k+2][j]-(A[1][j]*p); • } • } • for(j=3;j>1;j--) • { • A[2][j]=A[2][j]/A[2][2]; • } • q=A[0][2]/A[2][2]; • r=A[1][2]/A[2][2]; • for(k=0;k<1;k++) • { • for(j=2;j<4;j++) • { • A[k][j]=A[k][j]-(A[2][j]*q); • A[k+1][j]=A[k+1][j]-(A[2][j]*r); • } • } • for(i=0;i<3;i++) • { • for(j=0;j<4;j++) • { • printf(" %5.1f",A[i][j]); • } • printf("\n"); • } • }

  35. III.6 4원 연립방정식을 이용하여 LU 분해법의 알고리즘을 Doolittle 방법 및 Crout 방법에 대하여 유도하라.

  36. III.7 위의 3원 연립방정식의 해를 LU 분해법으로 계산기와 Excel을 이용하여 구하라.

More Related