'수치해석' 태그의 글 목록

신변잡기/공대딩 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;

%식세우기 끝!

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;

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;

if OK == FALSE

fprintf(1,'주어진 행렬이 LU분해가 가능하지 않습니다.\n');

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;

%식세우기 끝!

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;

OK = FALSE;

while OK == FALSE

fprintf(1,'최대 반복회수를 입력하세요 (유효숫자개수에 따른 반복회수를 구할때는 크게 입력하세요.)\n');

NN = input(' ');

if NN > 0

OK = TRUE;

ERRK = zeros(1,NN);

else

fprintf(1,'양의 정수로 입력하세요\n');

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;

% 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;

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;

fprintf(OUP, '\n\n');

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;

%식세우기 끝!

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;

OK = FALSE;

while OK == FALSE

fprintf(1,'최대 반복회수를 입력하세요 (유효숫자개수에 따른 반복회수를 구할때는 크게 입력하세요.)\n');

NN = input(' ');

if NN > 0

OK = TRUE;

ERRK = zeros(1,NN);

else

fprintf(1,'양의 정수로 입력하세요\n');

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;

% 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;

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;

fprintf(OUP, '\n\n');

end;

,

신변잡기/공대딩 2009. 4. 2. 20:03

수치해석 MATLAB으로 exponential 값을 근사하라.

clear all

close all

b = 2; %구간 설정

m = 40; %구간개수 설정

degree = 5; %전개할 테일러함수의 항 개수

degrees = 5; %마지막으로 출력할 테일러함수의 항 개수 (degree보다 크다)

a = 0; %f(x-a) 의 a 값

delta_x = 2*b/m; %구간 사이의 거리

func = exp(a); %테일러함수 0째 항

eps_t = 0; %오차값 초기화

f_true = 0; %실제값 초기화

j = 1; %테일러함수 1항으로부터 시작.

while degree <=degrees

fprintf ('\ndegree : %d\n',degree);

fprintf ('----x값\t테일러함수값\t--실제함수값\t\t오차\n');

x=-b; %x 는 -b로부터

while x<b+delta_x %b 까지 계산한다.

f_true = exp(x); %실제값은 e^x

y=funct(a,j,x,func,degree); %각 x값마다 계산한다.

eps_t=(f_true-y)/f_true; %오차값 계산

fprintf('%2.4f\t%2.7f\t%2.7f\t%2.7f\n', x,y,f_true,eps_t);

%x값, 테일러함수, 실제값, 오차 출력

x = x+delta_x; %x를 증분 delta_x만큼 증가.

end

degree=degree+1; %다음항 계산분까지 출력 증가.

end

funct.m 파일 (funct 함수를 이용하기 위해.)

function y=funct(a,j,x,func,degree)

y=func;

while j<=degree

y=y+exp(a)/factorial(j).*(x-a).^j;

j=j+1;

end

결과

degree:1

x값	테일러값	실제함수값	오차
-2.0000	-1.0000000	0.1353353	8.3890561
-1.9000	-0.9000000	0.1495686	7.0173050
-1.8000	-0.8000000	0.1652989	5.8397180
-1.7000	-0.7000000	0.1826835	4.8317632
-1.6000	-0.6000000	0.2018965	3.9718195
-1.5000	-0.5000000	0.2231302	3.2408445
-1.4000	-0.4000000	0.2465970	2.6220800
-1.3000	-0.3000000	0.2725318	2.1007890
-1.2000	-0.2000000	0.3011942	1.6640234
-1.1000	-0.1000000	0.3328711	1.3004166
-1.0000	0.0000000	0.3678794	1.0000000
-0.9000	0.1000000	0.4065697	0.7540397
-0.8000	0.2000000	0.4493290	0.5548918
-0.7000	0.3000000	0.4965853	0.3958742
-0.6000	0.4000000	0.5488116	0.2711525
-0.5000	0.5000000	0.6065307	0.1756394
-0.4000	0.6000000	0.6703200	0.1049052
-0.3000	0.7000000	0.7408182	0.0550988
-0.2000	0.8000000	0.8187308	0.0228778
-0.1000	0.9000000	0.9048374	0.0053462
0.0000	1.0000000	1.0000000	0.0000000
0.1000	1.1000000	1.1051709	0.0046788
0.2000	1.2000000	1.2214028	0.0175231
0.3000	1.3000000	1.3498588	0.0369363
0.4000	1.4000000	1.4918247	0.0615519
0.5000	1.5000000	1.6487213	0.0902040
0.6000	1.6000000	1.8221188	0.1219014
0.7000	1.7000000	2.0137527	0.1558050
0.8000	1.8000000	2.2255409	0.1912079
0.9000	1.9000000	2.4596031	0.2275176
1.0000	2.0000000	2.7182818	0.2642411
1.1000	2.1000000	3.0041660	0.3009707
1.2000	2.2000000	3.3201169	0.3373727
1.3000	2.3000000	3.6692967	0.3731769
1.4000	2.4000000	4.0552000	0.4081673
1.5000	2.5000000	4.4816891	0.4421746
1.6000	2.6000000	4.9530324	0.4750691
1.7000	2.7000000	5.4739474	0.5067545
1.8000	2.8000000	6.0496475	0.5371631
1.9000	2.9000000	6.6858944	0.5662510
2.0000	3.0000000	7.3890561	0.5939942

degree:2

----x값	테일러함수값	--실제함수값	오차
-2.0000	1.0000000	0.1353353	-6.3890561
-1.9000	0.9050000	0.1495686	-5.0507345
-1.8000	0.8200000	0.1652989	-3.9607109
-1.7000	0.7450000	0.1826835	-3.0780908
-1.6000	0.6800000	0.2018965	-2.3680620
-1.5000	0.6250000	0.2231302	-1.8010557
-1.4000	0.5800000	0.2465970	-1.3520160
-1.3000	0.5450000	0.2725318	-0.9997667
-1.2000	0.5200000	0.3011942	-0.7264608
-1.1000	0.5050000	0.3328711	-0.5171038
-1.0000	0.5000000	0.3678794	-0.3591409
-0.9000	0.5050000	0.4065697	-0.2420996
-0.8000	0.5200000	0.4493290	-0.1572813
-0.7000	0.5450000	0.4965853	-0.0974952
-0.6000	0.5800000	0.5488116	-0.0568289
-0.5000	0.6250000	0.6065307	-0.0304508
-0.4000	0.6800000	0.6703200	-0.0144408
-0.3000	0.7450000	0.7408182	-0.0056448
-0.2000	0.8200000	0.8187308	-0.0015503
-0.1000	0.9050000	0.9048374	-0.0001797
0.0000	1.0000000	1.0000000	0.0000000
0.1000	1.1050000	1.1051709	0.0001547
0.2000	1.2200000	1.2214028	0.0011485
0.3000	1.3450000	1.3498588	0.0035995
0.4000	1.4800000	1.4918247	0.0079263
0.5000	1.6250000	1.6487213	0.0143877
0.6000	1.7800000	1.8221188	0.0231153
0.7000	1.9450000	2.0137527	0.0341416
0.8000	2.1200000	2.2255409	0.0474226
0.9000	2.3050000	2.4596031	0.0628569
1.0000	2.5000000	2.7182818	0.0803014
1.1000	2.7050000	3.0041660	0.0995837
1.2000	2.9200000	3.3201169	0.1205129
1.3000	3.1450000	3.6692967	0.1428875
1.4000	3.3800000	4.0552000	0.1665023
1.5000	3.6250000	4.4816891	0.1911532
1.6000	3.8800000	4.9530324	0.2166415
1.7000	4.1450000	5.4739474	0.2427768
1.8000	4.4200000	6.0496475	0.2693789
1.9000	4.7050000	6.6858944	0.2962796
2.0000	5.0000000	7.3890561	0.3233236

degree : 3

x값	테일러함수값	실제함수값	오차
-2.0000	-0.3333333	0.1353353	3.4630187
-1.9000	-0.2381667	0.1495686	2.5923572
-1.8000	-0.1520000	0.1652989	1.9195464
-1.7000	-0.0738333	0.1826835	1.4041598
-1.6000	-0.0026667	0.2018965	1.0132081
-1.5000	0.0625000	0.2231302	0.7198944
-1.4000	0.1226667	0.2465970	0.5025621
-1.3000	0.1788333	0.2725318	0.3438074
-1.2000	0.2320000	0.3011942	0.2297329
-1.1000	0.2831667	0.3328711	0.1493203
-1.0000	0.3333333	0.3678794	0.0939061
-0.9000	0.3835000	0.4065697	0.0567422
-0.8000	0.4346667	0.4493290	0.0326315
-0.7000	0.4878333	0.4965853	0.0176243
-0.6000	0.5440000	0.5488116	0.0087674
-0.5000	0.6041667	0.6065307	0.0038976
-0.4000	0.6693333	0.6703200	0.0014720
-0.3000	0.7405000	0.7408182	0.0004296
-0.2000	0.8186667	0.8187308	0.0000783
-0.1000	0.9048333	0.9048374	0.0000045
0.0000	1.0000000	1.0000000	0.0000000
0.1000	1.1051667	1.1051709	0.0000038
0.2000	1.2213333	1.2214028	0.0000568
0.3000	1.3495000	1.3498588	0.0002658
0.4000	1.4906667	1.4918247	0.0007763
0.5000	1.6458333	1.6487213	0.0017516
0.6000	1.8160000	1.8221188	0.0033581
0.7000	2.0021667	2.0137527	0.0057535
0.8000	2.2053333	2.2255409	0.0090799
0.9000	2.4265000	2.4596031	0.0134587
1.0000	2.6666667	2.7182818	0.0189882
1.1000	2.9268333	3.0041660	0.0257418
1.2000	3.2080000	3.3201169	0.0337690
1.3000	3.5111667	3.6692967	0.0430955
1.4000	3.8373333	4.0552000	0.0537253
1.5000	4.1875000	4.4816891	0.0656425
1.6000	4.5626667	4.9530324	0.0788135
1.7000	4.9638333	5.4739474	0.0931894
1.8000	5.3920000	6.0496475	0.1087084
1.9000	5.8481667	6.6858944	0.1252978
2.0000	6.3333333	7.3890561	0.1428765

degree : 4

x값	테일러함수값	실제함수값	오차
-2.0000	0.3333333	0.1353353	-1.4630187
-1.9000	0.3048375	0.1495686	-1.0381113
-1.8000	0.2854000	0.1652989	-0.7265694
-1.7000	0.2741708	0.1826835	-0.5007967
-1.6000	0.2704000	0.2018965	-0.3393000
-1.5000	0.2734375	0.2231302	-0.2254619
-1.4000	0.2827333	0.2465970	-0.1465402
-1.3000	0.2978375	0.2725318	-0.0928541
-1.2000	0.3184000	0.3011942	-0.0571252
-1.1000	0.3441708	0.3328711	-0.0339463
-1.0000	0.3750000	0.3678794	-0.0193557
-0.9000	0.4108375	0.4065697	-0.0104972
-0.8000	0.4517333	0.4493290	-0.0053510
-0.7000	0.4978375	0.4965853	-0.0025216
-0.6000	0.5494000	0.5488116	-0.0010721
-0.5000	0.6067708	0.6065307	-0.0003960
-0.4000	0.6704000	0.6703200	-0.0001193
-0.3000	0.7408375	0.7408182	-0.0000260
-0.2000	0.8187333	0.8187308	-0.0000032
-0.1000	0.9048375	0.9048374	-0.0000001
0.0000	1.0000000	1.0000000	0.0000000
0.1000	1.1051708	1.1051709	0.0000001
0.2000	1.2214000	1.2214028	0.0000023
0.3000	1.3498375	1.3498588	0.0000158
0.4000	1.4917333	1.4918247	0.0000612
0.5000	1.6484375	1.6487213	0.0001721
0.6000	1.8214000	1.8221188	0.0003945
0.7000	2.0121708	2.0137527	0.0007855
0.8000	2.2224000	2.2255409	0.0014113
0.9000	2.4538375	2.4596031	0.0023441
1.0000	2.7083333	2.7182818	0.0036598
1.1000	2.9878375	3.0041660	0.0054353
1.2000	3.2944000	3.3201169	0.0077458
1.3000	3.6301708	3.6692967	0.0106630
1.4000	3.9974000	4.0552000	0.0142533
1.5000	4.3984375	4.4816891	0.0185759
1.6000	4.8357333	4.9530324	0.0236823
1.7000	5.3118375	5.4739474	0.0296148
1.8000	5.8294000	6.0496475	0.0364067
1.9000	6.3911708	6.6858944	0.0440814
2.0000	7.0000000	7.3890561	0.0526530

degree : 5

x값	테일러함수값	실제함수값	오차
-2.0000	0.0666667	0.1353353	0.5073963
-1.9000	0.0984959	0.1495686	0.3414667
-1.8000	0.1279360	0.1652989	0.2260323
-1.7000	0.1558494	0.1826835	0.1468885
-1.6000	0.1830187	0.2018965	0.0935026
-1.5000	0.2101563	0.2231302	0.0581450
-1.4000	0.2379147	0.2465970	0.0352085
-1.3000	0.2668964	0.2725318	0.0206779
-1.2000	0.2976640	0.3011942	0.0117207
-1.1000	0.3307499	0.3328711	0.0063723
-1.0000	0.3666667	0.3678794	0.0032967
-0.9000	0.4059168	0.4065697	0.0016059
-0.8000	0.4490027	0.4493290	0.0007262
-0.7000	0.4964369	0.4965853	0.0002988
-0.6000	0.5487520	0.5488116	0.0001087
-0.5000	0.6065104	0.6065307	0.0000334
-0.4000	0.6703147	0.6703200	0.0000080
-0.3000	0.7408173	0.7408182	0.0000013
-0.2000	0.8187307	0.8187308	0.0000001
-0.1000	0.9048374	0.9048374	0.0000000
0.0000	1.0000000	1.0000000	0.0000000
0.1000	1.1051709	1.1051709	0.0000000
0.2000	1.2214027	1.2214028	0.0000001
0.3000	1.3498578	1.3498588	0.0000008
0.4000	1.4918187	1.4918247	0.0000040
0.5000	1.6486979	1.6487213	0.0000142
0.6000	1.8220480	1.8221188	0.0000389
0.7000	2.0135714	2.0137527	0.0000900
0.8000	2.2251307	2.2255409	0.0001843
0.9000	2.4587583	2.4596031	0.0003435
1.0000	2.7166667	2.7182818	0.0005942
1.1000	3.0012584	3.0041660	0.0009679
1.2000	3.3151360	3.3201169	0.0015002
1.3000	3.6611119	3.6692967	0.0022306
1.4000	4.0422187	4.0552000	0.0032011
1.5000	4.4617188	4.4816891	0.0044560
1.6000	4.9231147	4.9530324	0.0060403
1.7000	5.4301589	5.4739474	0.0079994
1.8000	5.9868640	6.0496475	0.0103780
1.9000	6.5975124	6.6858944	0.0132192
2.0000	7.2666667	7.3890561	0.0165636

x	degree : 1	degree : 2	degree : 3	degree : 4	degree : 5	TRUE
-2	-1	1	-0.3333	0.3333	0.0667	0.1353
-1.9	-0.9	0.905	-0.2382	0.3048	0.0985	0.1496
-1.8	-0.8	0.82	-0.152	0.2854	0.1279	0.1653
-1.7	-0.7	0.745	-0.0738	0.2742	0.1558	0.1827
-1.6	-0.6	0.68	-0.0027	0.2704	0.183	0.2019
-1.5	-0.5	0.625	0.0625	0.2734	0.2102	0.2231
-1.4	-0.4	0.58	0.1227	0.2827	0.2379	0.2466
-1.3	-0.3	0.545	0.1788	0.2978	0.2669	0.2725
-1.2	-0.2	0.52	0.232	0.3184	0.2977	0.3012
-1.1	-0.1	0.505	0.2832	0.3442	0.3307	0.3329
-1	0	0.5	0.3333	0.375	0.3667	0.3679
-0.9	0.1	0.505	0.3835	0.4108	0.4059	0.4066
-0.8	0.2	0.52	0.4347	0.4517	0.449	0.4493
-0.7	0.3	0.545	0.4878	0.4978	0.4964	0.4966
-0.6	0.4	0.58	0.544	0.5494	0.5488	0.5488
-0.5	0.5	0.625	0.6042	0.6068	0.6065	0.6065
-0.4	0.6	0.68	0.6693	0.6704	0.6703	0.6703
-0.3	0.7	0.745	0.7405	0.7408	0.7408	0.7408
-0.2	0.8	0.82	0.8187	0.8187	0.8187	0.8187
-0.1	0.9	0.905	0.9048	0.9048	0.9048	0.9048
0	1	1	1	1	1	1
0.1	1.1	1.105	1.1052	1.1052	1.1052	1.1052
0.2	1.2	1.22	1.2213	1.2214	1.2214	1.2214
0.3	1.3	1.345	1.3495	1.3498	1.3499	1.3499
0.4	1.4	1.48	1.4907	1.4917	1.4918	1.4918
0.5	1.5	1.625	1.6458	1.6484	1.6487	1.6487
0.6	1.6	1.78	1.816	1.8214	1.822	1.8221
0.7	1.7	1.945	2.0022	2.0122	2.0136	2.0138
0.8	1.8	2.12	2.2053	2.2224	2.2251	2.2255
0.9	1.9	2.305	2.4265	2.4538	2.4588	2.4596
1	2	2.5	2.6667	2.7083	2.7167	2.7183
1.1	2.1	2.705	2.9268	2.9878	3.0013	3.0042
1.2	2.2	2.92	3.208	3.2944	3.3151	3.3201
1.3	2.3	3.145	3.5112	3.6302	3.6611	3.6693
1.4	2.4	3.38	3.8373	3.9974	4.0422	4.0552
1.5	2.5	3.625	4.1875	4.3984	4.4617	4.4817
1.6	2.6	3.88	4.5627	4.8357	4.9231	4.953
1.7	2.7	4.145	4.9638	5.3118	5.4302	5.4739
1.8	2.8	4.42	5.392	5.8294	5.9869	6.0496
1.9	2.9	4.705	5.8482	6.3912	6.5975	6.6859
2	3	5	6.3333	7	7.2667	7.3891

,

신변잡기/공대딩 2009. 4. 2. 19:53

수치해석 MATLAB으로 테일러 전개로 sin(x)의 값을 근사하라.

f=sin(pi/4); %구하고 싶은 함수. 참값.

x=1; %func값을 알고있을 때의 x값

j=0; %while 구문을 0부터 실행한다.

n=3; %유효숫자 개수

eps_s=0.5*10^(2-n);

known=1;

unknown=pi/4;

func=0; %테일러함수 func를 초기화한다.

func_before=0; %이전 테일러함수 func_before를 초기화한다.

fprintf ('terms\t근사오차\t　　상대오차\t　　　함수값\n')

while j<100 %초항 0부터 100항까지 반복

if mod(j,4)==0

plusfunc=sin(x); %0,4,8.. 항은 sinx

elseif mod(j,4)==1

plusfunc=cos(x); %1,5,9.. 항은 cosx

elseif mod(j,4)==2

plusfunc=-sin(x); %2,6,10.. 항은 -sinx

else

plusfunc=-cos(x); %3,7,11.. 항은 -cosx

end

func=func+plusfunc/factorial(j)*((unknown-known)^j); %테일러급수전개.

eps_a=(func-func_before)/func*100; %근사오차 = (현재값-이전값)/현재값

eps_t=(f-func)/f*100; %상대오차 = (참값-현재값)/참값

fprintf('%2.0f\t%2.7f\t%2.7f\t%2.7f\n', j,eps_a,eps_t,func)

%각 회수마다 상대오차, 근사오차, 함수값을 출력한다.

if abs(eps_a) < eps_s

break %근사오차가 백분율허용치보다 작아지면 멈춘다.

end

func_before=func; %이전값을 현재 값으로 교체한다.

j=j+1; %회수를 1회 증가한다.

end

fprintf('\n실제값 : %2.10f, ',f)

fprintf('근사값 : %2.10f',func)

결과

terms	근사오차	상대오차	함수값
0	100.0000000	-19.0019679	0.8414710
1	-15.9815978	-2.6041804	0.7255211
2	-2.7439891	0.1360748	0.7061446
3	0.1258767	0.0102110	0.7070346
4	0.0105166	-0.0003056	0.7071089

매트랩은 처음 다뤄봐서 코딩이 다소 지저분한 듯도 하다.
어쨌든 결과값은 제대로 출력 되는듯 :)

,

모놀로9

검색결과 리스트

수치해석에 해당되는 글 3건

글

수치해석 MATLAB 열전달 Gauss Seidel, Jacobi, LU Decomposition

트랙백

댓글

글

수치해석 MATLAB으로 exponential 값을 근사하라.

트랙백

댓글

글

수치해석 MATLAB으로 테일러 전개로 sin(x)의 값을 근사하라.

트랙백

댓글

사이드 메뉴

CATEGORY

RECENT POSTS

TAG

RECENT TRACKBACK

ARCHIVE

RECENT COMMENT

CALENDAR

MY LINK

COUNTER

티스토리툴바