Error this answer is NaN, in my runge-kutta code

3 views (last 30 days)
Muhammad Giovani
Muhammad Giovani on 27 Nov 2021
Answered: the cyclist on 27 Nov 2021
hello guys
I'm compiling a code in research, code like a
close
clear
clc
C1 = 1;
C2 = 1;
A1 = 1;
A2 = 1;
S10 = 1000;
S20 = 100;
I10 = 100;
I20 = 60;
R0 = 30;
V0 = 35;
K0 = 20;
a = 0.8;
alpha = 0.5;
miu1 = 0.5;
miu2 = 0.15;
gamma = 0.0111;
n = 0.2;
delta = 0.01;
tau = 0.2;
test = -1;
delta1 = 0.1;
T = 12;
N = 1000;
t = linspace(0,T,N+1);
h = T/N;
h2 = h/2;
u1 = zeros(1,N+1);
u2 = zeros(1,N+1);
S1 = zeros(1,N+1); S1(1) = S10;
S2 = zeros(1,N+1); S2(1) = S20;
I1 = zeros(1,N+1); I1(1) = I10;
I2 = zeros(1,N+1); I2(1) = I20;
V = zeros(1,N+1); V(1) = V0;
R = zeros(1,N+1); R(1) = R0;
K = zeros(1,N+1); K(1) = K0;
Lambda1 = zeros(1,N+1); Lambda1(N+1) = 0;
Lambda2 = zeros(1,N+1); Lambda2(N+1) = 0;
Lambda3 = zeros(1,N+1); Lambda3(N+1) = 0;
Lambda4 = zeros(1,N+1); Lambda4(N+1) = 0;
Lambda5 = zeros(1,N+1); Lambda5(N+1) = 0;
Lambda6 = zeros(1,N+1); Lambda6(N+1) = 0;
Lambda7 = zeros(1,N+1); Lambda7(N+1) = 0;
S11 = zeros(1,N+1); S11(1) = S10;
S21 = zeros(1,N+1); S21(1) = S20;
I11 = zeros(1,N+1); I11(1) = I10;
I21 = zeros(1,N+1); I21(1) = I20;
V1 = zeros(1,N+1); V1(1) = V0;
R1 = zeros(1,N+1); R1(1) = R0;
K1 = zeros(1,N+1); K1(1) = K0;
while(test < 0)
oldu1 = u1;
oldu2 = u2;
oldS1 = S1;
oldS2 = S2;
oldI1 = I1;
oldI2 = I2;
oldV = V;
oldR = R;
oldK = K;
oldLambda1 = Lambda1;
oldLambda2 = Lambda2;
oldLambda3 = Lambda3;
oldLambda4 = Lambda4;
oldLambda5 = Lambda5;
oldLambda6 = Lambda6;
oldLambda7 = Lambda7;
for j = 1:N
k11 = a - n*S1(j) - u1(j)*S1(j) - alpha*S1(j)*(I1(j)+I2(j)) - miu1*S1(j);
k12 = n*S1(j) - alpha*S2(j)*(I2(j)+I1(j)) - miu1*S2(j);
k13 = alpha*S1(j)*(I1(j)+I2(j)) + gamma*V(j)*(I1(j)+I2(j)) - (delta*I1(j)+u2(j)*K(j)) - (miu1+miu2)*I1(j);
k14 = alpha*S2(j)*(I2(j)+I1(j)) - (delta*I2(j)+u2(j)*K(j)) - (miu1+miu2)*I2(j);
k15 = u1(j)*S1(j) - gamma*V(j)*(I1(j)+I2(j)) - miu1*V(j);
k16 = (delta*I1(j)+u2(j)*K(j)) + (delta*I2(j)+u2(j)*K(j)) - tau*R(j) - miu1*R(j);
k17 = tau*R(j) - miu1*K(j);
k21 = a - n*(S1(j)+h2*k11) - (0.5*(u1(j)+u1(j+1)))*(S1(j)+h2*k11) - alpha*(S1(j)+h2*k11)*((I1(j)+h2*k11)+(I2(j)+h2*k11)) - miu1*(S1(j)+h2*k11);
k22 = n*(S1(j)+h2*k12) - alpha*(S2(j)+h2*k12)*((I2(j)+h2*k12)+(I1(j)+h2*k12)) - miu1*(S2(j)+h2*k12);
k23 = alpha*(S1(j)+h2*k13)*((I1(j)+h2*k13)+(I2(j)+h2*k13)) + gamma*(V(j)+h2*k13)*((I1(j)+h2*k13)+(I2(j)+h2*k13)) - (delta*(I1(j)+h2*k13)+(0.5*(u2(j)+u2(j+1)))*(K(j)+h2*k13)) - (miu1+miu2)*(I1(j)+h2*k13);
k24 = alpha*(S2(j)+h2*k14)*((I2(j)+h2*k14)+(I1(j)+h2*k14)) - (delta*(I2(j)+h2*k14)+(0.5*(u2(j)+u2(j+1)))*(K(j)+h2*k14)) - (miu1+miu2)*(I2(j)+h2*k14);
k25 = (0.5*(u1(j)+u1(j+1)))*(S1(j)+h2*k15) - gamma*(V(j)+h2*k15)*((I1(j)+h2*k15)+(I2(j)+h2*k15)) - miu1*(V(j)+h2*k15);
k26 = (delta*(I1(j)+h2*k16)+(0.5*(u2(j)+u2(j+1)))*(K(j)+h2*k16)) + (delta*(I2(j)+h2*k16)+(0.5*(u2(j)+u2(j+1)))*(K(j)+h2*k16)) - tau*(R(j)+h2*k16) - miu1*(R(j)+h2*k16);
k27 = tau*(R(j)+h2*k17) - miu1*(K(j)+h2*k17);
k31 = a - n*(S1(j)+h2*k21) - (0.5*(u1(j)+u1(j+1)))*(S1(j)+h2*k21) - alpha*(S1(j)+h2*k21)*((I1(j)+h2*k21)+(I2(j)+h2*k21)) - miu1*(S1(j)+h2*k21);
k32 = n*(S1(j)+h2*k22) - alpha*(S2(j)+h2*k22)*((I2(j)+h2*k22)+(I1(j)+h2*k22)) - miu1*(S2(j)+h2*k22);
k33 = alpha*(S1(j)+h2*k23)*((I1(j)+h2*k23)+(I2(j)+h2*k23)) + gamma*(V(j)+h2*k23)*((I1(j)+h2*k23)+(I2(j)+h2*k23)) - (delta*(I1(j)+h2*k23)+(0.5*(u2(j)+u2(j+1)))*(K(j)+h2*k23)) - (miu1+miu2)*(I1(j)+h2*k23);
k34 = alpha*(S2(j)+h2*k24)*((I2(j)+h2*k24)+(I1(j)+h2*k24)) - (delta*(I2(j)+h2*k24)+(0.5*(u2(j)+u2(j+1)))*(K(j)+h2*k24)) - (miu1+miu2)*(I2(j)+h2*k24);
k35 = (0.5*(u1(j)+u1(j+1)))*(S1(j)+h2*k25) - gamma*(V(j)+h2*k25)*((I1(j)+h2*k25)+(I2(j)+h2*k25)) - miu1*(V(j)+h2*k25);
k36 = (delta*(I1(j)+h2*k26)+(0.5*(u2(j)+u2(j+1)))*(K(j)+h2*k26)) + (delta*(I2(j)+h2*k26)+(0.5*(u2(j)+u2(j+1)))*(K(j)+h2*k26)) - tau*(R(j)+h2*k26) - miu1*(R(j)+h2*k26);
k37 = tau*(R(j)+h2*k27) - miu1*(K(j)+h2*k27);
k41 = a - n*(S1(j)+h2*k31) - (0.5*(u1(j)+u1(j+1)))*(S1(j)+h2*k31) - alpha*(S1(j)+h2*k31)*((I1(j)+h2*k31)+(I2(j)+h2*k31)) - miu1*(S1(j)+h2*k31);
k42 = n*(S1(j)+h2*k32) - alpha*(S2(j)+h2*k32)*((I2(j)+h2*k32)+(I1(j)+h2*k32)) - miu1*(S2(j)+h2*k32);
k43 = alpha*(S1(j)+h2*k33)*((I1(j)+h2*k33)+(I2(j)+h2*k33)) + gamma*(V(j)+h2*k33)*((I1(j)+h2*k33)+(I2(j)+h2*k33)) - (delta*(I1(j)+h2*k33)+(0.5*(u2(j)+u2(j+1)))*(K(j)+h2*k33)) - (miu1+miu2)*(I1(j)+h2*k33);
k44 = alpha*(S2(j)+h2*k34)*((I2(j)+h2*k34)+(I1(j)+h2*k34)) - (delta*(I2(j)+h2*k34)+(0.5*(u2(j)+u2(j+1)))*(K(j)+h2*k34)) - (miu1+miu2)*(I2(j)+h2*k34);
k45 = (0.5*(u1(j)+u1(j+1)))*(S1(j)+h2*k35) - gamma*(V(j)+h2*k35)*((I1(j)+h2*k35)+(I2(j)+h2*k35)) - miu1*(V(j)+h2*k35);
k46 = (delta*(I1(j)+h2*k36)+(0.5*(u2(j)+u2(j+1)))*(K(j)+h2*k36)) + (delta*(I2(j)+h2*k36)+(0.5*(u2(j)+u2(j+1)))*(K(j)+h2*k36)) - tau*(R(j)+h2*k36) - miu1*(R(j)+h2*k36);
k47 = tau*(R(j)+h2*k37) - miu1*(K(j)+h2*k37);
S1(j+1) = S1(j) + (h/6)*(k11 + 2*k21 + 2*k31 + k41);
S2(j+1) = S2(j) + (h/6)*(k12 + 2*k22 + 2*k32 + k42);
I1(j+1) = I1(j) + (h/6)*(k13 + 2*k23 + 2*k33 + k43);
I2(j+1) = I2(j) + (h/6)*(k14 + 2*k24 + 2*k34 + k44);
V(j+1) = V(j) + (h/6)*(k15 + 2*k25 + 2*k35 + k45);
R(j+1) = R(j) + (h/6)*(k16 + 2*k26 + 2*k36 + k46);
K(j+1) = K(j) + (h/6)*(k17 + 2*k27 + 2*k37 + k47);
end
for j = 1:N
i = N+2-j;
k51 = Lambda1(i)*n + Lambda1(i)*u1(i) + Lambda1(i)*alpha*(I1(i)+I2(i)) + Lambda1(1)*miu1 - Lambda2(i)*n - Lambda3(i)*(alpha*(I1(i)+I2(i))) - Lambda5(i)*u1(i);
k52 = Lambda2(i)*(alpha*((I2(i)+I1(i)))+miu1) - Lambda4(i)*(alpha*(I2(i)+I1(i)));
k53 = - A1*I1(i) + Lambda1(i)*(alpha*S1(i)) + Lambda2(i)*(alpha*S2(i)) - Lambda3(i)*((alpha*S1(i))+(gamma*V(i))-delta-(miu1+miu2)) - Lambda4(i)*((alpha*S2(i))) + Lambda5(i)*(gamma*V(i)) - Lambda6(i)*(delta);
k54 = - A2*I2(i) + Lambda1(i)*(alpha*S1(i)) + Lambda2(i)*(alpha*S2(i)) - Lambda3(i)*((alpha*S1(i))+(gamma*V(i))) - Lambda4(i)*((alpha*S2(i)) - delta - (miu1+miu2)) + Lambda5(i)*(gamma*V(i)) - Lambda6(i)*delta;
k55 = - Lambda3(i)*(gamma*(I1(i)+I2(i))) + Lambda5(i)*(gamma*((I1(i)+I2(i))+miu1));
k56 = Lambda6(i)*(tau+miu1) - Lambda7(i)*(tau);
k57 = Lambda3(i)*u2(i) + Lambda4(i)*u2(i) - Lambda6(i)*2*u2(i) + Lambda7(i)*miu1;
k61 = (Lambda1(i)-h2*k51)*n + (Lambda1(i)-h2*k51)*(0.5*(u1(i)+u1(i-1))) + (Lambda1(i)-h2*k51)*alpha*((I1(i)-h2*k51)+(I2(i)-h2*k51)) + (Lambda1(i)-h2*k51)*miu1 - (Lambda2(i)-h2*k51)*n - (Lambda3(i)-h2*k51)*(alpha*((I1(i)-h2*k51)+(I2(i)-h2*k51))) - (Lambda5(i)-h2*k51)*(0.5*(u1(i)+u1(i-1)));
k62 = (Lambda2(i)-h2*k52)*(alpha*(((I2(i)-h2*k52)+(I1(i)-h2*k52)))+miu1) - (Lambda4(i)-h2*k52)*(alpha*((I2(i)-h2*k52)+(I1(i)-h2*k52)));
k63 = - A1*(I1(i)-h2*k53) + (Lambda1(i)-h2*k53)*(alpha*(S1(i)-h2*k53)) + (Lambda2(i)-h2*k53)*(alpha*(S2(i)-h2*k53)) - (Lambda3(i)-h2*k53)*((alpha*(S1(i)-h2*k53))+(gamma*(V(i)-h2*k53))-delta-(miu1+miu2)) - (Lambda4(i)-h2*k53)*((alpha*(S2(i)-h2*k53))) + (Lambda5(i)-h2*k53)*(gamma*(V(i)-h2*k53)) - (Lambda6(i)-h2*k53)*(delta);
k64 = - A2*(I2(i)-h2*k54) + (Lambda1(i)-h2*k54)*(alpha*(S1(i)-h2*k54)) + (Lambda2(i)-h2*k54)*(alpha*(S2(i)-h2*k54)) - (Lambda3(i)-h2*k54)*((alpha*(S1(i)-h2*k54))+(gamma*(V(i)-h2*k54))) - (Lambda4(i)-h2*k54)*((alpha*(S2(i)-h2*k54)) - delta - (miu1+miu2)) + (Lambda5(i)-h2*k54)*(gamma*(V(i)-h2*k54)) - (Lambda6(i)-h2*k54)*delta;
k65 = - (Lambda3(i)-h2*k55)*(gamma*((I1(i)-h2*k55)+(I2(i)-h2*k55))) + (Lambda5(i)-h2*k55)*(gamma*(((I1(i)-h2*k55)+(I2(i)-h2*k55))+miu1));
k66 = (Lambda6(i)-h2*k56)*(tau+miu1) - (Lambda7(i)-h2*k56)*(tau);
k67 = (Lambda3(i)-h2*k57)*(0.5*(u2(i)+u2(i-1))) + (Lambda4(i)-h2*k57)*(0.5*(u2(i)+u2(i-1))) - (Lambda6(i)-h2*k57)*2*(0.5*(u2(i)+u2(i-1))) + (Lambda7(i)-h2*k57)*miu1;
k71 = (Lambda1(i)-h2*k61)*n + (Lambda1(i)-h2*k61)*(0.5*(u1(i)+u1(i-1))) + (Lambda1(i)-h2*k61)*alpha*((I1(i)-h2*k61)+(I2(i)-h2*k61)) - (Lambda1(i)-h2*k61)*miu1 - (Lambda2(i)-h2*k61)*n - (Lambda3(i)-h2*k61)*(alpha*((I1(i)-h2*k61)+(I2(i)-h2*k61))) - (Lambda5(i)-h2*k61)*(0.5*(u1(i)+u1(i-1)));
k72 = (Lambda2(i)-h2*k62)*(alpha*(((I2(i)-h2*k62)+(I1(i)-h2*k62)))+miu1) - (Lambda4(i)-h2*k62)*(alpha*((I2(i)-h2*k62)+(I1(i)-h2*k62)));
k73 = - A1*(I1(i)-h2*k63) + (Lambda1(i)-h2*k63)*(alpha*(S1(i)-h2*k63)) + (Lambda2(i)-h2*k63)*(alpha*(S2(i)-h2*k63)) - (Lambda3(i)-h2*k63)*((alpha*(S1(i)-h2*k63))+(gamma*(V(i)-h2*k63))-delta-(miu1+miu2)) - (Lambda4(i)-h2*k63)*((alpha*(S2(i)-h2*k63))) + (Lambda5(i)-h2*k63)*(gamma*(V(i)-h2*k63)) - (Lambda6(i)-h2*k63)*(delta);
k74 = - A2*(I2(i)-h2*k64) + (Lambda1(i)-h2*k64)*(alpha*(S1(i)-h2*k64)) + (Lambda2(i)-h2*k64)*(alpha*(S2(i)-h2*k64)) - (Lambda3(i)-h2*k64)*((alpha*(S1(i)-h2*k64))+(gamma*(V(i)-h2*k64))) - (Lambda4(i)-h2*k64)*((alpha*(S2(i)-h2*k64)) - delta - (miu1+miu2)) + (Lambda5(i)-h2*k64)*(gamma*(V(i)-h2*k64)) - (Lambda6(i)-h2*k64)*delta;
k75 = - (Lambda3(i)-h2*k65)*(gamma*((I1(i)-h2*k65)+(I2(i)-h2*k65))) + (Lambda5(i)-h2*k65)*(gamma*(((I1(i)-h2*k65)+(I2(i)-h2*k65))+miu1));
k76 = (Lambda6(i)-h2*k66)*(tau+miu1) - (Lambda7(i)-h2*k66)*(tau);
k77 = (Lambda3(i)-h2*k67)*(0.5*(u2(i)+u2(i-1))) + (Lambda4(i)-h2*k67)*(0.5*(u2(i)+u2(i-1))) - (Lambda6(i)-h2*k67)*2*(0.5*(u2(i)+u2(i-1))) + (Lambda7(i)-h2*k67)*miu1;
k81 = (Lambda1(i)-h2*k71)*n + (Lambda1(i)-h2*k71)*(0.5*(u1(i)+u1(i-1))) + (Lambda1(i)-h2*k71)*alpha*((I1(i)-h2*k71)+(I2(i)-h2*k71)) - (Lambda1(i)-h2*k71)*miu1 - (Lambda2(i)-h2*k71)*n - (Lambda3(i)-h2*k71)*(alpha*((I1(i)-h2*k71)+(I2(i)-h2*k71))) - (Lambda5(i)-h2*k71)*(0.5*(u1(i)+u1(i-1)));
k82 = (Lambda2(i)-h2*k72)*(alpha*(((I2(i)-h2*k72)+(I1(i)-h2*k72)))+miu1) - (Lambda4(i)-h2*k72)*(alpha*((I2(i)-h2*k72)+(I1(i)-h2*k72)));
k83 = - A1*(I1(i)-h2*k73) + (Lambda1(i)-h2*k73)*(alpha*(S1(i)-h2*k73)) + (Lambda2(i)-h2*k73)*(alpha*(S2(i)-h2*k73)) - (Lambda3(i)-h2*k73)*((alpha*(S1(i)-h2*k73))+(gamma*(V(i)-h2*k73))-delta-(miu1+miu2)) - (Lambda4(i)-h2*k73)*((alpha*(S2(i)-h2*k73))) + (Lambda5(i)-h2*k73)*(gamma*(V(i)-h2*k73)) - (Lambda6(i)-h2*k73)*(delta);
k84 = - A2*(I2(i)-h2*k74) + (Lambda1(i)-h2*k74)*(alpha*(S1(i)-h2*k74)) + (Lambda2(i)-h2*k74)*(alpha*(S2(i)-h2*k74)) - (Lambda3(i)-h2*k74)*((alpha*(S1(i)-h2*k74))+(gamma*(V(i)-h2*k74))) - (Lambda4(i)-h2*k74)*((alpha*(S2(i)-h2*k74)) - delta - (miu1+miu2)) + (Lambda5(i)-h2*k74)*(gamma*(V(i)-h2*k74)) - (Lambda6(i)-h2*k74)*delta;
k85 = - (Lambda3(i)-h2*k75)*(gamma*((I1(i)-h2*k75)+(I2(i)-h2*k75))) + (Lambda5(i)-h2*k75)*(gamma*(((I1(i)-h2*k75)+(I2(i)-h2*k75))+miu1));
k86 = (Lambda6(i)-h2*k76)*(tau+miu1) - (Lambda7(i)-h2*k76)*(tau);
k87 = (Lambda3(i)-h2*k77)*(0.5*(u2(i)+u2(i-1))) + (Lambda4(i)-h2*k77)*(0.5*(u2(i)+u2(i-1))) - (Lambda6(i)-h2*k77)*2*(0.5*(u2(i)+u2(i-1))) + (Lambda7(i)-h2*k77)*miu1;
Lambda1(i-1) = Lambda1(i) - (h/6)*(k51 + 2*k61 + 2*k71 + k81);
Lambda2(i-1) = Lambda2(i) - (h/6)*(k52 + 2*k62 + 2*k72 + k82);
Lambda3(i-1) = Lambda3(i) - (h/6)*(k53 + 2*k63 + 2*k73 + k83);
Lambda4(i-1) = Lambda4(i) - (h/6)*(k54 + 2*k64 + 2*k74 + k84);
Lambda5(i-1) = Lambda5(i) - (h/6)*(k55 + 2*k65 + 2*k75 + k85);
Lambda6(i-1) = Lambda6(i) - (h/6)*(k56 + 2*k66 + 2*k76 + k86);
Lambda7(i-1) = Lambda7(i) - (h/6)*(k57 + 2*k67 + 2*k77 + k87);
end
for k = 1:N
k111 = a - n*S1(k) - alpha*S1(k)*(I1(k)+I2(k)) - miu1*S1(k);
k112 = n*S1(k) - alpha*S2(k)*(I2(k)+I1(k)) - miu1*S2(k);
k113 = alpha*S1(k)*(I1(k)+I2(k)) + gamma*V(k)*(I1(k)+I2(k)) - (delta*I1(k)) - (miu1+miu2)*I1(k);
k114 = alpha*S2(k)*(I2(k)+I1(k)) - (delta*I2(k)) - (miu1+miu2)*I2(k);
k115 = - gamma*V(k)*(I1(k)+I2(k)) - miu1*V(k);
k116 = (delta*I1(k)) + (delta*I2(k)) - tau*R(k) - miu1*R(k);
k117 = tau*R(k) - miu1*K(k);
k121 = a - n*(S1(k)+h2*k111) - alpha*(S1(k)+h2*k111)*((I1(k)+h2*k111)+(I2(k)+h2*k111)) - miu1*(S1(k)+h2*k111);
k122 = n*(S1(k)+h2*k112) - alpha*(S2(k)+h2*k112)*((I2(k)+h2*k112)+(I1(k)+h2*k112)) - miu1*(S2(k)+h2*k112);
k123 = alpha*(S1(k)+h2*k113)*((I1(k)+h2*k113)+(I2(k)+h2*k113)) + gamma*(V(k)+h2*k113)*((I1(k)+h2*k113)+(I2(k)+h2*k113)) - (delta*(I1(k)+h2*k113)) - (miu1+miu2)*(I1(k)+h2*k113);
k124 = alpha*(S2(k)+h2*k114)*((I2(k)+h2*k114)+(I1(k)+h2*k114)) - (delta*(I2(k)+h2*k114)) - (miu1+miu2)*(I2(k)+h2*k114);
k125 = - gamma*(V(k)+h2*k115)*((I1(k)+h2*k115)+(I2(k)+h2*k115)) - miu1*(V(k)+h2*k115);
k126 = (delta*(I1(k)+h2*k116)) + (delta*(I2(k)+h2*k116)) - tau*(R(k)+h2*k116) - miu1*(R(k)+h2*k116);
k127 = tau*(R(k)+h2*k117) - miu1*(K(k)+h2*k117);
k131 = a - n*(S1(k)+h2*k121) - alpha*(S1(k)+h2*k121)*((I1(k)+h2*k121)+(I2(k)+h2*k121)) - miu1*(S1(k)+h2*k121);
k132 = n*(S1(k)+h2*k122) - alpha*(S2(k)+h2*k122)*((I2(k)+h2*k122)+(I1(k)+h2*k122)) - miu1*(S2(k)+h2*k122);
k133 = alpha*(S1(k)+h2*k123)*((I1(k)+h2*k123)+(I2(k)+h2*k123)) + gamma*(V(k)+h2*k123)*((I1(k)+h2*k123)+(I2(k)+h2*k123)) - (delta*(I1(k)+h2*k123)) - (miu1+miu2)*(I1(k)+h2*k123);
k134 = alpha*(S2(k)+h2*k124)*((I2(k)+h2*k124)+(I1(k)+h2*k124)) - (delta*(I2(k)+h2*k124)) - (miu1+miu2)*(I2(k)+h2*k124);
k135 = - gamma*(V(k)+h2*k125)*((I1(k)+h2*k125)+(I2(k)+h2*k125)) - miu1*(V(k)+h2*k125);
k136 = (delta*(I1(k)+h2*k126)) + (delta*(I2(k)+h2*k126)) - tau*(R(k)+h2*k126) - miu1*(R(k)+h2*k126);
k137 = tau*(R(k)+h2*k127) - miu1*(K(k)+h2*k127);
k141 = a - n*(S1(k)+h2*k131) - alpha*(S1(k)+h2*k131)*((I1(k)+h2*k131)+(I2(k)+h2*k131)) - miu1*(S1(k)+h2*k131);
k142 = n*(S1(k)+h2*k132) - alpha*(S2(k)+h2*k132)*((I2(k)+h2*k132)+(I1(k)+h2*k132)) - miu1*(S2(k)+h2*k132);
k143 = alpha*(S1(k)+h2*k133)*((I1(k)+h2*k133)+(I2(k)+h2*k133)) + gamma*(V(k)+h2*k133)*((I1(k)+h2*k133)+(I2(k)+h2*k133)) - (delta*(I1(k)+h2*k133)) - (miu1+miu2)*(I1(k)+h2*k133);
k144 = alpha*(S2(k)+h2*k134)*((I2(k)+h2*k134)+(I1(k)+h2*k134)) - (delta*(I2(k)+h2*k134)) - (miu1+miu2)*(I2(k)+h2*k134);
k145 = - gamma*(V(k)+h2*k135)*((I1(k)+h2*k135)+(I2(k)+h2*k135)) - miu1*(V(k)+h2*k135);
k146 = (delta*(I1(k)+h2*k136)) + (delta*(I2(k)+h2*k136)) - tau*(R(k)+h2*k136) - miu1*(R(k)+h2*k136);
k147 = tau*(R(k)+h2*k137) - miu1*(K(k)+h2*k137);
S11(k+1) = S11(k) + (h/6)*(k111 + 2*k121 + 2*k131 + k141);
S21(k+1) = S21(k) + (h/6)*(k112 + 2*k122 + 2*k132 + k142);
I11(k+1) = I11(k) + (h/6)*(k113 + 2*k123 + 2*k133 + k143);
I21(k+1) = I21(k) + (h/6)*(k114 + 2*k124 + 2*k134 + k144);
V1(k+1) = V1(k) + (h/6)*(k115 + 2*k125 + 2*k135 + k145);
R1(k+1) = R1(k) + (h/6)*(k116 + 2*k126 + 2*k136 + k146);
K1(k+1) = K1(k) + (h/6)*(k117 + 2*k127 + 2*k137 + k147);
end
temp1 = ((Lambda1-Lambda5).*S1)./(C1);
ua1 = min(max(temp1,0),1);
u1 = 0.5*(ua1 + oldu1);
temp2 = ((Lambda3+Lambda4-2.*Lambda6).*k)./(C1);
ua2 = min(max(temp2,0),1);
u2 = 0.5*(ua2 + oldu2);
temp3 = delta1*sum(abs(u1))-sum(abs(u1-oldu1));
temp4 = delta1*sum(abs(u2))-sum(abs(u2-oldu2));
temp5 = delta1*sum(abs(S1))-sum(abs(S1-oldS1));
temp6 = delta1*sum(abs(S2))-sum(abs(S2-oldS2));
temp7 = delta1*sum(abs(I1))-sum(abs(I1-oldI1));
temp8 = delta1*sum(abs(I2))-sum(abs(I2-oldI2));
temp9 = delta1*sum(abs(V))-sum(abs(V-oldV));
temp10 = delta1*sum(abs(R))-sum(abs(R-oldR));
temp11 = delta1*sum(abs(K))-sum(abs(K-oldK));
temp12 = delta1*sum(abs(Lambda1))-sum(abs(Lambda1-oldLambda1));
temp13 = delta1*sum(abs(Lambda2))-sum(abs(Lambda2-oldLambda2));
temp14 = delta1*sum(abs(Lambda3))-sum(abs(Lambda3-oldLambda3));
temp15 = delta1*sum(abs(Lambda4))-sum(abs(Lambda4-oldLambda4));
temp16 = delta1*sum(abs(Lambda5))-sum(abs(Lambda5-oldLambda5));
temp17 = delta1*sum(abs(Lambda6))-sum(abs(Lambda6-oldLambda6));
temp18 = delta1*sum(abs(Lambda7))-sum(abs(Lambda7-oldLambda7));
test = min(temp1,min(temp2,min(temp3,min(temp4,min(temp5,min(temp6,min(temp7,min(temp8,min(temp9,min(temp10,min(temp11,min(temp12,min(temp13,min(temp14,min(temp15,min(temp16,min(temp17,temp18)))))))))))))))));
end
y1(1,:) = t;
y1(2,:) = u1;
y1(3,:) = u2;
y1(4,:) = S1;
y1(5,:) = S2;
y1(6,:) = I1;
y1(7,:) = I2;
y1(8,:) = R;
y1(9,:) = V;
y1(10,:) = K;
y1(11,:) = S11;
y1(12,:) = S21;
y1(13,:) = I11;
y1(14,:) = I21;
y1(15,:) = R1;
y1(16,:) = V1;
y1(17,:) = K1;
y1(18,:) = Lambda1;
y1(19,:) = Lambda2;
y1(20,:) = Lambda3;
y1(21,:) = Lambda4;
y1(22,:) = Lambda5;
y1(23,:) = Lambda6;
y1(24,:) = Lambda7;
figure(1)
plot(y1(1,:),y1(2,:),'blue',y1(1,:),y1(3,:),'--green','LineWidth',2)
xlabel('time (month)')
ylabel('Optimal Control')
legend('Value u1','Value u2')
grid on;
but in workspace window the ans k11 and etc is a NaN, So how can I solve this? If any body could please suggest me if possible,
Thanks

Answers (1)

the cyclist
the cyclist on 27 Nov 2021
Your code is long and complicated. It also has lots of unnecessary complication in it. For example, instead of
Lambda1 = zeros(1,N+1); Lambda1(N+1) = 0;
Lambda2 = zeros(1,N+1); Lambda2(N+1) = 0;
Lambda3 = zeros(1,N+1); Lambda3(N+1) = 0;
Lambda4 = zeros(1,N+1); Lambda4(N+1) = 0;
Lambda5 = zeros(1,N+1); Lambda5(N+1) = 0;
Lambda6 = zeros(1,N+1); Lambda6(N+1) = 0;
Lambda7 = zeros(1,N+1); Lambda7(N+1) = 0;
you could have defined a single variable
Lambda = zeros(7,N+1);
and used indexing instead of seven different variables.
Combined with the fact that we do not have any idea of what these variables mean, or what you are calculating, it becomes extremely difficult for an "outsider" to debug this code for you. You may get some kind soul who is willing to put in that time, but I will instead just point you to the MATLAB tools for debugging this yourself.
Using these tools, you should be able to trace backward from your final output, to see where the NaN values are starting to get produced.

Products


Release

R2015a

Community Treasure Hunt

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

Start Hunting!