MF_point_fat

PHOTO EMBED

Tue May 24 2022 15:44:26 GMT+0000 (Coordinated Universal Time)

Saved by @tino #python

clc; clear; close all;
% known parameter
c=0.95;
siga=0.02;
sigt=siga/(1-c);
D=1/(3*sigt);
nu_sigf=0.0165;
N=100;
H=100;
h=H/N;
% calculate parameter
sigd=D/h^2;
ave_siga=siga;
ave_nu_sigf=nu_sigf;
alpha_l=10^7; %black
alpha_r=alpha_l; %black
% set matrix for M and F (linear system)
M=zeros(N+1);
F=zeros(N+1);
M(1,1)=sigd+(alpha_l/h)+(0.5*siga);
M(N+1,N+1)=M(1,1);
F(1,1)=0.5*nu_sigf;
F(N+1,N+1)=F(1,1);
for i=2:N
M(i-1,i)=-sigd;
M(i,i-1)=-sigd;
M(i,i)=2*sigd+ave_siga;
F(i,i)=ave_nu_sigf;
end

%%% solve matrix
%%% using eig and inv
A=inv(M)*F;
A
[flux,k]=eig(A)
flux_eig1=flux(:,1);
k_eig1=k(1,1);
flux_eig2=flux(:,2);
k_eig2=k(2,2);
flux_eig3=flux(:,3);
k_eig3=k(3,3);
k_eig=[k_eig1 k_eig2 k_eig3]
content_copyCOPY