원자로 수치해석 과제 #2

: 1-D, 1-G Neutron Diffusion Eigenvalue Problem

제출 일자 : 2010/03/26

2008-11833. 반영석

Part 1. Discretization and Direct Solution

1. Input processing

주어진 값들을 파일로부터 입력받아 global 변수에 저장하는 함수이다. .inp 파일에는 line 1부터 line 6 까지 6개의 블록으로 되어있는데 순서대로 다음 정보가 들어있다.

#1. material의 수를 나타내는 M과 그 크기를 나타내는 [eqn: xi]

#2. 4쌍의 D,[eqn: SIGMA_a], [eqn: nu SIGMA_f]

#3. #2.에서의 물질의 구성순서를 나타내는 숫자

#4. Boundary condition (albedo)

#5. 하나의 domain이 n개의 mesh로 나누어 져있음을 나타내는 n

#6. eigenvalue와 fission source의 convergence criteria

수업자료 게시판에서 조교님이 올려주신 코드를 참고하여 한줄 씩 읽어가면서 변수들을 global 형식으로 저장해두고 읽은 결과를 출력해주는 함수(readinp.m) 를 작성하였다. line 2에서는 composition 의 첫줄로부터 composition의 종류의 수를 입력받고 그 수만큼 한줄 씩 D,[eqn: SIGMA_a], [eqn: nu SIGMA_f]를 읽어나가는 방법을 택하였다. 이 수들은 임시로 저장되고 실제의 m번째 domain에서의 값은 line 3에서 type을 읽으면서 저장하기로 하였다. line 3에서는 구성을 읽기 위해서 10진수를 2진수로 바꾸는 것처럼 제일 높은 자릿수부터 한자리씩 빼 가면서 읽는 방법을 택하였다. 그 외의 line 들은 딱히 특징 없이 한 줄씩 FGETL을 해주면서 값을 저장해 나갔다. 아래는 그 코드이다.

function readinp

global M xi H D sig_a sig_f al_l al_r n ntot conv_eig conv_sol

% open input file

%--------------------# line 1 ---------------------

fid = fopen('nd1d_1.inp');

% read one line and display

% the first line is ignored

tline = fgetl(fid);

disp(tline);

% read one line and display

tline = fgetl(fid);

% read formatted data from string

tmparr = sscanf(tline, '%f');

M=tmparr(1);xi=tmparr(2);H=M*xi;

c=zeros(1,M);

str = sprintf('M = %f,  xi = %f, H = %f\n', M, xi, H);

disp(str);

%--------------------# line 2 ---------------------

tline=fgetl(fid);

tline=fgetl(fid);

comp = sscanf(tline, '%f');

for i=1:comp

tline=fgetl(fid);

tmparr = sscanf(tline, '%f');

tD(i)=tmparr(1);

tsig_a(i)=tmparr(2);

tsig_f(i)=tmparr(3);

end

%-------------------- # line 3 ---------------------

tline=fgetl(fid);

disp(tline);

tline=fgetl(fid);

tmp = sscanf(tline, '%d');

for i=1:M

c(i)=floor(tmp/10^(M-i));

tmp=tmp-c(i)*10^(M-i);

D(i)=tD(c(i));

sig_a(i)=tsig_a(c(i));

sig_f(i)=tsig_f(c(i));

str = sprintf('D(%d) = %f,  Sig_a(%d) = %f, nuSig_f(%d) = %f\n',i,D(i),i,sig_a(i),i,sig_f(i));

disp(str);

end

%-------------------- # line 4 ---------------------

tline=fgetl(fid);

disp(tline);

tline=fgetl(fid);

tmparr = sscanf(tline, '%e');

al_l=tmparr(1);

al_r=tmparr(2);

str=sprintf('albedo_left = %e , alebedo_right = %e\n',al_l,al_r);

disp(str);

%-------------------- # line 5 ---------------------

tline=fgetl(fid);

disp(tline);

tline=fgetl(fid);

n=sscanf(tline,'%d');

ntot=M*n;

str=sprintf('n = %d , tot# = %d\n',n,ntot);

disp(str);

%-------------------- # line 6 ---------------------

tline=fgetl(fid);

disp(tline);

tline=fgetl(fid);

tmparr = sscanf(tline, '%e');

conv_eig=tmparr(1);

conv_sol=tmparr(2);

str=sprintf('E_eig < %e , E_sol < %e\n', conv_eig, conv_sol);

disp(str);

end

2. Deriving Nodal balance equation

지난번 숙제와는 달리 이번에는 box scheme을 이용한다. Nodal balance equation의 유도는 아래의 diffusion equation부터 시작을 한다.

[eqn: - d over dx D(x) d over dx phi(x) + SIGMA_a (x) phi(x) = 1 over k nu SIGMA_f (x) phi(x)], [eqn: x in [0,H]] ... (1)

위의 식을 구간 [eqn: [x_i , x_i+1 ]]에서 적분을 해주는데 Fick's law [eqn: J(x) = -D(x) d over dx phi(x)]의 관계를 사용해주면 식은 다음과 같이 변형이 된다.

[eqn: int _ {x_i } ^{x_i+1 } {[ d over dx J(x) + SIGMA_a (x) phi(x) ] dx = int_ {x_i } ^{x_i+1 } {[ 1 over k nu SIGMA_f (x)phi(x)] dx]

[eqn: J_i+1 - J_i + SIGMA_a,i bar phi_i h_i = 1 over k nu SIGMA_f,i bar phi_i h_i] ... (2)

우리가 구하고자하는 값은 flux값 [eqn: phi]와 eigenvalue [eqn: lambda](=[eqn: 1 over k])이므로 current의 값 [eqn: J]를 다시 [eqn: phi]로 나타내줄 필요가 있다.

위의 그림은 current와 flux 사이의 notation을 나타낸 그림이다. 그림에서 보듯이 current [eqn: J]는 두 box mesh의 경계에서 정의가 된다. 그 값은 또 다시 경계에서 정의된 flux인 [eqn: phi_s]를 사용하여 다음과 같이 나타낼 수 있다.

[eqn: J_i+1 = - D_i {phi_s - bar phi_i }over{ h_i over 2} = - D_i+1 { {bar phi}_i+1 - phi_s } over {h_i+1 over 2}] ... (3)

i 번째 mesh로부터 구한 current와 i+1 번째 mesh로부터 구한 current는 같은 [eqn: J_i+1]이므로 [eqn: phi_s]에 대해서 정리하면 다음과 같다. (편의상 앞으로 [eqn: bar phi = phi]로 bar 표시를 생략한다.)

[eqn: phi_s = {beta_i phi_i + beta_i+1 phi_i+1 } over {beta_i + beta_i+1}] (단, [eqn: beta_i = D_i over h_i]) ... (4)

이를 다시 [eqn: J_i+1]식에 대입해주면 다음과 같다.

[eqn: J_i+1 = - 2 {beta_i beta_i+1 } over { beta_i + beta_i+1 } (phi _i+1 - phi_i )#```````````````````=- {tilde D} _i+1 (phi_i+1 - phi_i )], (여기서 [eqn: {tilde D } _ i+1 = {beta_i beta_i+1 } over { beta_i + beta_i+1 }]) ... (5)

이를 다시 식 (2)에 대입을 해주고 각각 [eqn: phi]에 대해서 정리를 하면 (단, 우변의 source항은 이항하지 않는다.)

[eqn: -{ tildeD }_i phi_i-1 + ( { tilde D} _i + {tilde D} _i+1 + SIGMA_a,i h_i ) phi_i - { tildeD }_i phi_i+1 = 1 over k nu SIGMA_f,i phi_i h_i] ... (6)

이것이 [eqn: i= 2, cdots , N-1]에서의 nodal balance equation이다. 왼쪽 끝과 오른쪽 끝에서는 [eqn: phi_N+1]와 [eqn: phi_0]가 존재하지 않으므로 식이 달라지는데 오른쪽 경계를 먼저 살펴본다. 다음은 오른쪽 경계 node에서의 current이다.

[eqn: J_N+1 = -D_N { phi_s - phi_N } over {h_N over 2}] ... (7)

양쪽 끝에서는 boundary condition들이 주어지는데 일반적으로는 albedo [eqn: alpha== J over phi]값으로 주어진다고 한다. 이 경우에는 [eqn: J_N+1 = alpha_R phi_s]로 B.C가 주어졌으므로 이와 식 (7)을 비교하여 다시 [eqn: phi_s]에 대해 정리를 해보면 다음과 같다.

[eqn: phi_s = {2beta_N } over {alpha_R +2 beta_N } phi _N] ... (8)

식 (8)을 식(7)에 다시 대입하면 다음과 같다.

[eqn: J_N+1 = {2alpha_R beta_N } over {alpha_R +2 beta_N } phi _N = {tildeD}_N+1 phi_N] (여기서 [eqn: {tilde D}_N+1 = {2alpha_R beta_N } over {alpha_R +2 beta_N }]) ... (9)

이를 식 (2)에 대입하면 식 (6)과 모양은 같으나 [eqn: phi_N+1]에 관한 항이 하나 빠진 다음과 같은 모양의 식이 된다.

[eqn: -{ tildeD }_N phi_N-1 + ( { tilde D} _N + {tilde D} _N+1 + SIGMA_a,N h_N ) phi_N = 1 over k nu SIGMA_f,N phi_N h_N] ... (10)

왼쪽경계에서도 오른쪽 경계와 마찬가지 방법으로 nodal balance equation을 구할 수 있다.

[eqn: ( { tilde D} _1 + {tilde D} _2 + SIGMA_a,1 h_1 ) phi_1 -{ tildeD }_2 phi_2 = 1 over k nu SIGMA_f,2 phi_2 h_2] ... (11)

이것으로 우리가 필요한 N개의 linear equation이 완성되었다.

3. Constructing linear system

2.에서 구한 N개의 식으로 [eqn: rm M PHI = 1 over k F PHI] 꼴의 eigenvalue equation을 만들 수 있다. 여기에는 M과 F는 N by N matrix이고 또 M은 tri-diagonal하다는 특징을 갖고 F는 diagonal하다.  그의 구성은 아래와 같다.

표  M

[eqn: {tilde D}_1 +{tilde D}_2 + SIGMA_a,1 h_1]	[eqn: -{tilde D}_2]
[eqn: -{tilde D}_2]	[eqn: {tilde D}_2 +{tilde D}_3 +SIGMA_a,2 h_2]	[eqn: ddots]
	[eqn: DDOTS]	[eqn: cdots]	[eqn: -{tilde D}_N-1]
		[eqn: -{tilde D}_N-1]	[eqn: {tilde D}_N-1 +{tilde D}_N + SIGMA_a,N-1 h_N-1]	[eqn: -{tilde D}_N]
			[eqn: -{tilde D}_N]	[eqn: {tilde D}_N +{tilde D}_N+1 +SIGMA_a,N h_N]

표  F

[eqn: nu SIGMA_f,1 h_1]
	[eqn: cdots]
		[eqn: nu SIGMA_f,N h_N]

빈칸은 모두 0이고 표의 크기는 달라보여도 둘 다 N by N 정사각행렬이다. 행렬 [eqn: M]의 기본적인 구성 방법을 소개하자면 먼저 각 성분에 들어갈 값들인 [eqn: beta]와 [eqn: tilde D]를 먼저 계산해 두고 경계조건이 필요하지 않은 2~(n-1) 번째 linear equation에 대한 구성요소들을 입력한 다음 첫 번째와 마지막 equation에서는 따로 각각 2개의 성분을 입력해주었다.

행렬 [eqn: F]는 diagonal 하기 때문에 [eqn: nu SIGMA_f h]만 계산해 두고 diagonal한 위치에 값을 입력해주면 되었다. 이를 구현한 코드는 다음과 같다.

function M=con_M

global M xi D sig_a al_l al_r n ntot

M=zeros(ntot);

h=xi/n;

%-------------------- ntot개의 점에서의 물성---------------------

%----- diffusivity beta와 sig_a*h 를 계산해둔다.-----

for i=1:ntot

j=ceil(i/n);

beta(i)=D(j)/h;

sig_ah(i)=sig_a(j)*h;

end

%----- tilde_D 를 계산한다.-----

for i=2:ntot

tild(i)=2*beta(i)*beta(i-1)/(beta(i)+beta(i-1));

end

%----- beta_0와 beta_N+1은 albedo로서 tilde_D를 따로 정의한다. -----

tild(1)=2*beta(1)*al_l/(2*beta(1)+al_l);

tild(ntot+1)=2*beta(ntot)*al_r/(2*beta(ntot)+al_r);

%--------------------------------------------------------------

%-------------------- ntot개의 equation들---------------------

%----- i=2 ... N-1 번째의 node에서의 equation -----

for i=2:ntot-1

M(i,i)=tild(i)+tild(i+1)+sig_ah(i);

M(i,i-1)=-tild(i);

M(i,i+1)=-tild(i+1);

end

%----- 첫번째와 N번째 node에서의 equation은 따로 정의한다.-----

M(1,1)=tild(1)+tild(2)+sig_ah(1);

M(1,2)=-tild(2);

M(ntot,ntot)=tild(ntot)+tild(ntot+1)+sig_ah(ntot);

M(ntot,ntot-1)=-tild(ntot);

%---------------------------------------------------------------

function F=con_F

global xi sig_f n ntot

F=zeros(ntot);

h=xi/n;

%-------------------- ntot개의 점에서의 물성---------------------

%----- sig_f*h 를 계산해둔다.-----

for i=1:ntot

j=ceil(i/n);

sig_fh(i)=sig_f(j)*h;

end

%-------------------- ntot개의 equation들---------------------

%----- diagonal한 성분을 채워 넣는다.-----

for i=1:ntot

F(i,i)=sig_fh(i);

end

4. LU - factorization

위의 tri-diagonal linear system을 LU-factorization 하는 함수를 작성하였다. LU factorization을 하면 주어진 행렬을 lower 삼각행렬, upper 삼각 행렬로 나누게 되는데 이는 [eqn: Ax=b]라는 행렬식을 풀 때, 계산수를 줄이는데 유용하다. 왜냐하면 [eqn: LUx=Ly=b]로부터 이미 elimination이 된 L을 통해 forward substitution을 하여 y를 쉽게 구할 수 있고 그를 통해 다시 [eqn: Ux=y]를 사용하여 backward substitution을 하여 x를 쉽게 구할 수 있기 때문이다.

안쪽의 loop에서는 M (multiplier) 를 구하고 대각행렬로 각 성분을 나누어 주면서 lower triangular matrix를 구하고  일을 하고 바깥쪽의 loop에서는 그로부터 행렬 A를 소거해가면서 upper triangular matrix를 만든다.

function [L,U]=luf(A)

n=length(A);

L=eye(n); % Lower triangle matrix

M=eye(n); % multiplyer

%-------------------- LU factorization ---------------------

for i=1:n

M=eye(n);

for j=i+1:n

L(j,i)=A(j,i)/A(i,i);

M(j,i)=-L(j,i);

end

A=M*A;  %소거되어가는 행렬 Upper triangle matrix

end

U=A;

%--------------------------------------------------------------

5. Solving by substitution

forward and backward substitution을 사용하여 1~4.를 통틀어서 part.1을 계산하는 함수를 작성하였다. 아래 함수의 실행 순서는 다음과 같다.

#1. .inp 파일을 읽는다.

#2. Linear system을 구축한다.

#3. 임의의 [eqn: phi^(0)]와 임의의 [eqn: k^(0)]를 설정한다.

#4. RHS ( [eqn: 1 over k F phi^(n)] )를 계산한다.

#5. [eqn: phi^(n+1)]을 구한다.

#6. [eqn: k^(n+1) = k^(n) { < phi ^(n+1) , phi^(n+1) > } over {< phi^(n) , phi^(n+1) >] 로 [eqn: k]를 업데이트 한다. 오차 또한 계산한다.

#7. 수렴조건을 판별하고 4~7을 반복한다.

function [k_n,phi_n,iter,rho]=solve

global M ntot conv_eig conv_sol

% #1 input

readinp();

% #2 construct Linear System

M=con_M();

F=con_F();

%--------------------initial value--------------------

% #3 intitial value

phi_o=ones(ntot,1); %initial eigen vector

k_o=1;        %initial eigen value

e_k=1;     % error of eigenvalue

e_s=1;     % error of flux(solution)

temp_k=1;   %eigenvalue error before one step

iter=0;   % # of iteration

%------------------------------------------------------

%-------------------- solving by LU & substitution--------------------

while e_k > conv_eig && e_s > conv_sol

iter=iter+1;

% #4 construct RHS ; B

B=F*phi_o/k_o;

% #5 LU-factorization & solve Ax=B

phi_n=solve_axb(M,B);

% #6 Update eigenvalue and flux

k_n=k_o*(phi_n'*phi_n)/(phi_o'*phi_n);

% #6-1 Calculate Error

e_k=abs(1-k_o/k_n);

phi_d=F*(phi_n-phi_o);

Q=F*phi_n;

e_s=abs((phi_d'*phi_d)/(Q'*Q))^0.5/ntot;

% #7 ready to do #4 to 7

k_o=k_n;

phi_o=phi_n;

rho=e_k/temp_k;

temp_k=e_k;

end

e_k

e_s

% solve AX=B by subsitution

function [X]=solve_axb(A,B)

n=length(A);

[L,U]=luf(A);

Y=zeros(n,1);

X=zeros(n,1);

Y(1)=B(1);

for i=2:n       %Y를 forward substitution으로 구함

s=0;

for j=1:i-1

s=s+L(i,j)*Y(j);

end

Y(i)=B(i)-s;

end

X(n)=Y(n)/U(n,n);

for i=n-1:-1:1  %X를 backward substitution으로 구함

s=0;

for j=i+1:n

s=s+U(i,j)*X(j);

end

X(i)=(Y(i)-s)/U(i,i);

end

nd1d_1.inp에 대한 실행결과는 다음과 같다.

e_k =

2.4666e-006

e_s =

8.9962e-006

rho =

0.8374

iter =

31

k =

1.0217

error는 eigenvalue의 error가 더 빨리 줄어들었음을 확인할 수 있었고 eigenvalue는 약 1.0217가 나왔다. iteration은 31번 반복되었고 dominance ratio [eqn: rho]는 0.84 정도가 나왔다.

Part 2. Power Iteration and Wielandt Shift

1. Power iteration

직접적인 계산을 통한 풀이가 아닌 반복을 통한 해에 접근해나가는 Power iteration을 하는 코드이다. Power iteration은 [eqn: Ax=lambda x]꼴의 eigenvalue problem을 풀 때, initial guess  solution에 행렬 [eqn: A]를 계속 곱해가면서 [eqn: x ^{(l+1)} =Ax ^{(l)} = lambda x ^{(l)}]로부터 eigenvalue를 구하는 것이다. 위의 식은 [eqn: l]이 매우 커질 때 수렴하는데 그 수렴은 다음과 같이 증명 할 수 있다.

초기값 [eqn: x^(0)]는 행렬 [eqn: A]의 eigen vector인 [eqn: u_i]들의 합으로 다음과 같이 나타낼 수 있다.

[eqn: x^(0) = c_1 u_1 +c_2 u_2 + cdots + c_N u_N]

여기에 행렬 [eqn: A]를 계속 곱해가고 있다고 생각 할 때 다음과 같이 되고

[eqn: x^(l+1) = A ^l phi^(0) = A^l ( c_1 u_1 + c_2 u_2 + cdots + c_N u_N )]

[eqn: = c_1 k_1 ^l u_1 + c_2 k_2 ^l u_2 + cdots + c_N k_N ^l u_N]

[eqn: = k_1 ^l ( c_1 u_1 + SUM_i=2 ^ N [ ( {k_i over k_1 )^l c_i u_i ] )] [eqn: SIMEQ k_1 ^l c_1 u_1]

[eqn: l]이 커질수록 즉, 반복횟수가 증가할수록 [eqn: k_1 > k_i] ([eqn: i>=2]) 이므로 [eqn: ( k_i over k_1 )^l]이 zero로 수렴해가므로 위의 식으로부터 [eqn: k_1] 즉, first eigenvalue를 구할 수 있다.

여기서는 Wielandt shift도 포함하고 있는데 일단 여기서는 shift값을 zero를 넣어주면 (즉 [eqn: k_s -> inf]) 그것이 그냥 power method가 된다. 이 함수의 실행순서는 다음과 같다.

#1. .inp 파일을 읽는다.

#2. Linear system을 구축한다. 이 때 [eqn: k_s]만큼 eigenvalue를 shift 한다. (지금은 inf값)

#3. 임의의 [eqn: phi^(0)]와 임의의 [eqn: k^(0)]를 설정한다.

#4. [eqn: phi^(n)]과 [eqn: k^(n)](scaling을 위한)으로 RHS [eqn: Q]를

#5. [eqn: M`' phi^(n+1) = Q]로부터 [eqn: phi^(n+1)]을 구한다.

#6. [eqn: k^(n+1) = k^(n) { < phi ^(n+1) , phi^(n+1) > } over {< phi^(n) , phi^(n+1) >] 로 [eqn: k]를 업데이트 한다. 오차 또한 계산한다.

#7. eigenvalue와 flux에 대해서 수렴조건을 판별하고 4~7을 반복한다.

function [k,e,phi_n,iter,rho]=pow_m()

global ntot conv_eig conv_sol

% #1 input

readinp();

%--------------------Wielandt Shift---------------------

% #2 construct Linear System

M=con_M();

F=con_F();

k_s=1.5;  %<- k_s 만큼의 shift를 준다. inf for power_method

%M=M-1/k_s*F; %<-----------------shifted------------------

%------------------------------------------------------

%--------------------initial value--------------------

% #3 intitial values

phi_o=ones(ntot,1);  %initial eigen vector

k_o=1;        %initial eigen value

temp_k=1;     %eigenvalue error before one step

iter=0;       %# of iteration

e_k=1;    % error of eigenvalue

e_s=1;    % error of flux(solution)

%-----------------------------------------------------

%--------------------Power Iteration---------------------

psi_o = F*phi_o;  % Fission source

while e_k > conv_eig || e_s > conv_sol

iter=iter+1;      % 반복 횟수

% #4 use psi_o to get Q

Q=1/k_o *psi_o; % Fission source scaling

% #5 M phi = Q => get new phi

phi_n=solve_axb(M,Q);

psi_n=F*phi_n;  % new psi

% #6 get new eigenvalue

k_n=k_o*((psi_n'*psi_n)/(psi_o'*psi_n));

% #6-1 Calculate Error

e_k=abs(1-k_o/k_n);

psi_d=F*(psi_n-psi_o);

e_s=abs((psi_d'*psi_d)/(psi_n'*psi_n))^0.5;

% #7 ready to do #4 to 7

k_o=k_n;

rho=e_k/temp_k;

temp_k=e_k;

psi_o=psi_n;

end

% #8 re-shift eigenvalue

k=k_n*k_s/(k_n+k_s);

x=1:ntot;

plot(x,phi_n);

e_k

e_s

iter

e=abs(1-k/analy);

2. Dependence of the eigenvalue error (compare with analytic eigenvalue) on the node size

nd1d_1.inp 파일에서 type을 22222로 바꾼 후에 mesh의 수를 2배씩 증가시키면서 해석적인 해 [eqn: k={ nu SIGMA_f } over {SIGMA_a + ( pi over H )^2 D]와의 오차율을 계산하고 그 경향성을 살펴보았다. 아래가 그 결과 표이다.

# of mesh in one domain	error ratio
5	1.28496e-005
10	3.19629e-006
20	7.22870e-007
40	1.10540e-007

그림  [eqn: n]-[eqn: e] 그래프

그림  [eqn: 1 over n^2]-[eqn: e] 그래프

지난숙제에서 discretization에 의해서 발생하는 오차는 mesh 크기의 제곱에 비례한다는 것을 보였다. 그 결과와 같이 이번 과제에서도 [eqn: 1over n^2 PROPTO (DELTA h )^2]와 [eqn: e](error)가 서로 비례함을 확인할 수 있었다.

3. Wielandt Shift

dominance ratio [eqn: sigma = lambda_2 overlambda_1]으로 정의되는데 power method의 수렴속도는 이 값에 의해 결정된다. 따라서 이 값을 줄이면 줄일수록 수렴속도가 매우 빨라지게 된다. Wielandt shift acceleration은 eigenvalue를 일정값만큼 평행이동 시켜서 이 [eqn: sigma]값을 줄이는 방법이다. 이동시키는 양 [eqn: lambda_s]는 [eqn: lambda_2]에 가까우면 가까울수록 [eqn: sigma_s = { lambda_2 - lambda_s } over { lambda_1 - lambda_s]가 0에 매우 가까워지게 된다. 그러나 [eqn: lambda_2]는 알고 있는 값이 아니기 때문에 적절한 [eqn: lambda_s]를 택하여 주는 것도 필요하다. 여기서는 [eqn: lambda]의 역수 [eqn: k]를 eigenvalue로서 구하고 있어서 표기가 살짝 다르다. Wielandt shift를 구현하는 방법은 아래와 같다.

[eqn: M phi = 1 over k F phi]라고 할 때, [eqn: k_s]만큼 eigenvalue를 shift 시켜준다고 하면 다음과 같다

[eqn: (M - 1 over k_s F ) phi = ( 1 over k - 1 over k_s ) F phi]

[eqn: M``' phi = 1 over k' F phi]

F값은 바뀌지 않고 M 행렬만 구성이 바뀐다. shift 전의 dominance ratio를 [eqn: sigma = k_2 over k_1]이라고 하면 shift 후의 값은 [eqn: sigma_s = { k_s - k_1 } over { k_s - k_2 } sigma]로 줄게 된다. shift 한 후에 [eqn: k ']을 구했으면 다시 [eqn: -k_s]만큼 평행 이동시켜 [eqn: k]를 얻게 된다. 여기서 shift 값 [eqn: 1 over k_s (= lambda_s )]를 0를 넣어주면 power method와 같다. 코드에서는 [eqn: k_s]로 shift 해주었기 때문에 [eqn: inf]값을 넣어주어야 그냥 power method가 된다.

4. Compare performance of the Wielandt shift and non-Wielandt shift based on dominance ratio

[eqn: k_s]	dominance ratio [eqn: rho]
inf (power method)	7.444260398466646e-001
5	6.946115658004977e-001
4	6.790024571266884e-001
3	6.491648222735723e-001
2	5.695117850106725e-001
1.5	4.451946025789574e-001
1.1	9.554227854039632e-002

[eqn: k_s]값이 커질수록 dominance ratio [eqn: rho]값은 점점 작아졌다. 1.1 까지 시도해보았는데 이 경우 일반적인 power method와 비교하면 약 8배 정도 빠름을 알 수 있었다.

5. Ploting flux

그림  inp_1

그림  inp_2

그림  inp_3

그림  inp_4

그림  inp_5

그림  inp_6