-
Notifications
You must be signed in to change notification settings - Fork 0
/
solve_splitted_ode.m~
107 lines (92 loc) · 3.19 KB
/
solve_splitted_ode.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
95
96
97
98
99
100
101
102
103
104
105
106
107
function [time_samples, transient_prob]= solve_splitted_ode(c1, c2, c3, x1_0, x2_0, interval, solver)
% Q_exact is the matrix of orignal Markov chain
% Q_splitted_T is the collection of matrices that each of them is the
% transpose of generator matrix of a splitted markov chain
Q_splitted_T = cell(x2_0, 1);
time_samples = cell(x2_0, 1);
tr_probs = cell(x2_0, 1);
for i = 0:x2_0 - 1,
x = [x1_0 x2_0 - i 0 i];
nState4ithStage = min([x1_0 x2_0 - i]) + 1;
nTransitions = (nState4ithStage - 1) * 2;
transition_index = 1;
state_index = 1;
I = zeros(1, nState4ithStage + nTransitions);
J = zeros(1, nState4ithStage + nTransitions);
NZ = zeros(1, nState4ithStage + nTransitions);
while x(1) >= 0 && x(2) >= 0,
exit_rate = 0;
% R1
rate = c1 * x(1) * x(2);
if rate ~= 0,
I(transition_index) = state_index;
J(transition_index) = state_index + 1;
NZ(transition_index) = rate;
exit_rate = exit_rate + rate;
transition_index = transition_index + 1;
end
% R2
rate = c2 * x(3);
if rate ~= 0,
I(transition_index) = state_index;
J(transition_index) = state_index - 1;
NZ(transition_index) = rate;
exit_rate = exit_rate + rate;
transition_index = transition_index + 1;
end
% Add diagonal element
if exit_rate ~= 0,
I(transition_index) = state_index;
J(transition_index) = state_index;
NZ(transition_index) = - exit_rate;
transition_index = transition_index + 1;
end
state_index = state_index + 1;
x(1) = x(1) - 1;
x(2) = x(2) - 1;
x(3) = x(3) + 1;
end
Q_splitted_T{i + 1} = (sparse(I, J, NZ, nState4ithStage, nState4ithStage))';
% solve ode for the currrent split
pi0 = zeros(nState4ithStage, 1);
pi0(1) = 1;
[time_sample_i, tr_prob_i] = solver (@deriva_i, interval, pi0);
time_samples{i + 1} = time_sample_i;
tr_probs{i + 1} = tr_prob_i;
% Test
[tt,expexp] = app_exp_of_x3_in_birth_death_process
% Test
end
pi0 = zeros(x2_0 + 1, 1);
pi0(1) = 1;
[time_samples, transient_prob] = solver(@deriva_main, interval, pi0);
function pdot = deriva_i(t, p)
pdot = Q_splitted_T{i + 1} * p;
end
function pdot = deriva_main(t, p)
pdot = zeros(x2_0 + 1, 1);
pdot(1) = - calc_lambda(t, 1) * p(1);
for j = 2:x2_0,
pdot(j) = calc_lambda(t, j - 1) * p(j - 1) - calc_lambda(t, j) * p(j);
end
pdot(x2_0 + 1) = calc_lambda(t, x2_0) * p(x2_0);
end
function lambda_i_t = calc_lambda(t, i)
% calculate outgoing lambda from state i at time t
if i < 1 || i > x2_0,
lambda_i_t = 0;
else
len = length(time_samples{i});
for k = 1:len - 1,
if t >= time_samples{i}(k) && t <= time_samples{i}(k + 1),
tr_probs_t = (tr_probs{i}(k + 1,:) - tr_probs{i}(k,:)) / (time_samples{i}(k + 1) - time_samples{i}(k)) * (t - time_samples{i}(k)) + tr_probs{i}(k,:);
lambda_i_t = c3 * sum((0:size(tr_probs{i}, 2) - 1).*tr_probs_t);
break;
end
end
if ~(k < len),
fprintf('ERROR, No interpolation interval found.\n');
end
end
end
end