how to make a for loop in an ODE 45

1 次查看（过去 30 天）
Zhukun Wang
Zhukun Wang 2021-3-18
回答： Alan Stevens 2021-3-18
data = [0 86; 1 86; 2 117; 6 120; 7 130; 8 135; 13 169; 16 179; 23 224; 27 230; 29 242; 36 234; 41 244; 50 245; 59 270; 63 271; 64 309; 69 354; 72 438; 77 476; 78 508; 85 528; 91 599; 99 759; 104 779; 105 844; 111 888; 113 964; 118 1048; 121 1093; 125 1201; 128 1322; 131 1437; 132 1617; 136 1766; 140 1835; 141 1963; 143 2115; 147 2225; 149 2458; 150 2599; 156 3052; 165 3944; 167 4269; 171 4366; 175 4963; 177 5325; 181 5843; 183 6242; 185 6553; 190 7157; 192 7470; 197 8011; 199 8376; 204 8973; 206 9191; 211 9911;214 10114; 218 13676; 220 13540; 225 13015; 227 13241; 232 14068; 234 14383; 239 15113; 241 15319; 246 15901; 248 16899; 253 17111; 260 17908; 267 18569; 274 19463; 281 20171; 288 20712; 295 21261; 302 21689; 309 22057; 316 22460; 323 22859; 330 23218; 337 23694; 344 23934; 351 24247; 358 24666; 365 24872; 372 25178; 379 25515; 386 25791; 393 26044; 400 26277; 407 26593; 414 26724; 421 26933; 428 27013; 435 27145; 442 27237; 449 27305; 456 27443; 463 27514; 470 27573; 477 27642; 484 27705; 491 27748; 498 27862; 505 27929; 512 27952; 519 28005; 527 28073; 534 28147; 541 28220; 548 28295; 555 28388; 562 28421; 569 28454; 576 28476; 583 28539; 590 28571; 596 28599; 603 28598; 610 28601; 617 28601; 624 28601; 631 28604; 638 28601; 645 28601; 652 28601; 659 28601];
function Math462hw3Q5parta
alpha=0.05;
beta=0.2;
gamma=0.125;
N=50000;
S=1-I-E;
for i= 0:1:659
I=data(i,2)/N;
end
E=2*I;
R=0;
y=N*alpha*E;
vector=[S,E,I,R,y];
tspan=0:1:659;
[t,soln]=ode45(@(t,vector)dotfunctions(t,vector,alpha,beta,gamma),tspan,vector);
%plot function
plot(t,soln);
xlabel('time');
ylabel('value');
title('SIERy model');
function ddt=dotfunctions(~,siery,alpha,beta,gamma)
S=siery(1);
I=siery(2);
R=siery(3);
E=siery(4);
y=siery(5);
Sdot=-beta*S*I;
Edot=beta*S*I-alpha*E;
Idot=alpha*E-gamma*I;
Rdot=gamma*I;
ydot=N*alpha*E;
ddt=[Sdot;Idot;Rdot;Edot;ydot];
end
end
This is my code for simulating a model using ODE45, here are two questions:
  1. Can I write a for loop in the first few lines just to loop through all values of i
  2. Just wondering can someone help me figure out why this code is not running.

请先登录，再回答此问题。

回答（1 个）

Alan Stevens
Alan Stevens 2021-3-18
More like this I think:
data = [0 86; 1 86; 2 117; 6 120; 7 130; 8 135; 13 169; 16 179; 23 224; ...
27 230; 29 242; 36 234; 41 244; 50 245; 59 270; 63 271; 64 309; ...
69 354; 72 438; 77 476; 78 508; 85 528; 91 599; 99 759; 104 779; ...
105 844; 111 888; 113 964; 118 1048; 121 1093; 125 1201; 128 1322; ...
131 1437; 132 1617; 136 1766; 140 1835; 141 1963; 143 2115; 147 2225; ...
149 2458; 150 2599; 156 3052; 165 3944; 167 4269; 171 4366; 175 4963; ...
177 5325; 181 5843; 183 6242; 185 6553; 190 7157; 192 7470; 197 8011; ...
199 8376; 204 8973; 206 9191; 211 9911;214 10114; 218 13676; 220 13540;...
225 13015; 227 13241; 232 14068; 234 14383; 239 15113; 241 15319; ...
246 15901; 248 16899; 253 17111; 260 17908; 267 18569; 274 19463; ...
281 20171; 288 20712; 295 21261; 302 21689; 309 22057; 316 22460; ...
323 22859; 330 23218; 337 23694; 344 23934; 351 24247; 358 24666; ...
365 24872; 372 25178; 379 25515; 386 25791; 393 26044; 400 26277; ...
407 26593; 414 26724; 421 26933; 428 27013; 435 27145; 442 27237; ...
449 27305; 456 27443; 463 27514; 470 27573; 477 27642; 484 27705; ...
491 27748; 498 27862; 505 27929; 512 27952; 519 28005; 527 28073; ...
534 28147; 541 28220; 548 28295; 555 28388; 562 28421; 569 28454; ...
576 28476; 583 28539; 590 28571; 596 28599; 603 28598; 610 28601; ...
617 28601; 624 28601; 631 28604; 638 28601; 645 28601; 652 28601; 659 28601];
[nrows, ncols] = size(data);
alpha=0.05;
beta=0.2;
gamma=0.125;
N=50000;
for i= 1:nrows % indices start with 1 in Matlab, not 0
I=data(i,2)/N;
end
E=2*I;
S=1-I-E; % needs to come after I and E are defined
R=0;
y=N*alpha*E;
vector=[S,E,I,R,y];
tspan=0:1:659;
[t,soln]=ode45(@(t,vector)dotfunctions(t,vector,alpha,beta,gamma,N),tspan,vector);
%plot function
subplot(2,1,1) % separate plots as y scale swamps the others
plot(t,soln(:,1:4)), grid
xlabel('time');
ylabel('value');
legend('S','E','I','R')
title('SIERy model');
subplot(2,1,2)
plot(t,soln(:,5)), grid
xlabel('time');
ylabel('y value');
function ddt=dotfunctions(~,siery,alpha,beta,gamma,N) % need to pass in N as well
S=siery(1);
I=siery(2);
% R=siery(3);
E=siery(4);
% y=siery(5);
Sdot=-beta*S*I;
Edot=beta*S*I-alpha*E;
Idot=alpha*E-gamma*I;
Rdot=gamma*I;
ydot=N*alpha*E;
ddt=[Sdot;Idot;Rdot;Edot;ydot];
end

类别

在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by