Solving ODE with big dimension
조회 수: 1 (최근 30일)
이전 댓글 표시
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
답변 (1개)
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
댓글 수: 0
참고 항목
카테고리
Help Center 및 File Exchange에서 Verification, Validation, and Test에 대해 자세히 알아보기
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!