글
신변잡기/공대딩 2009. 7. 1. 20:16수치해석 MATLAB 열전달 Gauss Seidel, Jacobi, LU Decomposition
결과그래프
프로그램 코드
LU Decomposition.
Lu Decomposition 프로그램코드
clc;
syms('AA', 'NAME', 'INP', 'OK', 'N', 'I', 'J', 'A');
syms('FLAG', 'XL', 'X1', 'M', 'KK', 'S', 'K', 'JJ');
syms('SS', 'OUP', 's');
syms('DATAHX','DATAHY','DATAK','DATAQ','DATAX','DATAY','DATAS','DATAF','DATAH');
syms('NORX','NORY','TWX','TWY','DELTAX','DELTAY','NORQ');
TRUE = 1;
FALSE = 0;
fprintf(1,'LU분해 해법.\n');
fprintf(1,'주어진 h를 입력하시오.\n');
DATAHX = input(' ');
DATAHY = DATAHX;
fprintf(1,'주어진 k를 입력하시오.\n');
DATAK = input(' ');
fprintf(1,'주어진 q를 입력하시오.\n');
DATAQ = input(' ');
fprintf(1,'주어진 유체의 온도 T_fluid를 입력하시오.\n');
DATAF = input(' ');
fprintf(1,'바의 x방향 길이를 입력하시오.(0.005)\n');
DATAX = input(' ');
fprintf(1,'바의 y방향 길이를 입력하시오.(0.02)\n');
DATAY = input(' ');
fprintf(1,'미소면적의 개수 (n X n)일 때, n를 입력하시오.(10)\n');
DATAS = input(' ');
fprintf(1,'1. X축방향으로만 계산함(Y축단열)\n2. Y축방향으로만 계산함(X축단열)\n3. 양방향 모두 계산함\n');
DATAH = input(' ');
if DATAH == 1
DATAHY = 0;
elseif DATAH == 2
DATAHX = 0;
end;
N = DATAS*DATAS; % 슬라이스 한 수 제곱만큼 행렬 항 개수.
A = zeros(N,N+1); % 데이터들이 들어가고 마지막 열에는 상수
XL = zeros(1,N); % 풀이할때 써먹는 임시적인거
X1 = zeros(1,N); % 값 들어갈거
%이제 식을 세우는거야!
% 1 | 2 | 3
%---------------------
% 4 | 5 | 6
%---------------------
% 7 | 8 | 9
DELTAX=DATAX/DATAS;
DELTAY=DATAY/DATAS;
TWX= (2*DATAK*DELTAY/DELTAX)*((DATAHX/DATAK*DELTAX/2)/((DATAHX/DATAK*DELTAX/2)+1)); %X방향 벽 계수 (좌,우)
TWY= (2*DATAK*DELTAX/DELTAY)*((DATAHY/DATAK*DELTAY/2)/((DATAHY/DATAK*DELTAY/2)+1)); %Y방향 벽 계수 (상,하)
NORX= DATAK*DELTAY/DELTAX; %벽이 아닐때 좌우 계수
NORY= DATAK*DELTAX/DELTAY; %벽이 아닐때 상하 계수
NORQ= DATAQ*DELTAX*DELTAY; %내부에서 발생하는 열
for I = 1 : N
%1번 면적일때
if I <= DATAS && mod(I,DATAS) == 1
A(I,I)=-(TWX+TWY+NORX+NORY); %좌,상,우,하
A(I,I+1)=NORX; %오른쪽거
A(I,I+DATAS)=NORY; %아래쪽거
A(I,N+1)=-(NORQ+TWX*DATAF+TWY*DATAF); %내부발생열과 좌측,상단의 대류열이 상수임.
%2번 면적일때
elseif I<= DATAS && mod (I,DATAS) >=2 && mod(I,DATAS) <= DATAS-1
A(I,I)=-(NORX+TWY+NORX+NORY);
A(I,I-1)=NORX;%왼쪽
A(I,I+1)=NORX;%오른쪽
A(I,I+DATAS)=NORY;%아래쪽
A(I,N+1)=-(NORQ+TWY*DATAF); %내부발생열과 상단대류열
%3번 면적
elseif I<= DATAS && mod(I,DATAS) ==0
A(I,I)=-(NORX+TWY+TWX+NORY);
A(I,I-1)=NORX;%왼쪽
A(I,I+DATAS)=NORY;%아래쪽
A(I,N+1)=-(NORQ+TWX*DATAF+TWY*DATAF);
%4번 면적
elseif I <= DATAS*(DATAS-1) && I > DATAS && mod(I,DATAS) == 1
A(I,I)=-(TWX+NORY+NORX+NORY);
A(I,I+1)=NORX;%오른쪽
A(I,I-DATAS)=NORY;%위쪽
A(I,I+DATAS)=NORY;%아래쪽
A(I,N+1)=-(NORQ+TWX*DATAF);
%5번 면적
elseif I <= DATAS*(DATAS-1) && I > DATAS && mod (I,DATAS) >=2 && mod(I,DATAS) <= DATAS-1
A(I,I)=-(NORX+NORY+NORX+NORY);
A(I,I-1)=NORX;%왼쪽
A(I,I+1)=NORX;%오른쪽
A(I,I-DATAS)=NORY;%위쪽
A(I,I+DATAS)=NORY;%아래쪽
A(I,N+1)=-(NORQ);
%6번 면적
elseif I <= DATAS*(DATAS-1) && I > DATAS && mod(I,DATAS) ==0
A(I,I)=-(NORX+NORY+TWX+NORY);
A(I,I-1)=NORX;%왼쪽
A(I,I-DATAS)=NORY;%위쪽
A(I,I+DATAS)=NORY;%아래쪽
A(I,N+1)=-(NORQ+TWX*DATAF);
%7번 면적
elseif I > DATAS*(DATAS-1) && mod(I,DATAS) == 1
A(I,I)=-(TWX+NORY+NORX+TWY);
A(I,I+1)=NORX;%오른쪽
A(I,I-DATAS)=NORY;%위쪽
A(I,N+1)=-(NORQ+TWX*DATAF+TWY*DATAF);
%8번 면적
elseif I > DATAS*(DATAS-1) && mod (I,DATAS) >=2 && mod(I,DATAS) <= DATAS-1
A(I,I)=-(NORX+NORY+NORX+TWY);
A(I,I-1)=NORX;%왼쪽
A(I,I+1)=NORX;%오른쪽
A(I,I-DATAS)=NORY;%위쪽
A(I,N+1)=-(NORQ+TWY*DATAF);
%9번 면적
else
A(I,I)=-(NORX+NORY+TWX+TWY);
A(I,I-1)=NORX;%왼쪽
A(I,I+1)=NORX;%오른쪽
A(I,I-DATAS)=NORY;%위쪽
A(I,N+1)=-(NORQ+TWX*DATAF+TWY*DATAF);
end;
end;
%식세우기 끝!
OK=TRUE;
for I = 1 : N
for J = 1 : N+1
fprintf(1, '%3.4f\t', A(I,J));
end;
fprintf ('\n');
end;
if OK == TRUE
for I = 1 : N
XL(I) = 1;
end;
% STEP 1
if abs(A(1,1)) <= 1.0e-20
OK = FALSE;
else
A(1,1) = A(1,1)/XL(1);
% STEP 2
for J = 2 : N
A(1,J) = A(1,J)/XL(1);
A(J,1) = A(J,1)/A(1,1);
end;
% STEP 3
M = N-1;
I = 2;
while I <= M && OK == TRUE
% STEP 4
KK = I-1;
S = 0;
for K = 1 : KK
S = S-A(I,K)*A(K,I);
end;
A(I,I) = (A(I,I)+S)/XL(I);
if abs(A(I,I)) <= 1.0e-20
OK = FALSE;
else
% STEP 5
JJ = I+1;
for J = JJ : N
SS = 0;
S = 0;
for K = 1 : KK
SS = SS-A(I,K)*A(K,J);
S = S-A(J,K)*A(K,I);
end;
A(I,J) = (A(I,J)+SS)/XL(I);
A(J,I) = (A(J,I)+S)/A(I,I);
end;
end;
I = I+1;
end;
if OK == TRUE
S = 0;
for K = 1 : M
S = S-A(N,K)*A(K,N);
end;
A(N,N) = (A(N,N)+S)/XL(N);
end;
end;
if OK == FALSE
fprintf(1,'주어진 행렬이 LU분해가 가능하지 않습니다.\n');
end;
end;
% 여기부터 후진대입
I=2;
while I <= N
SS=0;
for S = 1 : I-1;
SS=SS-A(S,N+1)*A(I,S);
end;
A(I,N+1)=A(I,N+1)+SS;
I=I+1;
end;
% 여기부터 전진대입
X1(1,N)=A(N,N+1)/A(N,N);
I=1;
while I <= N-1
SS=0;
for S = 1 : I;
SS=SS+A(N-I,N-S+1)*X1(1,N-S+1);
end;
X1(N-I)=(A(N-I,N+1)-SS)/A(N-I,N-I);
I=I+1;
end;
%식을 출력하자! 끝났어!!
OUP = 1;
fprintf(OUP, 'LU-Decomposition 해법.\n\n');
fprintf(OUP, '솔루션은 다음과 같다. :\n');
CHCH=1;
if DATAH == 1
for I = 1 : N/DATAS
fprintf(OUP, '%11.8f\t', X1(I));
end;
elseif DATAH == 2
for I = 1 : N/DATAS
K=I*10;
fprintf(OUP, '%11.8f\t', X1(K));
end;
else
for I = 1 : N
fprintf(OUP, '%11.2f\t', X1(I));
if mod(CHCH,10) == 0
fprintf(OUP, '\n');
end;
CHCH=CHCH+1;
end;
Gauss Seidel 해법
Gauss-Seidel Method 프로그램 코드
syms('AA', 'OK', 'NAME', 'INP', 'N', 'I', 'J', 'A', 'X1');
syms('TOL', 'NN', 'K', 'ERR', 'S', 'FLAG', 'OUP','ERRK','CHCH','XTRUE','TRUEERR','NEWERR');
TRUE = 1;
FALSE = 0;
fprintf(1,'Gauss Seidel해법.\n');
fprintf(1,'주어진 h를 입력하시오.\n');
DATAHX = input(' ');
DATAHY = DATAHX;
fprintf(1,'주어진 k를 입력하시오.\n');
DATAK = input(' ');
fprintf(1,'주어진 q를 입력하시오.\n');
DATAQ = input(' ');
fprintf(1,'주어진 유체의 온도 T_fluid를 입력하시오.\n');
DATAF = input(' ');
fprintf(1,'바의 x방향 길이를 입력하시오.(0.005)\n');
DATAX = input(' ');
fprintf(1,'바의 y방향 길이를 입력하시오.(0.02)\n');
DATAY = input(' ');
fprintf(1,'미소면적의 개수 (n X n)일 때, n를 입력하시오.(10)\n');
DATAS = input(' ');
N = DATAS*DATAS; % 슬라이스 한 수 제곱만큼 행렬 항 개수.
A = zeros(N,N+1); % 데이터들이 들어가고 마지막 열에는 상수
XL = zeros(1,N); % 풀이할때 써먹는 임시적인거
X1 = zeros(1,N); % 값 들어갈거
XTRUE = zeros(N);
%이제 식을 세우는거야!
% 1 | 2 | 3
%---------------------
% 4 | 5 | 6
%---------------------
% 7 | 8 | 9
DELTAX=DATAX/DATAS;
DELTAY=DATAY/DATAS;
TWX= (2*DATAK*DELTAY/DELTAX)*((DATAHX/DATAK*DELTAX/2)/((DATAHX/DATAK*DELTAX/2)+1)); %X방향 벽 계수 (좌,우)
TWY= (2*DATAK*DELTAX/DELTAY)*((DATAHY/DATAK*DELTAY/2)/((DATAHY/DATAK*DELTAY/2)+1)); %Y방향 벽 계수 (상,하)
NORX= DATAK*DELTAY/DELTAX; %벽이 아닐때 좌우 계수
NORY= DATAK*DELTAX/DELTAY; %벽이 아닐때 상하 계수
NORQ= DATAQ*DELTAX*DELTAY; %내부에서 발생하는 열
for I = 1 : N
%1번 면적일때
if I <= DATAS && mod(I,DATAS) == 1
A(I,I)=-(TWX+TWY+NORX+NORY); %좌,상,우,하
A(I,I+1)=NORX; %오른쪽거
A(I,I+DATAS)=NORY; %아래쪽거
A(I,N+1)=-(NORQ+TWX*DATAF+TWY*DATAF); %내부발생열과 좌측,상단의 대류열이 상수임.
%2번 면적일때
elseif I<= DATAS && mod (I,DATAS) >=2 && mod(I,DATAS) <= DATAS-1
A(I,I)=-(NORX+TWY+NORX+NORY);
A(I,I-1)=NORX;%왼쪽
A(I,I+1)=NORX;%오른쪽
A(I,I+DATAS)=NORY;%아래쪽
A(I,N+1)=-(NORQ+TWY*DATAF); %내부발생열과 상단대류열
%3번 면적
elseif I<= DATAS && mod(I,DATAS) ==0
A(I,I)=-(NORX+TWY+TWX+NORY);
A(I,I-1)=NORX;%왼쪽
A(I,I+DATAS)=NORY;%아래쪽
A(I,N+1)=-(NORQ+TWX*DATAF+TWY*DATAF);
%4번 면적
elseif I <= DATAS*(DATAS-1) && I > DATAS && mod(I,DATAS) == 1
A(I,I)=-(TWX+NORY+NORX+NORY);
A(I,I+1)=NORX;%오른쪽
A(I,I-DATAS)=NORY;%위쪽
A(I,I+DATAS)=NORY;%아래쪽
A(I,N+1)=-(NORQ+TWX*DATAF);
%5번 면적
elseif I <= DATAS*(DATAS-1) && I > DATAS && mod (I,DATAS) >=2 && mod(I,DATAS) <= DATAS-1
A(I,I)=-(NORX+NORY+NORX+NORY);
A(I,I-1)=NORX;%왼쪽
A(I,I+1)=NORX;%오른쪽
A(I,I-DATAS)=NORY;%위쪽
A(I,I+DATAS)=NORY;%아래쪽
A(I,N+1)=-(NORQ);
%6번 면적
elseif I <= DATAS*(DATAS-1) && I > DATAS && mod(I,DATAS) ==0
A(I,I)=-(NORX+NORY+TWX+NORY);
A(I,I-1)=NORX;%왼쪽
A(I,I-DATAS)=NORY;%위쪽
A(I,I+DATAS)=NORY;%아래쪽
A(I,N+1)=-(NORQ+TWX*DATAF);
%7번 면적
elseif I > DATAS*(DATAS-1) && mod(I,DATAS) == 1
A(I,I)=-(TWX+NORY+NORX+TWY);
A(I,I+1)=NORX;%오른쪽
A(I,I-DATAS)=NORY;%위쪽
A(I,N+1)=-(NORQ+TWX*DATAF+TWY*DATAF);
%8번 면적
elseif I > DATAS*(DATAS-1) && mod (I,DATAS) >=2 && mod(I,DATAS) <= DATAS-1
A(I,I)=-(NORX+NORY+NORX+TWY);
A(I,I-1)=NORX;%왼쪽
A(I,I+1)=NORX;%오른쪽
A(I,I-DATAS)=NORY;%위쪽
A(I,N+1)=-(NORQ+TWY*DATAF);
%9번 면적
else
A(I,I)=-(NORX+NORY+TWX+TWY);
A(I,I-1)=NORX;%왼쪽
A(I,I+1)=NORX;%오른쪽
A(I,I-DATAS)=NORY;%위쪽
A(I,N+1)=-(NORQ+TWX*DATAF+TWY*DATAF);
end;
end;
%식세우기 끝!
fprintf(1,'Part1 에서 구한 답을 파일로 입력받는다. (파일명 : xtrue.txt)\n');
fprintf(1,'Drive:\\Filename.txt 형식으로 입력.\n');
fprintf(1,'예를들어 : A:\\xtrue.txt\n');
NAME = input(' ','s');
INP = fopen(NAME,'rt');
for I = 1 : N
XTRUE(I) = fscanf(INP, '%f',1);
end;
for I = 1 : N
X1(I) = 0;
end;
OK = FALSE;
while OK == FALSE
fprintf(1,'유효숫자와 반복회수 둘중 하나의 값은 꼭 입력해야합니다.\n');
fprintf(1,'원하는 유효숫자의 개수를 양의 정수로 입력하세요(사용하지 않을려면 0)\n');
SIGNIFI = input(' ');
TOL = 0.5*10^(-SIGNIFI);
if TOL > 0
OK = TRUE;
else
fprintf(1,'양수를 입력하세요\n');
end;
end;
OK = FALSE;
while OK == FALSE
fprintf(1,'최대 반복회수를 입력하세요 (유효숫자개수에 따른 반복회수를 구할때는 크게 입력하세요.)\n');
NN = input(' ');
if NN > 0
OK = TRUE;
ERRK = zeros(1,NN);
else
fprintf(1,'양의 정수로 입력하세요\n');
end;
end;
if OK == TRUE
% STEP 1
K = 1;
OK = FALSE;
% STEP 2
while OK == FALSE && K <= NN
% ERR is used to test accuracy - it measures the infinity-norm
ERR = 0;
TRUEERR=0;
NEWERR=0;
% STEP 3
for I = 1 : N
S = 0;
for J = 1 : N
S = S-A(I,J)*X1(J);
end;
S = (S+A(I,N+1))/A(I,I);
if abs(S) > ERR
ERR = abs(S);
end;
X1(I) = X1(I) + S;
TRUEERR=(XTRUE(I)-X1(I))/XTRUE(I);
if TRUEERR >= NEWERR
NEWERR=TRUEERR;
end;
end;
% STEP 4
if SIGNIFI >= 1
if ERR <= TOL
OK = TRUE;
% process is complete
else
% STEP 6 - is not used since only one vector is required
end;
end;
if K == NN
OK = TRUE;
end;
K = K+1;
ERRK(K)=NEWERR;
end;
if OK == FALSE
fprintf(1,'최대 반복 횟수를 넘어섰습니다.\n');
% STEP 7
% procedure completed unsuccessfully
else
OUP = 1;
fprintf(OUP, 'Gauss - Seidel 반복 해법.\n\n');
fprintf(OUP, '솔루션은 다음과 같다. :\n');
CHCH=1;
for I = 1 : N
fprintf(OUP, '%11.8f\t', X1(I));
if mod(CHCH,10) == 0
fprintf(OUP, '\n');
end;
CHCH=CHCH+1;
end;
if SIGNIFI >= 1
fprintf(1, '유효숫자 %d 개를 구하는 시행회수는 %d이다.\n', SIGNIFI,K);
else
K=K-1;
fprintf(1, '\n%d 회의 반복을 시행했다.\n', K);
fprintf(OUP, '각 반복의 오차는 다음과 같다.\n');
CHCH=1;
for I = 1 : K
fprintf(OUP, '%d회 반복 : \t%11.10f\n',I, ERRK(I));
CHCH=CHCH+1;
end;
end;
fprintf(OUP, '\n\n');
end;
end;
Jacobi 해법
Jacobi Method 프로그램 코드
syms('AA', 'OK', 'NAME', 'INP', 'N', 'I', 'J', 'A', 'X1');
syms('TOL', 'NN', 'K', 'ERR', 'S', 'X2', 'FLAG', 'OUP','ERRK','CHCH','XTRUE','TRUEERR','NEWERR');
TRUE = 1;
FALSE = 0;
fprintf(1,'Jacobi 해법.\n');
fprintf(1,'주어진 h를 입력하시오.\n');
DATAHX = input(' ');
DATAHY = DATAHX;
fprintf(1,'주어진 k를 입력하시오.\n');
DATAK = input(' ');
fprintf(1,'주어진 q를 입력하시오.\n');
DATAQ = input(' ');
fprintf(1,'주어진 유체의 온도 T_fluid를 입력하시오.\n');
DATAF = input(' ');
fprintf(1,'바의 x방향 길이를 입력하시오.(0.005)\n');
DATAX = input(' ');
fprintf(1,'바의 y방향 길이를 입력하시오.(0.02)\n');
DATAY = input(' ');
fprintf(1,'미소면적의 개수 (n X n)일 때, n를 입력하시오.(10)\n');
DATAS = input(' ');
N = DATAS*DATAS; % 슬라이스 한 수 제곱만큼 행렬 항 개수.
A = zeros(N,N+1); % 데이터들이 들어가고 마지막 열에는 상수
XL = zeros(1,N); % 풀이할때 써먹는 임시적인거
X1 = zeros(1,N); % 값 들어갈거
X2 = zeros(N);
XTRUE = zeros(N);
%이제 식을 세우는거야!
% 1 | 2 | 3
%---------------------
% 4 | 5 | 6
%---------------------
% 7 | 8 | 9
DELTAX=DATAX/DATAS;
DELTAY=DATAY/DATAS;
TWX= (2*DATAK*DELTAY/DELTAX)*((DATAHX/DATAK*DELTAX/2)/((DATAHX/DATAK*DELTAX/2)+1)); %X방향 벽 계수 (좌,우)
TWY= (2*DATAK*DELTAX/DELTAY)*((DATAHY/DATAK*DELTAY/2)/((DATAHY/DATAK*DELTAY/2)+1)); %Y방향 벽 계수 (상,하)
NORX= DATAK*DELTAY/DELTAX; %벽이 아닐때 좌우 계수
NORY= DATAK*DELTAX/DELTAY; %벽이 아닐때 상하 계수
NORQ= DATAQ*DELTAX*DELTAY; %내부에서 발생하는 열
for I = 1 : N
%1번 면적일때
if I <= DATAS && mod(I,DATAS) == 1
A(I,I)=-(TWX+TWY+NORX+NORY); %좌,상,우,하
A(I,I+1)=NORX; %오른쪽거
A(I,I+DATAS)=NORY; %아래쪽거
A(I,N+1)=-(NORQ+TWX*DATAF+TWY*DATAF); %내부발생열과 좌측,상단의 대류열이 상수임.
%2번 면적일때
elseif I<= DATAS && mod (I,DATAS) >=2 && mod(I,DATAS) <= DATAS-1
A(I,I)=-(NORX+TWY+NORX+NORY);
A(I,I-1)=NORX;%왼쪽
A(I,I+1)=NORX;%오른쪽
A(I,I+DATAS)=NORY;%아래쪽
A(I,N+1)=-(NORQ+TWY*DATAF); %내부발생열과 상단대류열
%3번 면적
elseif I<= DATAS && mod(I,DATAS) ==0
A(I,I)=-(NORX+TWY+TWX+NORY);
A(I,I-1)=NORX;%왼쪽
A(I,I+DATAS)=NORY;%아래쪽
A(I,N+1)=-(NORQ+TWX*DATAF+TWY*DATAF);
%4번 면적
elseif I <= DATAS*(DATAS-1) && I > DATAS && mod(I,DATAS) == 1
A(I,I)=-(TWX+NORY+NORX+NORY);
A(I,I+1)=NORX;%오른쪽
A(I,I-DATAS)=NORY;%위쪽
A(I,I+DATAS)=NORY;%아래쪽
A(I,N+1)=-(NORQ+TWX*DATAF);
%5번 면적
elseif I <= DATAS*(DATAS-1) && I > DATAS && mod (I,DATAS) >=2 && mod(I,DATAS) <= DATAS-1
A(I,I)=-(NORX+NORY+NORX+NORY);
A(I,I-1)=NORX;%왼쪽
A(I,I+1)=NORX;%오른쪽
A(I,I-DATAS)=NORY;%위쪽
A(I,I+DATAS)=NORY;%아래쪽
A(I,N+1)=-(NORQ);
%6번 면적
elseif I <= DATAS*(DATAS-1) && I > DATAS && mod(I,DATAS) ==0
A(I,I)=-(NORX+NORY+TWX+NORY);
A(I,I-1)=NORX;%왼쪽
A(I,I-DATAS)=NORY;%위쪽
A(I,I+DATAS)=NORY;%아래쪽
A(I,N+1)=-(NORQ+TWX*DATAF);
%7번 면적
elseif I > DATAS*(DATAS-1) && mod(I,DATAS) == 1
A(I,I)=-(TWX+NORY+NORX+TWY);
A(I,I+1)=NORX;%오른쪽
A(I,I-DATAS)=NORY;%위쪽
A(I,N+1)=-(NORQ+TWX*DATAF+TWY*DATAF);
%8번 면적
elseif I > DATAS*(DATAS-1) && mod (I,DATAS) >=2 && mod(I,DATAS) <= DATAS-1
A(I,I)=-(NORX+NORY+NORX+TWY);
A(I,I-1)=NORX;%왼쪽
A(I,I+1)=NORX;%오른쪽
A(I,I-DATAS)=NORY;%위쪽
A(I,N+1)=-(NORQ+TWY*DATAF);
%9번 면적
else
A(I,I)=-(NORX+NORY+TWX+TWY);
A(I,I-1)=NORX;%왼쪽
A(I,I+1)=NORX;%오른쪽
A(I,I-DATAS)=NORY;%위쪽
A(I,N+1)=-(NORQ+TWX*DATAF+TWY*DATAF);
end;
end;
%식세우기 끝!
fprintf(1,'Part1 에서 구한 답을 파일로 입력받는다. (파일명 : xtrue.txt)\n');
fprintf(1,'Drive:\\Filename.txt 형식으로 입력.\n');
fprintf(1,'예를들어 : A:\\xtrue.txt\n');
NAME = input(' ','s');
INP = fopen(NAME,'rt');
for I = 1 : N
XTRUE(I) = fscanf(INP, '%f',1);
end;
for I = 1 : N
X1(I) = 0;
end;
OK = FALSE;
while OK == FALSE
fprintf(1,'유효숫자와 반복회수 둘중 하나의 값은 꼭 입력해야합니다.\n');
fprintf(1,'원하는 유효숫자의 개수를 양의 정수로 입력하세요(사용하지 않을려면 0)\n');
SIGNIFI = input(' ');
TOL = 0.5*10^(-SIGNIFI);
if TOL > 0
OK = TRUE;
else
fprintf(1,'양수를 입력하세요\n');
end;
end;
OK = FALSE;
while OK == FALSE
fprintf(1,'최대 반복회수를 입력하세요 (유효숫자개수에 따른 반복회수를 구할때는 크게 입력하세요.)\n');
NN = input(' ');
if NN > 0
OK = TRUE;
ERRK = zeros(1,NN);
else
fprintf(1,'양의 정수로 입력하세요\n');
end;
end;
if OK == TRUE
% STEP 1
K = 1;
OK = FALSE;
% STEP 2
while OK == FALSE && K <= NN
% err is used to test accuracy - it measures the infinity-norm
ERR = 0;
TRUEERR=0;
NEWERR=0;
% STEP 3
for I = 1 : N
S = 0;
for J = 1 : N
S = S-A(I,J)*X1(J);
end;
S = (S+A(I,N+1))/A(I,I);
if abs(S) > ERR
ERR = abs(S);
end;
% use X2 for X
X2(I) = X1(I)+S;
TRUEERR=(XTRUE(I)-X2(I))/XTRUE(I);
if TRUEERR >= NEWERR
NEWERR=TRUEERR;
end;
end;
% STEP 4
if ERR <= TOL && SIGNIFI >=1
% process is complete
OK = TRUE;
else
% STEP 5
if K==NN
OK =TRUE;
end;
K = K+1;
ERRK(K)=NEWERR;
% STEP 6
for I = 1 : N
X1(I) = X2(I);
end;
end;
end;
if OK == FALSE
fprintf(1,'최대 반복 횟수를 넘어섰습니다.\n');
% STEP 7
% procedure completed unsuccessfully
else
OUP = 1;
fprintf(OUP, 'JACOBI 반복해법 \n\n');
fprintf(OUP, '솔루션은 다음과 같다. :\n');
CHCH=1;
for I = 1 : N
fprintf(OUP, '%11.8f\t', X1(I));
if mod(CHCH,10) == 0
fprintf(OUP, '\n');
end;
CHCH=CHCH+1;
end;
if SIGNIFI >= 1
fprintf(1, '유효숫자 %d 개를 구하는 시행회수는 %d이다.\n', SIGNIFI,K);
else
K=K-1;
fprintf(1, '\n%d 회의 반복을 시행했다.\n', K);
fprintf(OUP, '각 반복의 오차는 다음과 같다.\n');
CHCH=1;
for I = 1 : K
fprintf(OUP, '%d회 반복 : \t%11.10f\n',I, ERRK(I));
CHCH=CHCH+1;
end;
end;
fprintf(OUP, '\n\n');
end;
end;
RECENT COMMENT