Solving ODE with big dimension

1 次查看(过去 30 天)
Gbeminiyi
Gbeminiyi 2023-9-21
移动Torsten 2023-9-21
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.

回答(1 个)

Torsten
Torsten 2023-9-21
移动:Torsten 2023-9-21
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

类别

Help CenterFile Exchange 中查找有关 Verification, Validation, and Test 的更多信息

产品


版本

R2023b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by