Part II. 편미분 방정식의 해법
1장 유한차분법
1.1 서론
편미분 방정식은 다음 방정식과 같이 독립변수의 수가 2개 이상인 경우이다.
(1.1)
간편하게 미분식을 표현하기 위하여 1계, 2계 미분을 다음과 같이 나타낸다.
(1.2)
편미분 방정식의 일반적인 형태는 다음과 같다.
(1.3)
편미분 방정식 계수는 방정식내 가장 계수가 높은 미분항에 의해 결정된다.
변수가 여러개의 독립변수에 의존하는 경우 (), 미분값은 다음과 같은 방향별 미분을 고려한 전미분의 식으로 나타내어야 한다.
(1.4)
편미분 방정식은 다음과 같이 계수의 특성에 따라 선형, 반선형, 비선형으로 나누어진다.
(1.5)
선형 : a, b, c가 상수이거나 독립변수의 함수인 경우
반선형 : 계수가 가장 높은 미분항 선형인 경우
비선형 : 계수가 종속 변수의 함수인 경우
편미분 방정식은 특성식에 따라 쌍곡형(hyperbolic), 포물형(parabolic), 타원형(elliptic)으로 나누어지며, 이에 따라 수치해의 불안정성과 오차가 발생하는 특성이 다르기 때문에 서로 다른 수치해석 기법을 적용하여야 한다. 특성식과 유한차분법을 이해하기 위해서는 다음과 같이 차분 변수에 기초를 둔 유한 수학과 편미분의 기본 개념에 대한 이해를 필요로 한다.
1.2 유한 수학(finite mathematics)의 개념
차분법을 이해하기 위해서는 미분과 차분의 개념을 정확히 이해하여야 한다. 수치해석에서는 미분식을 차분화된 식으로 해석하는 유한 수학의 원리를 이용한다. 미분식과 차분식의 차이점은 다음의 표에서 알 수 있다. 이 표에서 알 수 있드시 차분식에서는 유한의 개념을 사용한다. 즉, 컴퓨터는 주어진 자료만을 처리하는 계산기이기 때문에 무한이 수를 계산할 수는 없고, 무수히 많은 수를 계산함으로서 근사치를 구하게 된다.
표 1.1 미분식과 차분식의 기본 개념의 차이
미분식
차분식
차분식은 미분식의 근사치를 구하는 것이며 그 오차는 다음과 같은 테일러 급수로서 평가된다.
(1.6)
따라서, 위의 표에 제시된 차분식으로 미분식을 평가하는 경우 오차는 다음과 같다.
(1.7)
이러한 오차를 절단오차(truncation error)라고 하며, 오차항 중에서 가장 계수가 큰 항이 1계이기 때문에 1계정확(fisrt order correct) 절단오차라 한다.
테일러 급수는 전개 방식에 따라 전방(forward) 테일러 급수, 후방(backward) 테일러 급수 등이 있다. 전후방 테일러 급수를 합치면 다음과 같이 2계 도함수의 차분식을 얻을 수 있으며 오차는 2계 정확 절단오차가 된다.
(1.8)
전방 테일러 급수에서 후방 테일러 급수를 빼주면, 다음과 같은 1계도함수의 차분식이 구해지는 데, 오차는 2계 정확 절단오차로서 표보다 오차가 줄어든다.
(1.9)
공학 및 자연과학 문제의 대부분의 편미분방정식에서는 2계이상의 도함수는 사용되지 않으므로 다음과 같은 1계 및 2계 도함수에 대한 차분식을 가지고 유한차분법을 수행할 수 있다.
표 1.2 유한차분법에서 많이 사용하는 1계, 2계 도함수의 차분식
1계 정확 (first-order correct)
2계 정확 (second-order correct)
2계 정확 (second-order correct)
수치해석 기법의 정확도 및 계산의 효율에 따라 차분식의 유형을 선택하게 된다.
1.3 편미분 방정식의 분류
편미분 방정식은 특성식에 따라 쌍곡형(hyperbolic), 포물형(parabolic), 타원형(elliptic)으로 나누어진다. 이러한 이름은 해석기학학의 분류 원리에 따라 정의되었다. 편미분방정식의 유형에 따라 적절한 수치해석 기법을 선택하여야 한다.
1.3.1 1계 준선형(1st order quasi-linear) 편미분 연립 방정식
다음의 준선형 1계 연립 편미분 방정식의 유형을 평가하자.
(1.10)
(1.11)
위의 준선형방정식이 유지되기 위해서는 모든 계수는 독립변수의 도함수()에 의존하면 안된다. 만약, u와 v가 연속적으로 미분가능하다면, 다음과 같은 방향도함수를 고려할 수 있다.
(1.12)
(1.13)
위의 4개의 식을 이용하여 다음과 같이 4개의 미지수에 대한 연립방정식을 구성할 수 있다.
(1.14)
위의 편미분 연립 방정식의 행렬식(determinant)이 0이 되어야만 무수히 많은 해가 존재하며, 경계조건 및 초기조건에 따라 주어진 문제의 해를 결정하게 된다. 만약에 행렬식이 0이 아니면 Crammer의 법칙에 따라 해를 다음과 같이 구하게 된다.
다음의 연립방정식의 해는 다음과 같다.
(1.15)
(1.16)
여거서, 행렬식은 다음과 같다.
, (1.17)
만약, 이면 의 의미가 없는 해(trivial solution)를 갖고,이면 무수히 많은 해을 갖는다. 이 무수히 많은 해중 경계조건과 초기 조건을 만족시키는 해를 구한다.
행렬식은 다음과 같이 계산된다.
(1.18)
(1.19)
(1.20)
위식의 항들을 정리하면 다음의 식을 얻을 수 있다.
(1.21)
dx와 dy의 어떤 특성이 행렬식을 0으로 만드는지 결정하기 위하여 위의 식을 dx의 자승으로 나누어 주면 다음의 식을 얻을 수 있다.
(1.22)
위 식의 해는 다음과 같이 구해진다.
(1.23)
(1.24)
여기서,
에 따라 실수나 복소수의 해를 얻는다. 해에 따라 다음과 같이 편미분 방정식을 분류한다.
인 경우 2개의 실근을 가지며 쌍곡형 편미분 방정식이 된다.
인 경우 1개의 실근을 가지며 포물형 편미분 방정식이 된다.
인 경우 2개의 허근을 가지며 타원형 편미분 방정식이 된다.
위 특성식의 해 나 를 특성방향이라 칭한다. 편미분방정식의 해는 이러한 방향을 따라 진행된다.
1.3.2 2계 준선형(2nd order quasi-linear) 편미분 방정식
다음과 같은 2계 준선형 편미분 방정식의 유형을 평가하자.
(1.25)
본 문제는 2계 미분 문제로서 다음과 같이 2개의 1계 미분 문제로 변환시킬 수 있다.
방향별 미분은 다음과 같다.
(1.26)
(1.27)
행렬 형태로 연립방정식을 구성하면 다음과 같다.
(1.28)
행렬식은 다음과 같다.
(1.29)
혹은 다음과 같이 표현할 수 있다.
(1.30)
A = a, B = -b, C = c라 하면
인 경우 2개의 실근을 가지며 쌍곡형 편미분 방정식이 된다.
인 경우 1개의 실근을 가지며 포물형 편미분 방정식이 된다.
인 경우 2개의 허근을 가지며 타원형 편미분 방정식이 된다.
1.3.3 선형 파동(1inar wave) 편미분 방정식
다음의 선형 파동 문제의 편미분 방정식의 분류는 다음과 같다.
(1.31)
, , 이므로 이다. 값은 항상 양수이므로 에 편미분 방정식 항상 쌍곡형의 유형이다. 시간에 따라 값이 변하는 경우 특성식의 해가 시간에 따라 달라지므로 수치해법에 있어서 이러한 특성선의 성질을 고려하여야 한다. 이러한 경우는 특성선을 추적하여야 하므로 입자추적법의나 Lagrangian 수치해법을 사용하게 된다.
1.4 포물형 편미분방정식의 해석
일반적으로 지하수 흐름의 지하수 유동 방정식 및 오염물 이동에 관계된 물질이동방정식은 포물형의 편미분 방정식 유형에 속한다. 이러한 일차원 물질이동방정식에 대하여 여러 수치해석 기법을 적용하여 일차원 모형을 개발하였다. Excel을 사용하여 모형을 개발함으로서 수치해석 알고리즘에 대하여 계산결과를 상세히 확인할 수 있었다. 이러한 계산결과는 저자가 기존에 개발한 BASIC 및 FORTRAN을 이용한 수치해석모형과 비교하여 수치해석 기법의 타당성을 검증하였다. 또한 모든 수치해석 기법의 결과는 이론해와 비교하였다. 수치해석 알고리즘으로 유한차분법 방법중 범용적 Crank-Nicholson 해법 (GCN : Generalized Crank-Nicholson Method)과 범용적 특성 평균법 (GCA : Generalized Characteristic Averaging Method)를 사용하였다. 이중 GCN 모형은 다음과 같이 서술할 수 있다.
1.4.1 범용적 Crank-Nicholson 해법
(H.W. and Exam. 2009) 물질이동식을 해석하기 위한 범용적 유한차분법 알고리즘을 유도하라.
1차원 물질이동식은 다음과 같다.
(1.32)
유한차분법에서는 위의 물질이동식을 수치해석적으로 해석하기 위하여 모든 도함수를 차분화한다. Taylor 급수를 이용하여 차분화하게 되는 데, 시간영역을 차분화하는 방법에 따라 전방차분법, 후방차분법, Crank-Nicholson법 등으로 구분된다. 범용적 Crank-Nicholson법에서는 시간가중함수를 이용하여 이러한 시간 차분 알고리즘을 일반화한 것이다. 시간가중함수 이 0인 경우에는 전방차분법, 1이면 후방차분법, 0.5이면 중앙차분법 혹은 Crank-Nicholson법이 된다. GCN의 수치해석 아날로그는 다음과 같다. 유속에 대한 항을 차분화하기 위하여 중앙차분법을 이용하였다.
(1.33)
위의 식을 간략히 표현하기 위하여 다음의 계수를 정의한다.
, , (1.34)
시간차분항을 정리한다.
(1.35)
위의 식을 주변수에 대하여 정리한다.
(1.36)
위의 식을 모든 절점에 대하여 합쳐주면 삼대각 행렬식이 구성된다.
(1.37)
(1.38)
(1.39)
삼대각 행렬은 다음과 같다.
(1.40)
이러한 삼대각 행렬이 구성되는 이유는 각 절점에 연결된 2개의 절점만 수치해석 아날로그에 포함되기 때문에 합쳐줄 때 각 절점별 식은 3개의 절점(확산항을 평가할 때 필요)을 필요로 하기 때문이다. 이러한 삼대각 행렬은 가우스소거법의 일종인 Thomas Algorithm을 이용하여 해를 구한다.
1.4.2 경계조건의 해석
경계조건에는 경계조건이 주변수로 주어지는 제1경계조건, 주변수의 도함수로 주어지는 제2경계조건, 주변수의 도함수가 경계지역 내외부의 주변수로부터 결정되는 제3경계조건이 있다. 주변수가 농도인 경우에는 다음과 같이 해석된다.
o 제1경계조건
제1경계조건에서는 좌측() 혹은 우측농도()로 주어진다. 이러한 주어진 농도값으로 해를 고정시키기 위하여 최종행렬식에서 좌측절점(1) 및 우측절점(NS)에 대한 식의 성분은 다음과 같이 결정된다.
좌측 경계조건인 경우 : (1.41)
우측 경계조건인 경우 : (1.42)
o 제2경계조건
제2경계조건에서는 다음과 같이 주어진 플럭스로 경계조건이 결정된다.
(*는 플럭스를 의미한다.) (1.43)
좌측경계조건식은 다음과 같이 해석된다.
(1.44)
최종행렬식의 각 성분은 다음과 같이 해석된다.
(1.45)
(1.46)
(1.47)
따라서, 을 에 을 에 더해준다.
(1.48)
우측경계조건은 다음과 같은 식으로 주어진다.
(1.49)
따라서, 다음과 같이 최종행렬식을 성분을 해석한다.
(1.50)
(1.51)
(1.52)
(1.53)
를 에, 에 를 더해준다.
(1.54)
o 제3경계조건
제3경계조건은 다음의 식으로 정의된다.
(1.55)
좌측경계조건은 다음의 식으로 정의된다.
(1.56)
여기서, 은 플럭스의 전도도 계수이고, 는 경계지역에서의 주변수의 값이다.
위의 식은 다음과 같이 전개된다.
(1.57)
(1.58)
따라서, 다음과 같이 더해준다.
(1.59)
(1.60)
우측경계조건은 다음의 식으로 정의된다.
(1.61)
따라서, 우측 절점에 대한 최종행렬식은 다음과 같다.
(1.62)
(1.63)
따라서, 다음과 같이 더해준다.
(1.64)
(1.65)
1.4.3 물질이동식의 이론해(Analytical Solution)
1) 오차함수(Error Function)를 이용한 이론해
일차원 물질이동식의 이론해는 경계조건에 따라서 Boltzman 변환을 사용하거나 Laplace 변환을 사용하여 다음의 초기조건 및 경계조건에 대하여 이론해를 얻을 수 있다.
초기조건:,
경계조건: ,
물질이동식에 Laplace 변환을 적용하면 다음과 같은 이론해를 얻을 수 있다.
(1.66)
여기서, , 이다. (1.67)
위의 오차함수는 다음과 같이 급수를 이용하여 수치해석적으로 구할 수 있다. (1.68)
오차함수를 해석하기 위한 프로그램은 다음과 같다.
--------------- F1.1 오차함수를 해석하기 위한 프로그램 -----------------
character*10 answer
write(*,*) 'stop computation?(y/n)'
read(*,*) answer
if (answer='y') then
goto 100
endif
10 write (*,*) 'input z, D, t, theta0, thetai='
read (*,*) z, D, t, theta0, theti
u = z / (4 * D * t) ** 0.5
call err(u,sumu)
subroutine err(u, sumu) (sub err)
sumu = 0
ifact = 1
do i = 2, 100 (for i =1 2 to 100)
ifact=ifact*(i-1)
sumu = sumu + (-1)**(i+1)*u**(2*i-1)/((2*i-1)*ifact)
enddo
sumu = sumu + u
sumu = sumu * (4 / 3.1415917) ** 0.5
return
end
theta = theta0 + (thetai - theta0) * sumu
write(*,*) 'soil moisture content =', theta
goto 10
100 stop
end
-------------------------------------------------------------------
2) Fourier 급수를 사용한 일반해
다음과 같은 Fourier 급수를 사용한 일반해를 가정한다.
(1.69)
여기서, 는 시간영역에서의 n번째 요소의 진동수, 는 공간영역에서의 진동수이다. 다음과 같이 n번째 해를 고려하여 물질이동식에 대입한다.
, (1.70)
다음과 같이 진동수 간의 관계를 구할 수 있다.
, (1.71)
Fourier 급수의 해의 n번째 해는 다음과 같다.
(1.72)
여기서, Fourier 급수의 각 성분에 있어서 는 이동을 나타내고, 는 크기의 변화를 나타낸다. 따라서, ∆t 시간 이수에는 n번째 진동파는 만큼 감소하고 만큼 이동할 것이다.
1.4.4 Basic을 이용한 GCN 모형의 개발
저자는 Quick Basic Compiler를 이용하여 다음과 같은 GCN 모형을 개발하였다. 본 모형은 Graphic 모듈을 내장하여 계산결과를 도시할 수 있는 장점이 있다. 다음에 본 프로그램의 주 모듈을 나타내었다.
--------- B1.1 GCN 알고리즘을 이용한 유한차분모형 --------------------
10 '
20 ' Program 2.BAS by Professor Joon Hyun Kim
30 '
40 ' Generalized Crank-Nicholson Method
50 '
60 ' 1 Dimensional Mass Transport Problem
70 '
80 ' CE284H, Winter 1987, Revision Winter 2007
90 '
100 '
110 CLS : SCREEN 0: KEY OFF: 'WIDTH "lpt1:", 130: ' LPRINT CHR$(15)
120 '
130 DIM X(100), Y(100), U(100), OLDU(100), UANAL(100)
132 DIM A(100), B(100), C(100), D(100), BETA(100), GAMMA(100), VEC(100)
140 '
150 ' INPUT "PRINT FILE, HOW MANY DATA SETS TO TRY="; FILE1$, DSET
152 DSET = 3
154 FILE1$ = "3.PL1": FILE2$ = "3.PL2": FILE3$ = "3.PL3"
160 OPEN FILE1$ FOR OUTPUT AS #1
170 OPEN FILE3$ FOR OUTPUT AS #3
180 '
190 KOUNT2 = 0
200 GOSUB 1000: 'Input Data
210 '
220 GOSUB 2000: 'A(I),B(I),C(I) of Tridiagonal Banded Matrix with B.C.
230 '
240 KOUNT1 = 0: TIME = 0: ITER = 0
250 TIME = TIME + DT: ITER = ITER + 1
260 IF TIME > TMAX THEN 360
270 '
280 GOSUB 3000: 'R.H.S LOAD VECTOR D(I) & B.C. & SOLVE MATRIX FOR U(I)
290 '
292 GOSUB 7000: 'COMPUTE ANALYTICAL SOLUTION UANAL(I)
294 '
300 IP = INT(ITER / IPRINT) * IPRINT
310 IF IP = ITER THEN GOSUB 4000: 'Store the result for output file
320 '
330 FOR I = 1 TO NS: OLDU(I) = U(I): NEXT I
340 GOTO 250: 'ITERATATION
350 '
360 TIME = TIME - DT: GOSUB 4000: CLOSE #2: 'print last time step
370 '
380 OPEN FILE2$ FOR INPUT AS #2
390 KOUNT = KOUNT1: INF = 2
400 ' GOSUB 5000: ' <- Plot every modeling results
410 CLOSE #2
420 '
422 '... store the simulation result at last time step into file #3 ...
444 '
430 KOUNT2 = KOUNT2 + 1
440 FOR I = 1 TO NS: PRINT #3, USING " ##.#### "; U(I); : NEXT I: PRINT #3, " "
445 KOUNT2 = KOUNT2 + 1
448 FOR I = 1 TO NS: PRINT #3, USING " ##.#### "; UANAL(I); : NEXT I: PRINT #3, " "
450 IF KOUNT2 = 2 * DSET THEN 460 ELSE 200
460 CLOSE #3
470 OPEN FILE3$ FOR INPUT AS #3
480 KOUNT = KOUNT2: INF = 3
490 GOSUB 5000: ' <- Plot the last time result
500 CLOSE #3
510 '
520 END
1000 '
1010 '========= 1. Input Data ====================================
1020 '
1030 OPEN FILE2$ FOR OUTPUT AS #2
1040 INPUT "DD="; DD
1050 IPRINT = 50: DT = .025: NS = 41: EP = .5
1060 TMAX = 1.25: RK = 0: V = .369
1070 DX = 1 / (NS - 1): U0 = 0: YMIN = 0: YMAX = 1.2
1080 IBL = 1: UBL = 1: IBR = 2: UBR = 0: KBL = 0: KBR = 0
1090 '
1100 FOR I = 1 TO NS: OLDU(I) = U0: NEXT I
1110 OLDU(1) = UBL
1120 FOR I = 1 TO NS: X(I) = DX * (I - 1): NEXT I
1130 XMIN = X(1): XMAX = X(NS)
1140 '
1150 PRINT #1, "INPUT DATA"
1160 PRINT #1, "NS="; NS; " dx="; DX; " dt="; DT; " tmax="; TMAX; " IPRINT ="; IPRINT
1170 PRINT #1, "D="; DD; " V="; V; " k="; RK; " ep="; EP; " U0 ="; U0
1180 PRINT #1, "IBL,UBL,IBR,UBR,KBL,KBR="; IBL; UBL; IBR; UBR; KBL; KBR
1190 '
1200 RETURN
2000 '
2010 '====== 2. L.H.S. Coefficient Matrix A(I), B(I), C(I) & B.C. ============
2020 '
2030 AL = DD * DT / DX ^ 2: BE = .5 * V * DT / DX: GA = RK * DT
2040 FOR I = 1 TO NS
2050 A(I) = -EP * (AL + BE)
2060 B(I) = 1 + EP * (2 * AL + GA)
2070 C(I) = -EP * (AL - BE)
2080 NEXT I
2090 '
2100 ' Left Boundary Conditions
2110 '
2120 IF IBL = 1 THEN 2150
2130 IF IBL = 2 THEN 2180
2140 IF IBL = 3 THEN 2190
2150 '
2160 B(1) = 1: C(1) = 0: GOTO 2210
2170 '
2180 C(1) = C(1) + A(1): D(1) = D(1) + A(1) * 2 * DX * UBL: GOTO 3200
2190 '
2200 B(1) = B(1) - 2 * DX * KBL * A(1): C(1) = C(1) + A(1): GOTO 3200
2210 '
2220 ' Right Boundary Contions
2230 '
2240 IF IBL = 1 THEN 2270
2250 IF IBL = 2 THEN 2290
2260 IF IBL = 3 THEN 2310
2270 '
2280 A(NS) = 0: B(NS) = 1: GOTO 2340
2290 '
2300 A(NS) = A(NS) + C(NS): GOTO 2340
2310 '
2320 A(NS) = A(NS) + C(NS): B(NS) = B(NS) + 2 * DX * C(NS) * KBR: GOTO 2340
2330 '
2340 PRINT #1, "DD*dt/dx^2, V*dt/2/dx, k*dt=", USING " ###.#### ###.#### ###.####"; AL; BE; GA
2350 PRINT #1, "a(1),b(1),c(1)=", USING "###.#### ###.#### ###.####"; A(1); B(1); C(1)
2360 '
2370 PE = V * DX / DD: CR = 2 * BE: CRPE = CR / PE
2380 PRINT #1, "Peclet No.="; PE; " Courant No.="; CR; " Cr/Pe="; CRPE: PRINT #1, " "
2390 RETURN
3000 '
3010 '========== 3. R.H.S. LOAD VECTOR & B.C. ====================
3020 '
3030 FOR I = 2 TO NS - 1
3040 D(I) = (1 - EP) * (AL + BE) * OLDU(I - 1) + (1 - (1 - EP) * (2 * AL + GA)) * OLDU(I) + (1 - EP) * (AL - BE) * OLDU(I + 1)
3050 NEXT I
3060 '
3070 ' Left Boundary Conditions
3080 '
3090 IF IBL = 1 THEN 3120
3100 IF IBL = 2 THEN 3150
3110 IF IBL = 3 THEN 3180
3120 '
3130 D(1) = UBL: GOTO 3200
3140 '
3150 D(1) = (1 - EP) * (AL + BE) * (OLDU(2) - 2 * DX * UBL) + (1 - (1 - EP) * (2 * AL + GA)) * OLDU(1) + (1 - EP) * (AL - BE) * OLDU
3160 D(1) = D(1) + A(1) * 2 * DX * UBL: GOTO 3200
3170 '
3180 D(1) = (1 - EP) * (AL + BE) * (OLDU(2) - 2 * DX * KBL * (OLDU(1) - UBL)) + (1 - (1 - EP) * (2 * AL + GA)) * OLDU(1) + (1 - EP) * (AL - BE) * OLDU(2)
3190 D(1) = D(1) - A(1) * 2 * DX * KBL * U(1): GOTO 3200
3200 '
3210 ' RIGHT BOUNDARY CONDITION
3220 '
3230 IF IBR = 1 THEN 3260
3240 IF IBR = 2 THEN 3290
3250 IF IBR = 3 THEN 3320
3260 '
3270 D(NS) = UBR: GOTO 3340
3280 '
3290 D(NS) = (1 - EP) * (AL + BE) * OLDU(NS - 1) + (1 - (1 - EP) * (2 * AL + GA)) * OLDU(NS) + (1 - EP) * (AL - BE) * (OLDU(NS - 1) + 2 * DX * UBR)
3300 D(NS) = D(NS) - C / (NS) * 2 * DX * UBR: GOTO 3340
3310 '
3320 D(NS) = (1 - EP) * (AL + BE) * OLDU(NS - 1) + (1 - (1 - EP) * (2 * AL + GA)) * OLDU(NS) + (1 - EP) * (AL - BE) * (OLDU(NS - 1) + 2 * DX * KBR * (OLDU(NS) - UBR))
3330 D(NS) = D(NS) + 2 * DX * KBR * UBR: GOTO 3340
3340 '
3350 ' SOLVE TRIDIAGONAL MATRIX
3360 N = NS
3370 GOSUB 6000
3380 FOR I = 1 TO NS: U(I) = VEC(I): NEXT I
3390 '
3400 RETURN
4000 '
4010 ' ============== 4. PRINT OUT RESULT =======================
4020 '
4030 '
4040 PRINT #1, " ": PRINT #1, "TIME="; USING "##.####"; TIME
4050 KPRINT = 0: 'Print out with the format of 11 columns
4060 FOR I = 1 TO 11
4070 KPRINT = KPRINT + 1
4080 IF KPRINT > NS THEN 4130
4090 PRINT #1, USING " ##.####"; U(KPRINT);
4100 PRINT #2, USING " ##.####"; U(KPRINT);
4110 NEXT I
4120 PRINT #1, " ": GOTO 4060
4130 PRINT #1, " ": PRINT #2, " ": KPRINT = 0: KOUNT1 = KOUNT1 + 1
4140 FOR I = 1 TO 11
4150 KPRINT = KPRINT + 1
4160 IF KPRINT > NS THEN 4190
4165 PRINT #1, USING " ##.#### "; UANAL(KPRINT);
4170 PRINT #2, USING " ##.#### "; UANAL(KPRINT);
4185 NEXT I
4180 PRINT #1, " ": GOTO 4140
4190 PRINT #1, " ": PRINT #2, " ": KOUNT1 = KOUNT1 + 1
4200 '
4210 RETURN
5000 '
5010 ' ============5. PLOT THE RESULT ========================
5020 '
5030 KK = 0
5040 KK = KK + 1: PRINT #1, "KK,KOUNT="; KK; KOUNT
5050 IF KK > KOUNT THEN 5230
5060 FOR I = 1 TO NS: INPUT #INF, Y(I): PRINT #1, "I,Y="; I; Y(I): NEXT I
5070 '
5080 SCREEN 2
5090 VIEW (90, 10)-(634, 164)
5100 WINDOW (XMIN, YMIN)-(XMAX, YMAX)
5110 LOCATE 1, 12: PRINT " D="; DD; " V="; V; " DX="; DX; " DT="; DT; " EP="; EP
5120 LOCATE 2, 4: PRINT USING "##.##"; YMAX: LOCATE 21, 4: PRINT USING "##.##"; YMIN
5130 LOCATE 22, 10: PRINT USING " ##.##"; XMIN: LOCATE 22, 74: PRINT USING "##.##"; XMAX
5140 LOCATE 23, 20: PRINT "Distance in X direction at time="; USING "##.###"; TIME
5150 LINE (XMIN, YMIN)-(XMAX, YMAX), , B
5160 PSET (X(1), Y(1))
5170 FOR I = 2 TO NS
5180 LINE -(X(I), Y(I))
5190 NEXT I
5200 A$ = INPUT$(1)
5210 GOTO 5040
5220 '
5230 A$ = INPUT$(1): 'A$ = INPUT$(1)
5240 'SCREEN 0
5250 '
5260 RETURN
6000 '
6010 ' ======================6. Thomas Algorithm ====================
6020 '
6030 '
6040 ' Find Beta And Gamma
6050 '
6060 BETA(1) = B(1)
6070 GAMMA(1) = D(1) / B(1)
6080 FOR I = 2 TO N
6090 IM = I - 1
6100 BETA(I) = B(I) - A(I) * C(IM) / BETA(IM)
6110 GAMMA(I) = (D(I) - A(I) * GAMMA(IM)) / BETA(I)
6120 NEXT I
6130 '
6140 ' Solve For Vector
6150 '
6160 VEC(N) = GAMMA(N)
6170 NM = N - 1
6180 FOR I = 1 TO NM
6190 NMI = N - I
6200 NMI1 = NMI + 1
6210 VEC(NMI) = GAMMA(NMI) - C(NMI) * VEC(NMI1) / BETA(NMI)
6212 ' PRINT #1, "NMI,VEC(NMI)="; NMI, VEC(NMI)
6220 NEXT I
6230 '
6240 RETURN
7000 '
7010 '====================7. ANALYTIC SOLUTION===================
7020 '
7030 UANAL(1) = UBL
7040 DTT = DD * TIME
7050 DTT1 = SQR(DTT)
7060 FOR I = 2 TO NS
7070 UUU = .5 * (X(I) - V * TIME) / DTT1
7080 GOSUB 8000: '<-ERROR FUNCTION
7090 UANAL(I) = .5 * UBL * REC
7100 NEXT I
7110 '
7120 RETURN
8000 '
8010 '====================8. ERROR FUCTION========================
8020 '
8030 SUMC = 0
8040 BB = 1: UUU1 = UUU
8050 IF UUU < 0 THEN UUU = -UUU
8060 IF UUU > 2! THEN 8070 ELSE 8080
8070 RE = 1: GOTO 8140
8080 FOR KKK = 1 TO 10
8090 BB = BB * KKK
8100 CC = (-1) ^ KKK * UUU ^ (2 * KKK + 1) / (BB * (2 * KKK + 1))
8110 SUMC = CC + SUMC
8120 NEXT KKK
8130 RE = 1.12838 * (UUU + SUMC)
8140 IF UUU1 < 0 THEN RE = -RE
8150 REC = 1 - RE
8160 '
8170 RETURN
-------------------------------------------------------------------
위의 프로그램의 수치해석적 안정성 및 정확도를 분석하기 위하여, 확산이 지배적인 경우(D=0.01725, V=0.369), 유속이 지배적인 경우(D=0.0001725), 중간 경우(D=0.001725)에 대하여 모델링을 수행하여 다음과 같이 비교하였다.
그림 1.1 GCN-Basic 모형의 계산 결과
1.4.5 Excel을 이용한 GCN 모형의 개발
Excel로서 GCN 모형을 개발하였다. 그림에 나타난 바와 같이 Basic 모형인 GCN.BAS 프로그램의 결과와 동일하다.
그림 1.2 GCN-Excel 모형의 계산 과정
그림 1.3 GCN-Excel 모형의 계산 결과
1.4.6 범용적 특성 평균 해법 (GCA)
GCA방법에 의한 알고리즘은 다음과 같이 특성법(Characteristic Method), 중앙차분평균법(Centered Difference (Averaging) Method), 범용적(Generalized) Crank Nicholson 법의 결합에 의해 유도된다.
1) 특성법(Characteristic Method)
편미분방정식은 특성식의 속성에 따라 다음과 같이 분류된다.
타원형 : 정상상태의 확산식 : (1.73)
포물선형 : 물질이동방정식 : (1.74)
쌍곡선형 : 연속방정식 : (1.75)
위의 편미분방정식중 쌍곡선형인 경우에는 수치해석 불안정성과 오차가 심하기 때문에 격자망을 작게 사용하거나, 특성법을 사용하여야 한다. 특성법은 다음과 같이 설명된다.
연속방정식에 전방차분법을 적용하면 다음과 같다.
(1.76)
위의 전방차분 아날로그에 대하여 특성법에서는 주변수가 유속을 따라 이동한다고 가정한다. 즉. 시간과 공간격자의 크기를 특성선상의 유속에 따라 다음과 같이 결정한다.
(1.77)
따라서, 다음의 식으로 전개된다.
, (1.78)
즉, 주변수가 격자망을 따라서 유속에 의하여 이동하는 것을 알 수 있다.
물질이동식에 GCA의 알고리즘을 적용하기 위해서는 유속항은 중앙차분법을, 확산항은 특성법을 적용한다. 즉 분산이 유속에 따라 이동된 격자에서 일어난 것으로 가정한다. 따라서, 분산항은 다음과 같이 평가된다.
n 시각 : (1.79)
n+1 시각 : (1.80)
2) 중앙차분(평균)법 (Centered Difference (Averaging) Method)
물질이동식중 유속에 의한 항은 중앙차분법으로 평가된다. 즉, 시간도함수는 공간절점에 대하여 평균치를 취하고, 공간도함수는 시간절점에 대하여 평균치를 취한다. (1.81)
(1.82)
위의 식을 정리하면 다음과 같다.
(1.83)
여기서, , 이다.
만약 평균을 취하지 않는 다면, 위의 식은 다음과 같다 (전방차분특성법).
(1.84)
인 경우, 두 평균 차분 특성법이나 전방 차분 특성법 모두 다음과 같이 유속에 따라 주변수를 추적하는 알고리즘이 되므로 특성법을 입자추적법이라고도 한다.
(1.85)
3) GCN을 이용한 확산항의 평가
GCA에서는 확산항을 평가하기 위해서 다음과 같은 격자망을 사용하므로, 확산항에 대한 해석은 다음과 같다.
그림 1.4 GCA의 시간 및 공간의 격자점
o 시간에 대하여 평균화된 GCN
(1.86)
o 시간에 대하여 가중계수를 사용한 GCN
(1.87)
4) 평균특성법과 GCN의 결합
평균특성법 GCN은 다음과 같이 시간미분에 대하여 평균치를 사용하거나 가중계수를 사용하여 평가될 수 있다.
o 평균 GCA :
(1.88)
o 가중 GCA :
(1.89)
위의 GCA 아날로그는 파라미터를 추출하기 위하여 다음과 같이 정리된다.
o 평균 GCA
=
- (1.90)
o 가중 GCA
=
- (1.91)
계수는 GCN 모형과 유사하게 정의된다.
, (1.92)
, , (1.93)
위의 계수를 이용하여 식을 정리하면 다음과 같다.
o 평균 GCA
=
- (1.94)
o 가중 GCA
=
- (1.95)
위의 식을 현재와 장래의 주변수에 대해서 정리한다.
o 평균 GCA
(1.96)
, , (1.97)
(1.98)
o 가중 GCA
(1.99)
, , (1.100)
(1.101)
5) 최종행렬식의 구성
GCN 방법과 마찬가지로 최종행렬식은 다음과 같은 각 절점에 대한 식을 합침으로서 구성된다. 아래 식의 모든 절점을 조합하면 좌측행렬이 삼대각 행렬이 된다.
(1.102)
6) Basic을 이용한 GCA 모형의 개발
저자는 Quick Basic Compiler를 이용하여 다음과 같은 GCA 모형을 개발하였다. 본 모형은 Graphic 모듈을 내장하여 계산결과를 도시할 수 있는 장점이 있다.
---------- B1.2 GCA 알고리즘을 이용한 유한차분모형 --------------------
10 '
20 ' Program GCA.BAS by Professor Joon Hyun Kim
30 '
40 ' Characteristic Averaging Method
50 '
60 ' 1 Dimensional Mass Transport Problem
70 '
80 ' CE284H, Winter 1987, Revision Winter 2007
90 '
100 '
110 CLS : SCREEN 0: KEY OFF: 'WIDTH "lpt1:", 130: ' LPRINT CHR$(15)
120 '
130 DIM X(100), Y(100), U(100), OLDU(100), UANAL(100)
132 DIM A(100), B(100), C(100), D(100), BETA(100), GAMMA(100), VEC(100)
140 '
142 ' ...input output file name and nubmer of data sets...
144 '
150 ' INPUT "OUTPUT FILE, HOW MANY DATA SETS TO TRY="; FILE1$, DSET
152 DSET = 3
154 FILE1$ = "7.PL1": FILE2$ = "7.PL2": FILE3$ = "7.PL3"
160 OPEN FILE1$ FOR OUTPUT AS #1
170 OPEN FILE3$ FOR OUTPUT AS #3
180 '
190 KOUNT2 = 0
200 GOSUB 1000: 'Input Data
210 '
220 GOSUB 2000: 'A(I),B(I),C(I) of Tridiagonal Banded Matrix with B.C.
230 '
240 KOUNT1 = 0: TIME = 0: ITER = 0
250 TIME = TIME + DT: ITER = ITER + 1
260 IF TIME > TMAX THEN 360
270 '
280 GOSUB 3000: 'R.H.S LOAD VECTOR D(I) & B.C. & SOLVE MATRIX FOR U(I)
290 '
292 GOSUB 7000: 'COMPUTE ANALYTICAL SOLUTION UANAL(I)
294 '
196 ' ...store the solution at input interval...
198 '
300 IP = INT(ITER / IPRINT) * IPRINT
310 IF IP = ITER THEN GOSUB 4000
320 '
330 FOR I = 1 TO NS: OLDU(I) = U(I): NEXT I
340 GOTO 250: 'ITERATATION
350 '
352 ' ...store the solution at last time for time <= tmax to the file #2...
354 '
360 TIME = TIME - DT: GOSUB 4000: CLOSE #2
370 '
372 ' ...plot the result of file #2 (at each time step)...
374 '
380 OPEN FILE2$ FOR INPUT AS #2
390 KOUNT = KOUNT1: INF = 2
400 ' GOSUB 5000: ' <-Plot
410 CLOSE #2
420 '
422 ' ...comparison plotting of numerical and analytical solutions of file #3...
424 '
430 KOUNT2 = KOUNT2 + 1
440 FOR I = 1 TO NS: PRINT #3, USING " ##.#### "; U(I); : NEXT I
445 KOUNT2 = KOUNT2 + 1
448 FOR I = 1 TO NS: PRINT #3, USING " ##.#### "; UANAL(I); : NEXT I
450 IF KOUNT2 = 2 * DSET THEN 460 ELSE 200
460 CLOSE #3
470 OPEN FILE3$ FOR INPUT AS #3
480 KOUNT = KOUNT2: INF = 3
490 GOSUB 5000: ' <-Plot
500 CLOSE #3
510 '
520 END
1000 '
1010 '============ 1. Input Data =============================
1020 '
1030 OPEN FILE2$ FOR OUTPUT AS #2
1040 INPUT "DD="; DD
1050 NS = 41: TMAX = 1.25: EP = .5: DX = 1 / (NS - 1): '... simulation control data ...
1060 U0 = 0: RK = 0: V = .369: DT = DX / V: '... physical data ...
1070 IBL = 1: UBL = 1: IBR = 2: UBR = 0: KBL = 0: KBR = 0: '... boundary condition ...
1080 IPRINT = 50: YMIN = 0: YMAX = 1.2: '... printing and plotting control ...
1090 '
1100 FOR I = 1 TO NS: OLDU(I) = U0: NEXT I: '...initial concentration...
1110 OLDU(1) = UBL: '...boundary condition...
1120 FOR I = 1 TO NS: X(I) = DX * (I - 1): NEXT I: '...nodal points...
1130 XMIN = X(1): XMAX = X(NS): '...min and max value of x-axis...
1140 '
1150 PRINT #1, "INPUT DATA"
1160 PRINT #1, "NS="; NS; " dx="; DX; " dt="; DT; " tmax="; TMAX; " IPRINT ="; IPRINT
1170 PRINT #1, "D="; DD; " V="; V; " k="; RK; " ep="; EP; " U0 ="; U0
1180 PRINT #1, "IBL,UBL,IBR,UBR,KBL,KBR="; IBL; UBL; IBR; UBR; KBL; KBR
1190 '
1200 RETURN
2000 '
2010 '==== 2. L.H.S. COMPONENTS OF TRIDIAGONAL BANDED MATRIX A(I), B(I), C(I) & B.C. ======
2020 '
2030 AL = DD * DT / DX ^ 2: BE = .5 * V * DT / DX: GA = RK * DT
2032 PE = V * DX / DD: CR = 2 * BE: CRPE = CR / PE
2034 PRINT #1, "Peclet No.="; PE; " Courant No.="; CR; " Cr/Pe="; CRPE: PRINT #1, " "
2036 D1 = 1 / PE: D2 = 1 + 2 * BE - 2 / PE - .5 * GA: D3 = 1 - 2 * BE + 1 / PE - .5 * GA
2038 '
2040 FOR I = 1 TO NS
2050 A(I) = -1 / PE + .5 * GA
2060 B(I) = 1 + 2 * BE + 2 / PE + .5 * GA
2070 C(I) = -1 / PE
2080 NEXT I
2090 '
2100 ' Left Boundary Conditions
2110 '
2120 IF IBL = 1 THEN 2150
2130 IF IBL = 2 THEN 2180
2140 IF IBL = 3 THEN 2190
2150 '
2155 A(2) = -.5 * (AL + BE): B(2) = 1 + .5 * (2 * AL + GA): C(2) = -.5 * (AL - BE)
2160 B(1) = 1: C(1) = 0: GOTO 2210
2170 '
2180 C(1) = C(1) + A(1): D(1) = D(1) + A(1) * 2 * DX * UBL: GOTO 3200
2190 '
2200 B(1) = B(1) - 2 * DX * KBL * A(1): C(1) = C(1) + A(1): GOTO 2210
2210 '
2220 ' Right Boundary Contions
2230 '
2240 IF IBL = 1 THEN 2270
2250 IF IBL = 2 THEN 2290
2260 IF IBL = 3 THEN 2310
2270 '
2280 A(NS) = 0: B(NS) = 1: GOTO 2340
2290 '
2300 A(NS) = A(NS) + C(NS): GOTO 2340
2310 '
2320 A(NS) = A(NS) + C(NS): B(NS) = B(NS) + 2 * DX * C(NS) * KBR: GOTO 2340
2330 '
2340 PRINT #1, "DD*dt/dx^2, V*dt/2/dx, k*dt=", USING " ###.#### ###.#### ###.####"; AL; BE; GA
2350 PRINT #1, "a(1),b(1),c(1)=", USING "###.#### ###.#### ###.####"; A(1); B(1); C(1)
2360 '
2390 RETURN
3000 '
3010 '=========== 3. R.H.S. LOAD VECTOR & B.C. ================
3020 '
3030 FOR I = 3 TO NS
3040 D(I) = D1 * OLDU(I - 2) + D2 * OLDU(I - 1) + D3 * OLDU(I)
3050 NEXT I
3055 D(2) = -A(2) * OLDU(1) + (2 - B(2)) * OLDU(2) - C(2) * OLDU(3)
3060 '
3070 ' LEFT BOUNDARY CONDITION
3080 '
3090 IF IBL = 1 THEN 3120
3100 IF IBL = 2 THEN 3150
3110 IF IBL = 3 THEN 3180
3120 '
3130 D(1) = UBL: GOTO 3200
3140 '
3150 D(1) = D1 * (OLDU(2) - 2 * DX * UBL) + D2 * OLDU(1) + D3 * OLDU(2)
3160 D(1) = D(1) + A(1) * 2 * DX * UBL: GOTO 3200
3170 '
3180 D(1) = D1 * (OLDU(2) - 2 * DX * KBL * (OLDU(1) - UBL)) + D2 * OLDU(1) + D3 * OLDU(2)
3190 D(1) = D(1) - A(1) * 2 * DX * KBL * U(1): GOTO 3200
3200 '
3210 ' RIGHT BOUNDARY CONDITION
3220 '
3230 IF IBR = 1 THEN 3260
3240 IF IBR = 2 THEN 3290
3250 IF IBR = 3 THEN 3320
3260 '
3270 D(NS) = UBR: GOTO 3340
3280 '
3290 D(NS) = D1 * OLDU(NS - 1) + D2 * OLDU(NS) + D3 * (OLDU(NS - 1) + 2 * DX * UBR)
3300 D(NS) = D(NS) - C / (NS) * 2 * DX * UBR: GOTO 3340
3310 '
3320 D(NS) = D1 * OLDU(NS - 1) + D2 * OLDU(NS) + D3 * (OLDU(NS - 1) + 2 * DX * KBR * (OLDU(NS) - UBR))
3330 D(NS) = D(NS) + 2 * DX * KBR * UBR: GOTO 3340
3340 '
3350 ' SOLVE TRIDIAGONAL MATRIX
3360 N = NS
3370 GOSUB 6000
3380 FOR I = 1 TO NS: U(I) = VEC(I): NEXT I
3390 '
3400 RETURN
4000 '
4010 ' == 4. STORE THE RESULT IN FILE #1 AND #2 FOR PLOTTING ===
4020 '
4030 '
4040 PRINT #1, " ": PRINT #1, "TIME="; USING "##.####"; TIME
4050 KPRINT = 0: 'Print out with the format of 11 columns
4060 FOR I = 1 TO 11
4070 KPRINT = KPRINT + 1
4080 IF KPRINT > NS THEN 4130
4090 PRINT #1, USING " ##.####"; U(KPRINT);
4100 PRINT #2, USING " ##.####"; U(KPRINT);
4110 NEXT I
4120 PRINT #1, " ": GOTO 4060
4130 PRINT #1, " ": PRINT #2, " ": KPRINT = 0: KOUNT1 = KOUNT1 + 1
4140 FOR I = 1 TO 11
4150 KPRINT = KPRINT + 1
4160 IF KPRINT > NS THEN 4190
4165 PRINT #1, USING " ##.#### "; UANAL(KPRINT);
4170 PRINT #2, USING " ##.#### "; UANAL(KPRINT);
4185 NEXT I
4180 PRINT #1, " ": GOTO 4140
4190 PRINT #1, " ": PRINT #2, " ": KOUNT1 = KOUNT1 + 1
4200 '
4210 RETURN
5000 '
5010 ' =============5. PLOT THE RESULT =================
5020 '
5030 KK = 0
5040 KK = KK + 1: PRINT #1, "KK,KOUNT="; KK; KOUNT
5050 IF KK > KOUNT THEN 5230
5060 FOR I = 1 TO NS: INPUT #INF, Y(I): PRINT #1, "I,Y="; I; Y(I): NEXT I
5070 '
5080 SCREEN 2
5090 VIEW (90, 10)-(634, 164)
5100 WINDOW (XMIN, YMIN)-(XMAX, YMAX)
5110 LOCATE 1, 12: PRINT " D="; DD; " V="; V; " DX="; DX; " DT="; DT; " EP="; EP
5120 LOCATE 2, 4: PRINT USING "##.##"; XMAX: LOCATE 21, 4: PRINT USING "##.##"; IMIN
5130 LOCATE 22, 10: PRINT USING " ##.##"; XMIN: LOCATE 22, 74: PRINT USING "##.##"; XMAX
5140 LOCATE 23, 20: PRINT "Distance in X direction at time="; USING "##.###"; TIME
5150 LINE (XMIN, YMIN)-(XMAX, YMAX), , B
5160 PSET (X(1), Y(1))
5170 FOR I = 2 TO NS
5180 LINE -(X(I), Y(I))
5190 NEXT I
5200 A$ = INPUT$(1)
5210 GOTO 5040
5220 '
5230 A$ = INPUT$(1): 'A$ = INPUT$(1)
5240 'SCREEN 0
5250 '
5260 RETURN
6000 '
6010 ' ======6. THOMAS ALGORITHM FOR THE SOLUTION OF TRIDIAGONAL BANDED MATRIX =============
6020 '
6030 '
6040 ' Find Beta And Gamma
6050 '
6060 BETA(1) = B(1)
6070 GAMMA(1) = D(1) / B(1)
6080 FOR I = 2 TO N
6090 IM = I - 1
6100 BETA(I) = B(I) - A(I) * C(IM) / BETA(IM)
6110 GAMMA(I) = (D(I) - A(I) * GAMMA(IM)) / BETA(I)
6120 NEXT I
6130 '
6140 ' Solve For Vector
6150 '
6160 VEC(N) = GAMMA(N)
6170 NM = N - 1
6180 FOR I = 1 TO NM
6190 NMI = N - I
6200 NMI1 = NMI + 1
6210 VEC(NMI) = GAMMA(NMI) - C(NMI) * VEC(NMI1) / BETA(NMI)
6212 ' PRINT #1, "NMI,VEC(NMI)="; NMI, VEC(NMI)
6220 NEXT I
6230 '
6240 RETURN
7000 '
7010 '==================7. ANALYTIC SOLUTION================
7020 '
7030 UANAL(1) = UBL
7040 DTT = DD * TIME
7050 DTT1 = SQR(DTT)
7060 FOR I = 2 TO NS
7070 UUU = .5 * (X(I) - V * TIME) / DTT1
7080 GOSUB 8000: '<-ERROR FUNCTION
7090 UANAL(I) = .5 * UBL * REC
7100 NEXT I
7110 '
7120 RETURN
8000 '
8010 ' =====================8. ERROR FUCTION======================
8020 '
8030 SUMC = 0
8040 BB = 1: UUU1 = UUU
8050 IF UUU < 0 THEN UUU = -UUU
8060 IF UUU > 2! THEN 8070 ELSE 8080
8070 RE = 1: GOTO 8140
8080 FOR KKK = 1 TO 10
8090 BB = BB * KKK
8100 CC = (-1) ^ KKK * UUU ^ (2 * KKK + 1) / (BB * (2 * KKK + 1))
8110 SUMC = CC + SUMC
8120 NEXT KKK
8130 RE = 1.12838 * (UUU + SUMC)
8140 IF UUU1 < 0 THEN RE = -RE
8150 REC = 1 - RE
8160 '
8170 RETURN
-------------------------------------------------------------------
GCN 모형과 마찬가지로수치해석상의 안정도 및 정확도를 고찰하기 위하여 다음과 같이 확산이 지배적인 경우(D=0.01725, V=0.369), 유속이 지배적인 경우(D=0.0001725), 중간 경우(D=0.001725)에 대하여 모델링을 수행하여 다음과 같이 비교하였다.
그림 1.5 GCA-Basic 모형의 계산 결과
7) Excel을 이용한 GCA 모형의 개발
Excel로서 GCA 모형을 개발하였다. 그림에 나타난 바와 같이 GCA.BAS 프로그램의 결과와 동일하다.
그림 1.6 GCA-Excel 모형의 계산 과정
그림 1.7 GCA-Excel 모형의 계산 결과