Solving ODE with big dimension
1 Ansicht (letzte 30 Tage)
Ältere Kommentare anzeigen
I have an ODE code to run an SEIR type epidemiological system segemented by age (n=21), and social group (soc=10). I expect the output to be of the form Classes.S(:,:,i) where i=1:10. The code is runing but only the first row of my initail condition is repeated 10times.
This is the code snippet:
function [Classes] = ODE_social_TEST2(para,ICs,maxtime,Mix_H,Mix_W,Mix_S,Mix_O,n,soc,w)
opts = odeset('RelTol',1e-4,'AbsTol',1e-3); %tolerance level%tolerance level
tspan = [0:1:maxtime]; %tie span
%Shape of the initial condition
All = reshape([ICs.S'; ICs.E'; ICs.A'; ICs.J';ICs.Q';ICs.I';ICs.H'; ICs.R'; ICs.D'; ICs.Ir'], []*n,1);
[t , pop] = ode15s(@(t,y)diff_socialmodel(t,y,para),tspan,All,opts);
%% Define the structure of the output
Classes = struct('S',zeros(numel(t), n, 1),'E',zeros(numel(t), n, 1),'A',zeros(numel(t), n, 1),'J',zeros(numel(t), n, 1),'Q',zeros(numel(t), n, 1), ...
'I',zeros(numel(t), n, 1),'H',zeros(numel(t), n, 1),'R',zeros(numel(t), n, 1),'D',zeros(numel(t), n, 1),'Ir',zeros(numel(t), n, 1),'t',[]);
for i = 1:soc
Classes.S(:,:,i) = pop(:,1:1:n);
Classes.E(:,:,i) = pop(:,n+1:1:2*n);
Classes.A(:,:,i) = pop(:,(2*n)+1:1:3*n);
Classes.J(:,:,i) = pop(:,(3*n)+1:1:4*n);
Classes.Q(:,:,i) = pop(:,(4*n)+1:1:5*n);
Classes.I(:,:,i) = pop(:,(5*n)+1:1:6*n);
Classes.H(:,:,i) = pop(:,(6*n)+1:1:7*n);
Classes.R(:,:,i) = pop(:,(7*n)+1:1:8*n);
Classes.D(:,:,i) = pop(:,(8*n)+1:1:9*n);
Classes.Ir(:,:,i) = pop(:,(9*n)+1:1:10*n);
end
Classes.t = t;
%% Set up the population change
function dpop = diff_socialmodel(t,pop,para)
%% DEfine the age structured mixing
ageMix = Mix_H + Mix_W + Mix_S + Mix_O;
S = pop(1:1:n,:);
E = pop(n+1:1:2*n,:);
A = pop((2*n)+1:1:3*n,:);
J = pop((3*n)+1:1:4*n,:);
Q = pop((4*n)+1:1:5*n,:);
I = pop((5*n)+1:1:6*n,:);
H = pop((6*n)+1:1:7*n,:);
R = pop((7*n)+1:1:8*n,:);
D = pop((8*n)+1:1:9*n,:);
Ir = pop((9*n)+1:1:10*n,:);
% Set-up the size for the population
dpop = zeros(size(pop));
%% Run the logistic curve
y = GeneralisedRichard(maxtime);
for m =1:length(y)
[pi(m)] = y(m)/sum(para.N);
end
% Age mixing
AgeMat = (ageMix*(para.a.*para.h)');
kage=size(AgeMat); %Check the size
%% SET UP THE DIFFERENTIAL EQUATIONS
for j= 1:soc
for k= 1:soc
SocInf(j,:) = w(j,k)*(para.tau*J(k,:)).*(S(j,:)./para.N(j,:));
size(SocInf)
dpop(1:1:n,:) = - SocInf*AgeMat+para.epsilon * R;
dpop(n+1:1:2*n,:) = SocInf*AgeMat*pi(m).*E - (1-pi(m))*para.theta*E - (1-pi(m))*para.rho*E;
dpop((2*n)+1:1:3*n,:) = (1-pi(m))*para.rho*E - para.gamma*A;
dpop((3*n)+1:1:4*n,:) = (1-pi(m))*para.theta*E - para.delta.*J - para.gamma*J;
dpop((4*n)+1:1:5*n,:) = pi(m)*E - 1/para.PHI*para.p*Q - 1/para.PHI*(1-para.p)*para.gamma*Q;
dpop((5*n)+1:1:6*n,:) = 1/para.PHI*para.p*Q - para.gamma .*I -para.gamma*I;
dpop((6*n)+1:1:7*n,:) = para.gamma .*J+ para.gamma .*I - para.gamma.*H - para.gamma_2*H;
dpop((7*n)+1:1:8*n,:) = para.gamma*J + para.gamma*A +para.gamma_2*H +1/para.PHI*(1-para.p)*para.gamma*Q + para.gamma*I - para.epsilon*R;
dpop((8*n)+1:1:9*n,:)= para.gamma.*H;
dpop((9*n)+1:1:10*n,:) = 1/para.PHI*para.p*Q;
end
end
dpop = dpop(:);
end
end
The initial conditions are set up like this
ICs = struct('S',S,'E',E,'A',zeros(soc,n),'J',zeros(soc,n), ...
'Q',zeros(soc,n),'I',zeros(soc,n),'H',zeros(soc,n),'R',zeros(soc,n), ...
'D',zeros(soc,n),'Ir',zeros(soc,n));
However my output is giving just the re repetition of the first row for the number of times and in each social groups. What am I doing wrong please.
0 Kommentare
Antworten (1)
Torsten
am 21 Sep. 2023
Verschoben: Torsten
am 21 Sep. 2023
You write the same data in the Classes structure for each value of i:
for i = 1:soc
Classes.S(:,:,i) = pop(:,1:1:n);
Classes.E(:,:,i) = pop(:,n+1:1:2*n);
Classes.A(:,:,i) = pop(:,(2*n)+1:1:3*n);
Classes.J(:,:,i) = pop(:,(3*n)+1:1:4*n);
Classes.Q(:,:,i) = pop(:,(4*n)+1:1:5*n);
Classes.I(:,:,i) = pop(:,(5*n)+1:1:6*n);
Classes.H(:,:,i) = pop(:,(6*n)+1:1:7*n);
Classes.R(:,:,i) = pop(:,(7*n)+1:1:8*n);
Classes.D(:,:,i) = pop(:,(8*n)+1:1:9*n);
Classes.Ir(:,:,i) = pop(:,(9*n)+1:1:10*n);
end
0 Kommentare
Siehe auch
Kategorien
Mehr zu Verification, Validation, and Test finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!