-
Notifications
You must be signed in to change notification settings - Fork 1
/
create_left_reflecting_bc_matrix.m
95 lines (75 loc) · 2 KB
/
create_left_reflecting_bc_matrix.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
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
function bmat = create_left_reflecting_bc_matrix(diam,n,lambda,alpha,model)
% Create the iteration matrix for a one-sided, positive tempered fractional
% diffusion equation assuming reflecting BCs.
% n = number of gridpoints
% 1 < alpha <= 2
% model is either 'norm' for 2-term normalized TFD eqn or 'cent' for 3-term
% centered normalized TFD eqn
if (alpha >2 || alpha <= 1)
error('alpha must be greater than one and less than or equal to two')
end
% switch between normalized and centered normalized models
if (strcmp(model,'norm'))
cflg = 0;
elseif (strcmp(model,'cent'))
cflg = 1;
end
h = diam/n;
lh = lambda*h;
np1 = n +1;
gwgts = create_grunwald_weights(alpha,n);
alpham1 = alpha-1;
kappa = zeros(1,n+1);
for i = 0:n
i1 = i+1;
kappa(i1) = sum(gwgts(1:i1)'.*exp(-lh*(-1:i-1)));
end
bmat = zeros(np1);
for i = 0:n
for j = 1:(n-1)
i1 = i + 1;
j1 = j + 1;
if (i <= (j+1))
k = j-i+1;
k1 = k+1;
bmat(i1,j1) = gwgts(k1)*exp(-lh*(j-i));
end
if (i == j)
bmat(i1,j1) = gwgts(2) - exp(lh)*(1-exp(-lh))^alpha - ...
cflg*alpha*(lh)^alpham1;
end
if (i == (j-1))
bmat(i1,j1) = gwgts(3)*exp(-lh) + cflg*alpha*(lh)^alpham1;
end
end
if (i < n-1)
j = n;
i1 = i + 1;
j1 = j + 1;
k = n - i;
k1 = k+1;
bmat(i1,j1) = exp(lh)*(((1-exp(-lh))^alpha))-kappa(k1);
end
end
% special cases
i = 0; j = 0;
i1 = i+1; j1 = j+1;
bmat(i1,j1) = gwgts(2) - exp(lh)*(1-exp(-lh))^alpha + gwgts(1)*exp(lh) - ...
cflg*alpha*(lh)^alpham1;
i = 1; j = 0;
i1 = i+1; j1 = j+1;
bmat(i1,j1) = gwgts(1)*exp(lh);
i = 0; j = n;
i1 = i+1; j1 = j+1;
k = j-i;
k1 = k+1;
bmat(i1,j1) = exp(lh)*(((1-exp(-lh))^alpha)) - kappa(k1);
i = n-1; j = n;
i1 = i+1; j1 = j+1;
k = n-i;
k1 = k+1;
bmat(i1,j1) = exp(lh)*(1-exp(-lh))^alpha - kappa(k1) + cflg*alpha*(lh)^alpham1;
i = n; j = n;
i1 = i+1; j1 = j+1;
bmat(i1,j1) = -gwgts(1)*exp(lh);
end