-
Notifications
You must be signed in to change notification settings - Fork 0
/
q11.m
65 lines (51 loc) · 1.73 KB
/
q11.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
%% define constants
% total number of particles
N = 100 ;
a = 1 ;
lambda1 = 5 ;
% number of states
k = 25 ;
% vector of states used for plots later
states_vector = 0:k-1 ;
%% solve q10 odes
initial_condition = 4.* ones(1,k) ;
[t_meanfield, X_meanfield] = ode15s(@(t,X)RHS_meanfield_q10(t,X,lambda1,a,k),[0,1000],initial_condition) ;
index_stable = length(X_meanfield) ;
%% steady state in q9
k_0 = 0 ; % first state in which initial condition non-zero is n_0
q9steadystates = zeros(1,25) ;
% create a vector of theoretical distribution of states given in Question 9
for i = 0:k-1
q9steadystates(i + 1) = scaledpoisson(N, lambda1, a, i, k_0) ;
end
%% make plots
f1 = figure ;
figure(f1)
plot(states_vector, q9steadystates, 'linewidth', 5)
hold on
plot(states_vector, X_meanfield(index_stable , :), 'Color', 'r', 'LineStyle', '--', 'linewidth' , 3)
xlabel('Particle state')
ylabel('Number of particles')
title('Steady states: Theoretical distribution vs ODE solutions')
legend('Theoretical distribution', 'Solution of ODEs')
%% define scaled poisson eq'n
function steadystate = scaledpoisson(N,lambda1, a, k, k_0)
steadystate = N*exp(-lambda1/a)*(lambda1/a)^(k-k_0)*(1/factorial(k-k_0)) ;
end
% define mean-field equations from q10
function dX = RHS_meanfield_q10(t,X,lambda1,a,k)
dX(1) = a.*X(1).*khat(X,k) - lambda1.*X(1) ;
for i = 2:k-1
dX(i) = a.*X(i).*(khat(X,k) - (i-1)) + lambda1.*X(i-1) - lambda1.*X(i) ;
end
dX(k) = a.*X(k).*(khat(X,k) - (k-1)) + lambda1.*X(k-1) ;
dX = dX' ;
end
% define function that calculates k^
function khat = khat(X,k)
cumsum = 0 ;
for i = 1:k
cumsum = cumsum + (i-1).*X(i) ;
end
khat = cumsum / sum(X) ;
end