-
Notifications
You must be signed in to change notification settings - Fork 1
/
matrixs.m
58 lines (41 loc) · 1.36 KB
/
matrixs.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
function [varargout]= matrixs(p)
Ap=zeros(p.Np-1);
An=zeros(p.Nn-1);
Bp= zeros(p.Np-1,1);
Bn= zeros(p.Nn-1,1);
%% Positive and Negative solid particle matrices
sp=p.Ds_p/(p.delta_p^2);
sn=p.Ds_n/(p.delta_n^2);
for k=1:p.Np-1
%diagonal
Ap(k,k)=-2*sp;
An(k,k)=-2*sn;
%superdiagonal
if (k~= (p.Np-1))
Ap(k,k+1)= (k+1)/k.*sp;
An(k,k+1)= (k+1)/k.*sn;
end
%subdiagonal
if (k~= 1)
Ap(k,k-1)= (k-1)/k.*sp;
An(k,k-1)= (k-1)/k.*sn;
end
end
%1st order boundary discretization
% An(p.Nn-1,p.Nn-1)=(p.Ds_n/((p.Nn-1)*p.delta_n*(p.delta_n^2)))*(p.delta_n-(p.Nn-1)*p.delta_n); %equivalently sn*(2-p.Nn)/(p.Nn-1);
% Ap(p.Np-1,p.Np-1)=(p.Ds_p/((p.Np-1)*p.delta_p*(p.delta_p^2)))*(p.delta_p-(p.Np-1)*p.delta_p); %equivalently sp*(2-p.Np)/(p.Np-1);
%
% Bp(end,end)= -(p.Np)/(p.Np-1) /p.delta_p;
% Bn(end,end)= -(p.Nn)/(p.Nn-1) /p.delta_n;
%2nd order boundary discretization
An(p.Nn-1,p.Nn-1)=(p.Ds_n*(6-2*p.Nn))/(3*(p.Nn-1)*p.delta_n^2);
An(p.Nn-1,p.Nn-2)=-(p.Ds_n*(6-2*p.Nn))/(3*(p.Nn-1)*p.delta_n^2);
Ap(p.Np-1,p.Np-1)=(p.Ds_p*(6-2*p.Np))/(3*(p.Np-1)*p.delta_p^2);
Ap(p.Np-1,p.Np-2)=-(p.Ds_p*(6-2*p.Np))/(3*(p.Np-1)*p.delta_p^2);
Bp(end,end)= -(p.Np)/(p.Np-1) * (2)/(3*p.delta_p);
Bn(end,end)= -(p.Nn)/(p.Nn-1) * (2)/(3*p.delta_n);
varargout{1}=Ap;
varargout{2}=An;
varargout{3}=Bn;
varargout{4}=Bp;
end