Error: Unable to perform assignment because the left and right sides have a different number of elements.
6 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
Emman Obuya
am 8 Jan. 2021
Kommentiert: Alan Stevens
am 11 Jan. 2021
Am finding the optimal control for the set of equations but am getting the same error in line 203
Unable to perform assignment because the left and right sides have a different number of elements
lambda_S(j-1) = lambda_S(j)-(h/6)*(k1+2*k2+2*k3+k4);
Have tried to check for similar answers but i seem not to get a way out easily. Your help will be of good benefit and thanks in advance
function y =Try( )
clc;
A1=0.02;
A2=10;
L1=1.6;
L2=4.0;
L3=0.6;
L4=0.8;
Capital_lambda=1000/100;
Alpha=0.5/100;
lambda=0.005/100;
beta=50/100;
k=100000;
b=0.0027/100;
d=1.5/100;
mu_h=0.006/100;
mu_b=2/100;
Gamma=10/100;
r=20/100;
test=-1;
N=10000;
t=linspace(0,100,N+1); % create a matrix starting from 0 going to 100 with 1001 values
h=1/N;
h2=h/2;
U_note=zeros(1,N+1); % Creating Zero matrices
U_one=zeros(1,N+1);
U_two=zeros(1,N+1);
U_three=zeros(1,N+1);
S=zeros(1,N+1);%S
S0=40;
S(1)=S0;
S1=zeros(1,N+1);%S1
S10=40;
S1(1)=S10;
V=zeros(1,N+1);%V
V0=30;
V(1)=V0;
V1=zeros(1,N+1);%V1
V10=30;
V1(1)=V10;
I=zeros(1,N+1);%I
I0=20;
I(1)=I0;
I1=zeros(1,N+1);%I1
I10=20;
I1(1)=I10;
R=zeros(1,N+1);%R
R0=10;
R(1)=R0;
R1=zeros(1,N+1);%R
R10=10;
R1(1)=R10;
B=zeros(1,N+1);%B
B0=20;
B(1)=B0;
B1=zeros(1,N+1);%B
B10=20;
B1(1)=B10;
lambda_S=zeros(1,N+1);
%lambda1(100)=0;
lambda_V=zeros(1,N+1);
%lambda2(100)=0;
lambda_I=zeros(1,N+1);
%lambda3(100)=0;
lambda_R=zeros(1,N+1);
%lambda4(100)=0;
lambda_B=zeros(1,N+1);
%lambda5(100)=0;
while(test<0) %Creates a loop to be repeated when the conditions are true
oldU_note=U_note;
oldU_one=U_one;
oldU_two=U_two;
oldU_three=U_three;
oldS=S;
oldS1=S1;
oldV=V;
oldV1=V1;
oldI=I;
oldI1=I1;
oldR=R;
oldR1=R1;
oldB=B;
oldB1=B1;
oldlambda_S=lambda_S;
oldlambda_V=lambda_V;
oldlambda_I=lambda_I;
oldlambda_R=lambda_R;
oldlambda_B=lambda_B;
for i=1:N
k1=Capital_lambda+Gamma*R(i)+Alpha*V(i)-((1-U_one(i))*Alpha+mu_h+U_note(i))*S(i);
k2=Capital_lambda+Gamma*(R(i)+h2*k1)+Alpha*(V(i)+h2*k1)-((1-0.5*(U_one(i)+U_one(i+1)))*Alpha+mu_h+0.5*(U_note(i)+U_note(i+1))*(S(i)+h2*k1));
k3=Capital_lambda+Gamma*(R(i)+h2*k2)+Alpha*(V(i)+h2*k2)-((1-0.5*(U_one(i)+U_one(i+1)))*Alpha+mu_h+0.5*(U_note(i)+U_note(i+1))*(S(i)+h2*k2));
k4=Capital_lambda+Gamma*(R(i)+h*k3)+Alpha*(V(i)+h*k3)-((1-U_one(i+1))*Alpha+mu_h+U_note(i+1))*(S(i)+h*k3);
S(i+1)=S(i)+(h/6)*(k1+2*k2+2*k3+k4);
end
for i=1:N
k1=Capital_lambda+Gamma*R(i)+Alpha*V(i)-Alpha+mu_h;
k2=Capital_lambda+Gamma*(R(i)+h2*k1)+Alpha*(V(i)+h2*k1)-Alpha+mu_h;
k3=Capital_lambda+Gamma*(R(i)+h2*k2)+Alpha*(V(i)+h2*k2)-Alpha+mu_h;
k4=Capital_lambda+Gamma*(R(i)+h*k3)+Alpha*(V(i)+h*k3)-Alpha+mu_h;
S1(i+1)=S1(i)+(h/6)*(k1+2*k2+2*k3+k4);
end
for i=1:N
k1=U_note(i)*S(i)- mu_h*V(i)-Alpha*V(i);
k2=0.5*(U_note(i)+U_note(i+1))*(S(i)+h2*k1)- mu_h*(V(i)+h2*k1)-Alpha*(V(i)+h2*k1);
k3=0.5*(U_note(i)+U_note(i+1))*(S(i)+h2*k2)- mu_h*(V(i)+h2*k2)-Alpha*(V(i)+h2*k2);
k4=0.5*(U_note(i+1))*(S(i)+h*k3)- mu_h*(V(i)+h*k3)-Alpha*(V(i)+h*k3);
V(i+1)=V(i)+(h/6)*(k1+2*k2+2*k3+k4);
end
for i=1:N
k1=-mu_h*V(i)-Alpha*V(i);
k2=-mu_h*(V(i)+h2*k1)-Alpha*(V(i)+h2*k1);
k3=-mu_h*(V(i)+h2*k2)-Alpha*(V(i)+h2*k2);
k4=-mu_h*(V(i)+h*k3)-Alpha*(V(i)+h*k3);
V1(i+1)=V1(i)+(h/6)*(k1+2*k2+2*k3+k4);
end
for i=1:N
k1=(1-U_one(i))*lambda*S(i)-((1-U_one(i))*beta + mu_h + d+ U_two(i)+r)*I(i);
k2=(1-0.5*(U_one(i)+U_one(i+1)))*lambda*(S(i)+h2*k1)-((1-0.5*(U_one(i)+U_one(i+1)))*beta + mu_h + d+ 0.5*(U_two(i)+U_two(i+1))+r)*(I(i)+h2*k1);
k3=(1-0.5*(U_one(i)+U_one(i+1)))*lambda*(S(i)+h2*k2)-((1-0.5*(U_one(i)+U_one(i+1)))*beta + mu_h + d+ 0.5*(U_two(i)+U_two(i+1))+r)*(I(i)+h2*k2);
k4=(1-U_one(i+1))*lambda*(S(i)+h*k3)-((1-U_one(i+1))*beta + mu_h + d+U_two(i+1))+r*(I(i)+h*k3);
I(i+1)=I(i)+(h/6)*(k1+2*k2+2*k3+k4);
end
for i=1:N
k1=lambda*S(i)-beta + mu_h + d + r*I(i);
k2=lambda*(S(i)+h2*k1)-beta + mu_h + d + r*(I(i)+h2*k1);
k3=lambda*(S(i)+h2*k2)-beta + mu_h + d + r*(I(i)+h2*k2);
k4=lambda*(S(i)+h*k3)-beta + mu_h + d + r*(I(i)+h*k3);
I1(i+1)=I1(i)+(h/6)*(k1+2*k2+2*k3+k4);
end
%fprintf('sum1=%f',sum1);
for i=1:N
k1=(U_two(i)+r)*I(i)-(mu_h+Gamma)*R(i);
k2=(0.5*(U_two(i)+U_two(i+1))+r)*(I(i)+h2*k1)-(mu_h+Gamma)*(R(i)+h2*k1);
k3=(0.5*(U_two(i)+U_two(i+1))+r)*(I(i)+h2*k2)-(mu_h+Gamma)*(R(i)+h2*k2);
k4=(U_two(i+1))+r*(I(i)+h*k3)-(mu_h+Gamma)*(R(i)+h*k3);
R(i+1)=R(i)+(h/6)*(k1+2*k2+2*k3+k4);
end
for i=1:N
k1= r*I(i)-(mu_h+Gamma)*R(i);
k2= r*(I(i)+h2*k1)-(mu_h+Gamma)*(R(i)+h2*k1);
k3= r*(I(i)+h2*k2)-(mu_h+Gamma)*(R(i)+h2*k2);
k4= r*(I(i)+h*k3)-(mu_h+Gamma)*(R(i)+h*k3);
R1(i+1)=R1(i)+(h/6)*(k1+2*k2+2*k3+k4);
end
for i=1:N
k1=beta*(1-U_one(i))*I(i)+((1-B(i)/k)*b*B(i))-(U_three(i)+mu_b)*B(i);
k2=beta*((1-0.5*(U_one(i)+U_one(i+1)))*(I(i)+h2*k1)+((1-(B(i)+h2*k1)/k)*b*(B(i)+h2*k1))-(0.5*(U_three(i)+U_three(i+1))+mu_b)*(B(i)+h2*k1));
k3=beta*((1-0.5*(U_one(i)+U_one(i+1)))*(I(i)+h2*k2)+((1-(B(i)+h2*k2)/k)*b*(B(i)+h2*k2))-(0.5*(U_three(i)+U_three(i+1))+mu_b)*(B(i)+h2*k2));
k4=beta*((1-(U_one(i+1))*(I(i)+h*k3)+((1-(B(i)+h*k3)/k)*b*(B(i)+h*k3))-(U_three(i+1))+mu_b)*(B(i)+h*k3));
B(i+1)=B(i)+(h/6)*(k1+2*k2+2*k3+k4);
end
for i=1:N
k1=beta*I(i)+((1-B(i)/k)*b*B(i))-mu_b*B(i);
k2=beta*(I(i)+h2*k1)+((1-(B(i)+h2*k1)/k)*b*(B(i)+h2*k1))-mu_b*(B(i)+h2*k1);
k3=beta*(I(i)+h2*k2)+((1-(B(i)+h2*k2)/k)*b*(B(i)+h2*k2))-mu_b*(B(i)+h2*k2);
k4=beta*(I(i)+h*k3)+((1-(B(i)+h*k3)/k)*b*(B(i)+h*k3))-mu_b*(B(i)+h*k3);
B1(i+1)=B1(i)+(h/6)*(k1+2*k2+2*k3+k4);
end
for i=1:N
j=N+2-i;
k1=lambda_S(j)*((1-U_one(j))*lambda+mu_h+U_note(j))-A1-lambda_V*U_note(j)-lambda_I*lambda(1-U_one(j));
k2=(lambda_S(j)-h2*k1)*((1-0.5*(U_one(j)+U_one(j-1)))*lambda+mu_h+0.5*(U_note(j)+U_note(j-1)))-A1-(lambda_V-h2*k1)*0.5*(U_note(j)+U_note(j-1))-(lambda_I-h2*k1)*lambda(1-0.5*(U_one(j)+U_one(j-1)));
k3=(lambda_S(j)-h2*k2)*((1-0.5*(U_one(j)+U_one(j-1)))*lambda+mu_h+0.5*(U_note(j)+U_note(j-1)))-A1-(lambda_V-h2*k2)*0.5*(U_note(j)+U_note(j-1))-(lambda_I-h2*k2)*lambda(1-0.5*(U_one(j)+U_one(j-1)));
k4=(lambda_S(j)-h*k3)*((1-U_one(j-1)))*lambda+mu_h+U_note(j-1)-A1-(lambda_V-h*k3)*U_note(j-1)-(lambda_I-h*k3)*lambda(1-U_one(j-1));
lambda_S(j-1) = lambda_S(j)-(h/6)*(k1+2*k2+2*k3+k4);
end
for i=1:N
j=N+2-i;
k1=lambda_V(j)*(mu_h+Alpha)-lambda_S(j)*Alpha;
k2=(lambda_V(j)-h2*k1)*(mu_h+Alpha)-(lambda_S(j)-h2*k1)*Alpha;
k3=(lambda_V(j)-h2*k2)*(mu_h+Alpha)-(lambda_S(j)-h2*k2)*Alpha;
k4=(lambda_V(j)-h*k3)*(mu_h+Alpha)-(lambda_S(j)-h*k3)*Alpha;
lambda_V(j-1)=lambda_V(j)-(h/6)*(k1+2*k2+2*k3+k4);
end
for i=1:N
j=N+2-i;
k1=lambda_I(j)*(beta*(1-U_one(j))+mu_h+d+U_two(j))-A2-lambda_R(j)*Alpha*U_two(j)-lambda_B(j)*beta*(1-U_one(j));
k2=(lambda_I(j)-h2*k1)*(beta*(1-0.5*(U_one(j)+U_one(j-1)))+mu_h+d+(U_two(j)+U_two(j-1)))-A2-(lambda_R(j)-h2*k1)*Alpha*(U_two(j)+U_two(j-1))-(lambda_B(j)-h2*k1)*beta*(1-0.5*(U_one(j)+U_one(j-1)));
k3=(lambda_I(j)-h2*k2)*(beta*(1-0.5*(U_one(j)+U_one(j-1)))+mu_h+d+(U_two(j)+U_two(j-1)))-A2-(lambda_R(j)-h2*k2)*Alpha*(U_two(j)+U_two(j-1))-(lambda_B(j)-h2*k2)*beta*(1-0.5*(U_one(j)+U_one(j-1)));
k4=(lambda_I(j)-h*k3)*(beta*(1-0.5*(U_one(j)+U_one(j-1)))+mu_h+d+(U_two(j)+U_two(j-1)))-A2-(lambda_R(j)-h*k3)*Alpha*(U_two(j)+U_two(j-1))-(lambda_B(j)-h*k3)*beta*(1-0.5*(U_one(j)+U_one(j-1)));
lambda_I(j-1)=lambda_I(j)-(h/6)*(k1+2*k2(j)+2*k3(j)+k4(j));
end
for i=1:N
j=N+2-i;
k1=lambda_R(j)*(mu_h+Gamma)-lambda_S(j)*Gamma;
k2=(lambda_R(j)-h2*k1)*(mu_h+Gamma)-(lambda_S(j)-h2*k1)*Gamma;
k3=(lambda_R(j)-h2*k2)*(mu_h+Gamma)-(lambda_S(j)-h2*k2)*Gamma;
k4=(lambda_R(j)-h*k3)*(mu_h+Gamma)-(lambda_S(j)-h*k3)*Gamma;
lambda_R(j-1)=lambda_R(j)-(h/6)*(k1+2*k2+2*k3+k4);
end
for i=1:N
j=N+2-i;
k1=lambda_B*(mu_b+U_three);
k2=(lambda_B(j)-h2*k1)*(mu_b+0.5*(U_three(j)+U_three(j-1)));
k3=(lambda_B(j)-h2*k2)*(mu_b+0.5*(U_three(j)+U_three(j-1)));
k4=(lambda_B(j)-h*k3)*(mu_b+0.5*(U_three(j)+U_three(j-1)));
lambda_B(j-1)=lambda_B(j)-(h/6)*(k1+2*k2+2*k3+k4);
end
u0=((lambda_S-lambda_V)*S(i))/L1;
U_note=0.5*(u0+oldU_note);
u1=((lambda_I-lambda_S)*lambda*S(i)-lambda_I*beta)/L2;
U_one=0.5*(u1+oldU_one);
u2=((lambda_I-lambda_R)*I(i))/L3;
U_two=0.5*(u2+oldU_two);
u3=(lambda_B*B(i))/L4;
U_three=0.5*(u3+oldU_three);
temp1=0.001*sum(abs(U_note))-sum(abs(oldU_note-U_note));
temp2=0.001*sum(abs(U_one))-sum(abs(oldU_one-U_one));
temp3=0.001*sum(abs(U_two))-sum(abs(oldU_two-U_two));
temp4=0.001*sum(abs(U_three))-sum(abs(oldU_three-U_three));
temp5=0.001*sum(abs(S))-sum(abs(oldS-S));
temp6=0.001*sum(abs(V))-sum(abs(oldV-V));
temp7=0.001*sum(abs(I))-sum(abs(oldI-I));
temp8=0.001*sum(abs(R))-sum(abs(oldR-R));
temp9=0.001*sum(abs(B))-sum(abs(oldB-B));
temp10=0.001*sum(abs(lambda_S))-sum(abs(oldlambda_S-lambda_S));
temp11=0.001*sum(abs(lambda_V))-sum(abs(oldlambda_V-lambda_V));
temp12=0.001*sum(abs(lambda_I))-sum(abs(oldlambda_I-lambda_I));
temp13=0.001*sum(abs(lambda_R))-sum(abs(oldlambda_R-lambda_R));
temp14=0.001*sum(abs(lambda_B))-sum(abs(oldlambda_B-lambda_B));
temp15=0.001*sum(abs(S1))-sum(abs(oldS1-S1));
temp16=0.001*sum(abs(V1))-sum(abs(oldV1-V1));
temp17=0.001*sum(abs(I1))-sum(abs(oldI1-I1));
temp18=0.001*sum(abs(R1))-sum(abs(oldR1-R1));
temp19=0.001*sum(abs(B1))-sum(abs(oldB1-B1));
l1=min(temp2,temp3);
l2=min(temp4,temp5);
l3=min(temp6,temp7);
l4=min(temp8,temp9);
l5=min(temp10,temp11);
l6=min(temp12,temp13);
l7=min(temp14,temp15);
l8=min(temp16,temp17);
l9=min(temp18,temp19);
l11=min(l1,l2);
l12=min(l3,l4);
l13=min(l5,l6);
l14=min(l7,l8);
l15=min(l9,l10);
l16=min(l11,l12);
l17=min(l13,l14);
l18=min(l15,l16);
l19=min(l17,l18);
l20=min(temp1,l19);
test=min(l15,l20);
end
y(:,1)=t;
y(:,2)=S;
y(:,3)=V;
y(:,4)=I;
y(:,5)=R;
y(:,6)=B;
y(:,17)=S1;
y(:,18)=E1;
y(:,19)=I1;
y(:,20)=R1;
y(:,21)=B1;
y(:,14)=U_note;
y(:,15)=U_one;
y(:,16)=U_two;
y(:,17)=U_three;
y(:,8)=lambda_S;
y(:,9)=lambda_V;
y(:,10)=lambda_I;
y(:,11)=lambda_R;
y(:,12)=lambda_B;
plot(t,S);
xlabel('Time');
ylabel('Susceptible human');
plot(t,V);
xlabel('Time')
ylabel('Exposed human');
plot(t,R);
xlabel('Time')
ylabel('Exposed ');
2 Kommentare
Daniel Pollard
am 8 Jan. 2021
This code is extremely long, and as it's a function, I don't think I can easily run it as I don't know how to use it. Can you provide a minimum example which produces the same error?
I tried to make sense of it manually, but each variable is built from multiple others. "Different number of elements" means that lambda_S(j-1) has a different number of components to lambda_S(j)-(h/6)*(k1+2*k2+2*k3+k4);. I suspect the error is in the k's but I can't find it.
Akzeptierte Antwort
Alan Stevens
am 8 Jan. 2021
There are a lot of missing indices and definitions. In the following the indices are repaired; however, you will need to define l10 (see line 256 or thereabouts) to get the routine to work
A1=0.02;
A2=10;
L1=1.6;
L2=4.0;
L3=0.6;
L4=0.8;
Capital_lambda=1000/100;
Alpha=0.5/100;
lambda=0.005/100;
beta=50/100;
k=100000;
b=0.0027/100;
d=1.5/100;
mu_h=0.006/100;
mu_b=2/100;
Gamma=10/100;
r=20/100;
test=-1;
N=1000;
t=linspace(0,100,N+1); % create a matrix starting from 0 going to 100 with 1001 values
h=1/N;
h2=h/2;
U_note=zeros(1,N+1); % Creating Zero matrices
U_one=zeros(1,N+1);
U_two=zeros(1,N+1);
U_three=zeros(1,N+1);
S=zeros(1,N+1);%S
S0=40;
S(1)=S0;
S1=zeros(1,N+1);%S1
S10=40;
S1(1)=S10;
V=zeros(1,N+1);%V
V0=30;
V(1)=V0;
V1=zeros(1,N+1);%V1
V10=30;
V1(1)=V10;
I=zeros(1,N+1);%I
I0=20;
I(1)=I0;
I1=zeros(1,N+1);%I1
I10=20;
I1(1)=I10;
R=zeros(1,N+1);%R
R0=10;
R(1)=R0;
R1=zeros(1,N+1);%R
R10=10;
R1(1)=R10;
B=zeros(1,N+1);%B
B0=20;
B(1)=B0;
B1=zeros(1,N+1);%B
B10=20;
B1(1)=B10;
lambda_S=zeros(1,N+1);
%lambda1(100)=0;
lambda_V=zeros(1,N+1);
%lambda2(100)=0;
lambda_I=zeros(1,N+1);
%lambda3(100)=0;
lambda_R=zeros(1,N+1);
%lambda4(100)=0;
lambda_B=zeros(1,N+1);
%lambda5(100)=0;
while(test<0) %Creates a loop to be repeated when the conditions are true
oldU_note=U_note;
oldU_one=U_one;
oldU_two=U_two;
oldU_three=U_three;
oldS=S;
oldS1=S1;
oldV=V;
oldV1=V1;
oldI=I;
oldI1=I1;
oldR=R;
oldR1=R1;
oldB=B;
oldB1=B1;
oldlambda_S=lambda_S;
oldlambda_V=lambda_V;
oldlambda_I=lambda_I;
oldlambda_R=lambda_R;
oldlambda_B=lambda_B;
for i=1:N
k1=Capital_lambda+Gamma*R(i)+Alpha*V(i)-((1-U_one(i))*Alpha+mu_h+U_note(i))*S(i);
k2=Capital_lambda+Gamma*(R(i)+h2*k1)+Alpha*(V(i)+h2*k1)-((1-0.5*(U_one(i)+U_one(i+1)))*Alpha+mu_h+0.5*(U_note(i)+U_note(i+1))*(S(i)+h2*k1));
k3=Capital_lambda+Gamma*(R(i)+h2*k2)+Alpha*(V(i)+h2*k2)-((1-0.5*(U_one(i)+U_one(i+1)))*Alpha+mu_h+0.5*(U_note(i)+U_note(i+1))*(S(i)+h2*k2));
k4=Capital_lambda+Gamma*(R(i)+h*k3)+Alpha*(V(i)+h*k3)-((1-U_one(i+1))*Alpha+mu_h+U_note(i+1))*(S(i)+h*k3);
S(i+1)=S(i)+(h/6)*(k1+2*k2+2*k3+k4);
end
for i=1:N
k1=Capital_lambda+Gamma*R(i)+Alpha*V(i)-Alpha+mu_h;
k2=Capital_lambda+Gamma*(R(i)+h2*k1)+Alpha*(V(i)+h2*k1)-Alpha+mu_h;
k3=Capital_lambda+Gamma*(R(i)+h2*k2)+Alpha*(V(i)+h2*k2)-Alpha+mu_h;
k4=Capital_lambda+Gamma*(R(i)+h*k3)+Alpha*(V(i)+h*k3)-Alpha+mu_h;
S1(i+1)=S1(i)+(h/6)*(k1+2*k2+2*k3+k4);
end
for i=1:N
k1=U_note(i)*S(i)- mu_h*V(i)-Alpha*V(i);
k2=0.5*(U_note(i)+U_note(i+1))*(S(i)+h2*k1)- mu_h*(V(i)+h2*k1)-Alpha*(V(i)+h2*k1);
k3=0.5*(U_note(i)+U_note(i+1))*(S(i)+h2*k2)- mu_h*(V(i)+h2*k2)-Alpha*(V(i)+h2*k2);
k4=0.5*(U_note(i+1))*(S(i)+h*k3)- mu_h*(V(i)+h*k3)-Alpha*(V(i)+h*k3);
V(i+1)=V(i)+(h/6)*(k1+2*k2+2*k3+k4);
end
for i=1:N
k1=-mu_h*V(i)-Alpha*V(i);
k2=-mu_h*(V(i)+h2*k1)-Alpha*(V(i)+h2*k1);
k3=-mu_h*(V(i)+h2*k2)-Alpha*(V(i)+h2*k2);
k4=-mu_h*(V(i)+h*k3)-Alpha*(V(i)+h*k3);
V1(i+1)=V1(i)+(h/6)*(k1+2*k2+2*k3+k4);
end
for i=1:N
k1=(1-U_one(i))*lambda*S(i)-((1-U_one(i))*beta + mu_h + d+ U_two(i)+r)*I(i);
k2=(1-0.5*(U_one(i)+U_one(i+1)))*lambda*(S(i)+h2*k1)-((1-0.5*(U_one(i)+U_one(i+1)))*beta + mu_h + d+ 0.5*(U_two(i)+U_two(i+1))+r)*(I(i)+h2*k1);
k3=(1-0.5*(U_one(i)+U_one(i+1)))*lambda*(S(i)+h2*k2)-((1-0.5*(U_one(i)+U_one(i+1)))*beta + mu_h + d+ 0.5*(U_two(i)+U_two(i+1))+r)*(I(i)+h2*k2);
k4=(1-U_one(i+1))*lambda*(S(i)+h*k3)-((1-U_one(i+1))*beta + mu_h + d+U_two(i+1))+r*(I(i)+h*k3);
I(i+1)=I(i)+(h/6)*(k1+2*k2+2*k3+k4);
end
for i=1:N
k1=lambda*S(i)-beta + mu_h + d + r*I(i);
k2=lambda*(S(i)+h2*k1)-beta + mu_h + d + r*(I(i)+h2*k1);
k3=lambda*(S(i)+h2*k2)-beta + mu_h + d + r*(I(i)+h2*k2);
k4=lambda*(S(i)+h*k3)-beta + mu_h + d + r*(I(i)+h*k3);
I1(i+1)=I1(i)+(h/6)*(k1+2*k2+2*k3+k4);
end
%fprintf('sum1=%f',sum1);
for i=1:N
k1=(U_two(i)+r)*I(i)-(mu_h+Gamma)*R(i);
k2=(0.5*(U_two(i)+U_two(i+1))+r)*(I(i)+h2*k1)-(mu_h+Gamma)*(R(i)+h2*k1);
k3=(0.5*(U_two(i)+U_two(i+1))+r)*(I(i)+h2*k2)-(mu_h+Gamma)*(R(i)+h2*k2);
k4=(U_two(i+1))+r*(I(i)+h*k3)-(mu_h+Gamma)*(R(i)+h*k3);
R(i+1)=R(i)+(h/6)*(k1+2*k2+2*k3+k4);
end
for i=1:N
k1= r*I(i)-(mu_h+Gamma)*R(i);
k2= r*(I(i)+h2*k1)-(mu_h+Gamma)*(R(i)+h2*k1);
k3= r*(I(i)+h2*k2)-(mu_h+Gamma)*(R(i)+h2*k2);
k4= r*(I(i)+h*k3)-(mu_h+Gamma)*(R(i)+h*k3);
R1(i+1)=R1(i)+(h/6)*(k1+2*k2+2*k3+k4);
end
for i=1:N
k1=beta*(1-U_one(i))*I(i)+((1-B(i)/k)*b*B(i))-(U_three(i)+mu_b)*B(i);
k2=beta*((1-0.5*(U_one(i)+U_one(i+1)))*(I(i)+h2*k1)+((1-(B(i)+h2*k1)/k)*b*(B(i)+h2*k1))-(0.5*(U_three(i)+U_three(i+1))+mu_b)*(B(i)+h2*k1));
k3=beta*((1-0.5*(U_one(i)+U_one(i+1)))*(I(i)+h2*k2)+((1-(B(i)+h2*k2)/k)*b*(B(i)+h2*k2))-(0.5*(U_three(i)+U_three(i+1))+mu_b)*(B(i)+h2*k2));
k4=beta*((1-(U_one(i+1))*(I(i)+h*k3)+((1-(B(i)+h*k3)/k)*b*(B(i)+h*k3))-(U_three(i+1))+mu_b)*(B(i)+h*k3));
B(i+1)=B(i)+(h/6)*(k1+2*k2+2*k3+k4);
end
for i=1:N
k1=beta*I(i)+((1-B(i)/k)*b*B(i))-mu_b*B(i);
k2=beta*(I(i)+h2*k1)+((1-(B(i)+h2*k1)/k)*b*(B(i)+h2*k1))-mu_b*(B(i)+h2*k1);
k3=beta*(I(i)+h2*k2)+((1-(B(i)+h2*k2)/k)*b*(B(i)+h2*k2))-mu_b*(B(i)+h2*k2);
k4=beta*(I(i)+h*k3)+((1-(B(i)+h*k3)/k)*b*(B(i)+h*k3))-mu_b*(B(i)+h*k3);
B1(i+1)=B1(i)+(h/6)*(k1+2*k2+2*k3+k4);
end
for i=1:N
j=N+2-i;
k1=lambda_S(j)*((1-U_one(j))*lambda+mu_h+U_note(j))-A1-lambda_V(j)*U_note(j)-lambda_I(j)*lambda*(1-U_one(j));
k2=(lambda_S(j)-h2*k1)*((1-0.5*(U_one(j)+U_one(j-1)))*lambda+mu_h+0.5*(U_note(j)+U_note(j-1)))-A1-(lambda_V(j)-h2*k1)*0.5*(U_note(j)+U_note(j-1))-(lambda_I(j)-h2*k1)*lambda*(1-0.5*(U_one(j)+U_one(j-1)));
k3=(lambda_S(j)-h2*k2)*((1-0.5*(U_one(j)+U_one(j-1)))*lambda+mu_h+0.5*(U_note(j)+U_note(j-1)))-A1-(lambda_V(j)-h2*k2)*0.5*(U_note(j)+U_note(j-1))-(lambda_I(j)-h2*k2)*lambda*(1-0.5*(U_one(j)+U_one(j-1)));
k4=(lambda_S(j)-h*k3)*((1-U_one(j-1)))*lambda+mu_h+U_note(j-1)-A1-(lambda_V(j)-h*k3)*U_note(j-1)-(lambda_I(j)-h*k3)*lambda*(1-U_one(j-1));
lambda_S(j-1) = lambda_S(j)-(h/6)*(k1+2*k2+2*k3+k4);
end
for i=1:N
j=N+2-i;
k1=lambda_V(j)*(mu_h+Alpha)-lambda_S(j)*Alpha;
k2=(lambda_V(j)-h2*k1)*(mu_h+Alpha)-(lambda_S(j)-h2*k1)*Alpha;
k3=(lambda_V(j)-h2*k2)*(mu_h+Alpha)-(lambda_S(j)-h2*k2)*Alpha;
k4=(lambda_V(j)-h*k3)*(mu_h+Alpha)-(lambda_S(j)-h*k3)*Alpha;
lambda_V(j-1)=lambda_V(j)-(h/6)*(k1+2*k2+2*k3+k4);
end
for i=1:N
j=N+2-i;
k1=lambda_I(j)*(beta*(1-U_one(j))+mu_h+d+U_two(j))-A2-lambda_R(j)*Alpha*U_two(j)-lambda_B(j)*beta*(1-U_one(j));
k2=(lambda_I(j)-h2*k1)*(beta*(1-0.5*(U_one(j)+U_one(j-1)))+mu_h+d+(U_two(j)+U_two(j-1)))-A2-(lambda_R(j)-h2*k1)*Alpha*(U_two(j)+U_two(j-1))-(lambda_B(j)-h2*k1)*beta*(1-0.5*(U_one(j)+U_one(j-1)));
k3=(lambda_I(j)-h2*k2)*(beta*(1-0.5*(U_one(j)+U_one(j-1)))+mu_h+d+(U_two(j)+U_two(j-1)))-A2-(lambda_R(j)-h2*k2)*Alpha*(U_two(j)+U_two(j-1))-(lambda_B(j)-h2*k2)*beta*(1-0.5*(U_one(j)+U_one(j-1)));
k4=(lambda_I(j)-h*k3)*(beta*(1-0.5*(U_one(j)+U_one(j-1)))+mu_h+d+(U_two(j)+U_two(j-1)))-A2-(lambda_R(j)-h*k3)*Alpha*(U_two(j)+U_two(j-1))-(lambda_B(j)-h*k3)*beta*(1-0.5*(U_one(j)+U_one(j-1)));
lambda_I(j-1)=lambda_I(j)-(h/6)*(k1+2*k2+2*k3+k4);
end
for i=1:N
j=N+2-i;
k1=lambda_R(j)*(mu_h+Gamma)-lambda_S(j)*Gamma;
k2=(lambda_R(j)-h2*k1)*(mu_h+Gamma)-(lambda_S(j)-h2*k1)*Gamma;
k3=(lambda_R(j)-h2*k2)*(mu_h+Gamma)-(lambda_S(j)-h2*k2)*Gamma;
k4=(lambda_R(j)-h*k3)*(mu_h+Gamma)-(lambda_S(j)-h*k3)*Gamma;
lambda_R(j-1)=lambda_R(j)-(h/6)*(k1+2*k2+2*k3+k4);
end
for i=1:N
j=N+2-i;
k1=lambda_B(j)*(mu_b+U_three(j));
k2=(lambda_B(j)-h2*k1)*(mu_b+0.5*(U_three(j)+U_three(j-1)));
k3=(lambda_B(j)-h2*k2)*(mu_b+0.5*(U_three(j)+U_three(j-1)));
k4=(lambda_B(j)-h*k3)*(mu_b+0.5*(U_three(j)+U_three(j-1)));
lambda_B(j-1)=lambda_B(j)-(h/6)*(k1+2*k2+2*k3+k4);
end
u0=((lambda_S-lambda_V)*S(i))/L1;
U_note=0.5*(u0+oldU_note);
u1=((lambda_I-lambda_S)*lambda*S(i)-lambda_I*beta)/L2;
U_one=0.5*(u1+oldU_one);
u2=((lambda_I-lambda_R)*I(i))/L3;
U_two=0.5*(u2+oldU_two);
u3=(lambda_B*B(i))/L4;
U_three=0.5*(u3+oldU_three);
temp1=0.001*sum(abs(U_note))-sum(abs(oldU_note-U_note));
temp2=0.001*sum(abs(U_one))-sum(abs(oldU_one-U_one));
temp3=0.001*sum(abs(U_two))-sum(abs(oldU_two-U_two));
temp4=0.001*sum(abs(U_three))-sum(abs(oldU_three-U_three));
temp5=0.001*sum(abs(S))-sum(abs(oldS-S));
temp6=0.001*sum(abs(V))-sum(abs(oldV-V));
temp7=0.001*sum(abs(I))-sum(abs(oldI-I));
temp8=0.001*sum(abs(R))-sum(abs(oldR-R));
temp9=0.001*sum(abs(B))-sum(abs(oldB-B));
temp10=0.001*sum(abs(lambda_S))-sum(abs(oldlambda_S-lambda_S));
temp11=0.001*sum(abs(lambda_V))-sum(abs(oldlambda_V-lambda_V));
temp12=0.001*sum(abs(lambda_I))-sum(abs(oldlambda_I-lambda_I));
temp13=0.001*sum(abs(lambda_R))-sum(abs(oldlambda_R-lambda_R));
temp14=0.001*sum(abs(lambda_B))-sum(abs(oldlambda_B-lambda_B));
temp15=0.001*sum(abs(S1))-sum(abs(oldS1-S1));
temp16=0.001*sum(abs(V1))-sum(abs(oldV1-V1));
temp17=0.001*sum(abs(I1))-sum(abs(oldI1-I1));
temp18=0.001*sum(abs(R1))-sum(abs(oldR1-R1));
temp19=0.001*sum(abs(B1))-sum(abs(oldB1-B1));
l1=min(temp2,temp3);
l2=min(temp4,temp5);
l3=min(temp6,temp7);
l4=min(temp8,temp9);
l5=min(temp10,temp11);
l6=min(temp12,temp13);
l7=min(temp14,temp15);
l8=min(temp16,temp17);
l9=min(temp18,temp19);
l11=min(l1,l2);
l12=min(l3,l4);
l13=min(l5,l6);
l14=min(l7,l8);
l15=min(l9,l10); % l10 not defined
l16=min(l11,l12);
l17=min(l13,l14);
l18=min(l15,l16); %not defined
l19=min(l17,l18); % "
l20=min(temp1,l19);% "
test=min(l15,l20); % "
end
y(:,1)=t;
y(:,2)=S;
y(:,3)=V;
y(:,4)=I;
y(:,5)=R;
y(:,6)=B;
y(:,17)=S1;
y(:,18)=E1;
y(:,19)=I1;
y(:,20)=R1;
y(:,21)=B1;
y(:,14)=U_note;
y(:,15)=U_one;
y(:,16)=U_two;
y(:,17)=U_three;
y(:,8)=lambda_S;
y(:,9)=lambda_V;
y(:,10)=lambda_I;
y(:,11)=lambda_R;
y(:,12)=lambda_B;
plot(t,S);
xlabel('Time');
ylabel('Susceptible human');
plot(t,V);
xlabel('Time')
ylabel('Exposed human');
plot(t,R);
xlabel('Time')
ylabel('Exposed ');
4 Kommentare
Alan Stevens
am 8 Jan. 2021
Bearbeitet: Alan Stevens
am 8 Jan. 2021
Starting Immediately above the ‘while’ command, put
its = 0;
Then, change the ‘while’ command to, say,
while (test<0) && (its<100)
its = its + 1;
% ... rest of the code in the loop
end
However, you should still check your equations to see why test always seems to be negative.
Weitere Antworten (5)
Emman Obuya
am 10 Jan. 2021
1 Kommentar
Alan Stevens
am 10 Jan. 2021
To plot say U_note and U_one on the same graph, do something like
plot(t,U_note,t,U_one), grid
xlabel('time'), ylabel('Unote and Uone')
legend('U_note','U_one')
Emman Obuya
am 10 Jan. 2021
1 Kommentar
Alan Stevens
am 10 Jan. 2021
I'm not clear what you mean! You can plot other lines on the same graph by extending as follows:
plot(t,S,t,S1,t,U_note,t,U_one,t,U_three)
and then add labels in legend.
Try
doc plot
for more detail about using the plot command.
Emman Obuya
am 11 Jan. 2021
1 Kommentar
Alan Stevens
am 11 Jan. 2021
I'm still not clear what you are trying to plot! The commands
plot(t,S,t,V,t,I,t,R,t,B),grid
xlabel('t'),ylabel('SVIRB')
legend('S','V','I','R','B')
produce the following:
That's after 100 iterations, but I can't tell if it's what you are looking for.
Emman Obuya
am 11 Jan. 2021
1 Kommentar
Alan Stevens
am 11 Jan. 2021
Well, these plot commands at the end of your program will produce those style lines:
figure
plot(t,S,'b',t,S1,'r'),grid
xlabel('Time');
ylabel('Susceptible');
legend('S','S1')
figure
plot(t,I,'b',t,I1,'r'),grid
xlabel('Time')
ylabel('Infected');
legend('I','I1')
However, they are a different shape from the ones you show above.
Siehe auch
Kategorien
Mehr zu Multiobjective Optimization 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!