LINE Solver
MATLAB API documentation
Loading...
Searching...
No Matches
mam_transient2.m
1function V = mam_transient2(B, L, F, Lv, T, n, m, s)
2%MAM_TRANSIENT2 Laplace-domain transient V(s,n,m) for a closed (finite) QBD.
3% V = MAM_TRANSIENT2(B, L, F, Lv, T, n, m, s) returns the Laplace transform
4% (at complex argument s) of the transient transition-probability matrix from
5% level n to level m of a finite piecewise level-dependent QBD with regime
6% thresholds T; the top level is T(end).
7%
8% Block cell arrays follow the same convention as MAM_TRANSIENT2_OPEN, with
9% K = length(T)-1 regimes and Lv{K+1} the top boundary level.
10%
11% Ported from the transient-QBD research code (Horvath et al. formulation).
12%
13% Copyright (c) 2012-2026, Imperial College London
14% All rights reserved.
15 K = length(T) - 1;
16
17 Gs = cell(K, 1);
18 Rs = cell(K, 1);
19 Ghs = cell(K, 1);
20 Rhs = cell(K, 1);
21 for k = 1:K
22 if T(k+1) - T(k) == 1
23 Gs{k} = [];
24 Rs{k} = [];
25 Ghs{k} = [];
26 Rhs{k} = [];
27 else
28 Ik = eye(size(L{k}, 1));
29 [G, R] = qbd_fundmat(B{k}, L{k} - s*Ik, F{k}, 'GR');
30 Gs{k} = G; Rs{k} = R;
31 [G, R] = qbd_fundmat(F{k}, L{k} - s*Ik, B{k}, 'GR');
32 Ghs{k} = G; Rhs{k} = R;
33 end
34 end
35
36 SvHn = cell(K, 1);
37 SvH0 = cell(K, 1);
38 SvHhn = cell(K, 1);
39 SvHh0 = cell(K, 1);
40 for k = 1:K
41 NN = size(Lv{k}, 1);
42 if T(k+1) - T(k) > 1
43 d = T(k+1) - T(k);
44 Ik = eye(NN);
45 num = [mxpow(Ghs{k}, d-1), Gs{k}; Ghs{k}, mxpow(Gs{k}, d-1)];
46 den = [Ik, mxpow(Gs{k}, d);
47 mxpow(Ghs{k}, d), Ik];
48 SH = num / den;
49 SvHn{k} = SH(1:NN, 1:NN);
50 SvH0{k} = SH(1:NN, NN+1:2*NN);
51 SvHhn{k} = SH(NN+1:2*NN, 1:NN);
52 SvHh0{k} = SH(NN+1:2*NN, NN+1:2*NN);
53 else
54 NN1 = size(Lv{k+1}, 1);
55 SvH0{k} = zeros(NN1, NN);
56 SvHh0{k} = eye(NN);
57 SvHhn{k} = zeros(NN, NN1);
58 SvHn{k} = eye(NN1);
59 end
60 end
61
62 NNK = size(Lv{K+1}, 1);
63 SY = cell(K, 1);
64 SY{K} = SvH0{K} + SvHn{K} * ((s*eye(NNK) - Lv{K+1} - B{K}*SvHhn{K}) \ (B{K} * SvHh0{K}));
65 for k = K-1:-1:1
66 NNk1 = size(Lv{k+1}, 1);
67 SY{k} = SvH0{k} + SvHn{k} * ((s*eye(NNk1) - Lv{k+1} - F{k+1}*SY{k+1} - B{k}*SvHhn{k}) \ (B{k} * SvHh0{k}));
68 end
69
70 NN1 = size(Lv{1}, 1);
71 SYh = cell(K, 1);
72 SYh{1} = SvHhn{1} + SvHh0{1} * ((s*eye(NN1) - Lv{1} - F{1}*SvH0{1}) \ (F{1} * SvHn{1}));
73 for k = 2:K
74 NNk = size(Lv{k}, 1);
75 SYh{k} = SvHhn{k} + SvHh0{k} * ((s*eye(NNk) - Lv{k} - B{k-1}*SYh{k-1} - F{k}*SvH0{k}) \ (F{k} * SvHn{k}));
76 end
77
78 SV = cell(K+1, K+1);
79 for l = 0:K
80 if l == 0
81 SV{l+1, l+1} = (s*eye(NN1) - Lv{1} - F{1}*SY{1}) \ eye(NN1);
82 elseif l == K
83 SV{l+1, l+1} = (s*eye(NNK) - Lv{K+1} - B{K}*SYh{K}) \ eye(NNK);
84 else
85 NNl1 = size(Lv{l+1}, 1);
86 SV{l+1, l+1} = (s*eye(NNl1) - Lv{l+1} - F{l+1}*SY{l+1} - B{l}*SYh{l}) \ eye(NNl1);
87 end
88 for k = l+1:K-1
89 NNk1 = size(Lv{k+1}, 1);
90 SV{k+1, l+1} = (s*eye(NNk1) - Lv{k+1} - F{k+1}*SY{k+1} - B{k}*SvHhn{k}) \ (B{k} * SvHh0{k} * SV{k, l+1});
91 end
92 if l < K
93 SV{K+1, l+1} = (s*eye(NNK) - Lv{K+1} - B{K}*SvHhn{K}) \ (B{K} * SvHh0{K} * SV{K, l+1});
94 end
95 for k = l-1:-1:1
96 NNk1 = size(Lv{k+1}, 1);
97 SV{k+1, l+1} = (s*eye(NNk1) - Lv{k+1} - F{k+1}*SvH0{k+1} - B{k}*SYh{k}) \ (F{k+1} * SvHn{k+1} * SV{k+2, l+1});
98 end
99 if l > 0
100 SV{1, l+1} = (s*eye(NN1) - Lv{1} - F{1}*SvH0{1}) \ (F{1} * SvHn{1} * SV{2, l+1});
101 end
102 end
103
104 pos = find(T > n, 1);
105 if isempty(pos), kn = K+1; else, kn = pos - 1; end
106 pos = find(T > m, 1);
107 if isempty(pos), km = K+1; else, km = pos - 1; end
108
109 NN = size(Lv{kn}, 1);
110 II = eye(NN);
111
112 if T(kn) == n
113 if T(km) == m
114 V = SV{kn, km};
115 return;
116 end
117 Vu = SV{kn, km+1};
118 Vl = SV{kn, km};
119 Lu = T(km+1);
120 Ll = T(km);
121 else
122 d1 = T(kn+1) - n;
123 num = [mxpow(Ghs{kn}, d1-1), Gs{kn}; Ghs{kn}, mxpow(Gs{kn}, d1-1)];
124 den = [II, mxpow(Gs{kn}, d1); mxpow(Ghs{kn}, d1), II];
125 Tmp = num / den;
126 HTnn = Tmp(1:NN, 1:NN);
127 HTn0 = Tmp(1:NN, NN+1:2*NN);
128 HhTnn = Tmp(NN+1:2*NN, 1:NN);
129 HhTn0 = Tmp(NN+1:2*NN, NN+1:2*NN);
130
131 d2 = n - T(kn);
132 num = [mxpow(Ghs{kn}, d2-1), Gs{kn}; Ghs{kn}, mxpow(Gs{kn}, d2-1)];
133 den = [II, mxpow(Gs{kn}, d2); mxpow(Ghs{kn}, d2), II];
134 Tmp = num / den;
135 HnTn = Tmp(1:NN, 1:NN);
136 HnT0 = Tmp(1:NN, NN+1:2*NN);
137 HhnTn = Tmp(NN+1:2*NN, 1:NN);
138 HhnT0 = Tmp(NN+1:2*NN, NN+1:2*NN);
139
140 if kn == K
141 Yn = HTn0 + HTnn * ((s*eye(NNK) - Lv{K+1} - B{K}*HhTnn) \ (B{K} * HhTn0));
142 else
143 NNkn1 = size(Lv{kn+1}, 1);
144 Yn = HTn0 + HTnn * ((s*eye(NNkn1) - Lv{kn+1} - F{kn+1}*SY{kn+1} - B{kn}*HhTnn) \ (B{kn} * HhTn0));
145 end
146 if kn == 1
147 Yhn = HhnTn + HhnT0 * ((s*eye(NN1) - Lv{1} - F{1}*HnT0) \ (F{1} * HnTn));
148 else
149 NNkn = size(Lv{kn}, 1);
150 Yhn = HhnTn + HhnT0 * ((s*eye(NNkn) - Lv{kn} - B{kn-1}*SYh{kn-1} - F{kn}*HnT0) \ (F{kn} * HnTn));
151 end
152
153 Mkn = s*eye(size(L{kn},1)) - L{kn};
154 if T(km) < n
155 Vnl = (Mkn - B{kn}*HhnTn - F{kn}*Yn) \ (B{kn} * HhnT0 * SV{kn, km});
156 else
157 Vnl = (Mkn - F{kn}*HTn0 - B{kn}*Yhn) \ (F{kn} * HTnn * SV{kn+1, km});
158 end
159 if m == T(km)
160 V = Vnl;
161 return;
162 end
163 if T(km+1) < n
164 Vnu = (Mkn - B{kn}*HhnTn - F{kn}*Yn) \ (B{kn} * HhnT0 * SV{kn, km+1});
165 else
166 Vnu = (Mkn - F{kn}*HTn0 - B{kn}*Yhn) \ (F{kn} * HTnn * SV{kn+1, km+1});
167 end
168 Vnn = (s*II - L{kn} - B{kn}*Yhn - F{kn}*Yn) \ II;
169 if n == m
170 V = Vnn;
171 return;
172 end
173 if kn ~= km
174 Vu = Vnu; Vl = Vnl; Lu = T(km+1); Ll = T(km);
175 elseif n <= m
176 Vu = Vnu; Vl = Vnn; Lu = T(km+1); Ll = n;
177 else
178 Vu = Vnn; Vl = Vnl; Lu = n; Ll = T(km);
179 end
180 end
181
182 NN = size(Rs{km}, 1);
183 II = eye(NN);
184 Zden = [II, mxpow(Rs{km}, Lu-Ll); mxpow(Rhs{km}, Lu-Ll), II];
185 Znum = [mxpow(Rs{km}, m-Ll); mxpow(Rhs{km}, Lu-m)];
186 Z = Zden \ Znum;
187 V = Vl * Z(1:NN, :) + Vu * Z(NN+1:2*NN, :);
188end
the
Definition Station.m:245
Generated by doxygen 1.9.8