Can anyone help me to rectify the error?

main.m

clear
 epsilon = 1/12;
 gamma = 1/15;
 beta = 0.2;
 theta1 = 0.0082;
 theta2 = 0.0289;
 mu = 1.01/(70*365);
 Lambda = 0.379/(70*365);
 rho = 205;
 delta0 = 1050;
 r= 1/7;
 r0 = 0.04;
 delta1 = 0.0399;
 d = 1/10;
 psi = 0.0014;
 zeta1 = 0;
 zeta2 = 0.67;
alpha = 0.0238;
phi = 3;
 A1 = 1;
 A2 = 1;
 B1 = 50;
 B2 = 50;
 T = 365; %final time
 N = 200;
 x0 = [100; 20; 20; 10; 1000; 20; 30]; %initial value of x
 u0 = [0; 0]; %initial value of u
 lambdafinal =[-0.1;0;0;-0.1;0;0;0]; %final values of lambda
 t = linspace(0,T,N+1);
 h = T/N;
 x = zeros(7,N+1); %initialize x
 u = zeros(2,N+1); %initialize u
 lambda = zeros(7,N+1); %initialize x
 x (:,1) = x0; %assign the x0
 u (:,1) = u0; %assign the u0
 lambda(:,N+1) = lambdafinal; %assign the lambdafinal

 %R0
 %1/2*(zeta1 * Lambda*d/(k1*k2) + zeta2 * delta1*epsilon/(gamma*k4))...
 %+sqrt((1/2 *(zeta1*Lambda*d/(k1*k3)?zeta2*delta1*epsilon/(gamma*k4)))^2...
 %+phi^2 * (beta*epsilon*(theta1*d+theta2*k3)*(delta0/(gamma?delta1)))/(gamma*
%k1*k3*k4*(rho/(mu-Lambda))))

 parameters=[epsilon,gamma,beta,theta1,theta2,mu,Lambda,rho,delta0,r, r0, delta1, d, psi , zeta1, zeta2,alpha,phi,A1,A2,B1,B2];
k = 1;%counter of the iteration
 delta = 0.001; %error bound
 test = -1; %error
 while ( test<0 && k<1000)
 oldx = x;
 oldu = u;
 oldlambda = lambda;

 %forward Runge?Kutta 4th with 3 input algorithm for state
 for i = 1:N
 k1(1:7,1) = state(t(i), x(:,i) ,u(:,i) ,parameters);
 k2(1:7,1) = state(t(i)+h/2,x(:,i )+h*k1/2,(u(:,i)+u(:,i+1))/2,parameters);
 k3(1:7,1) = state(t(i)+h/2,x(:,i )+h*k2/2,(u(:,i)+u(:,i+1))/2,parameters);
 k4(1:7,1) = state(t(i)+h,x(:, i )+h*k3,u(:,i+1),parameters);
 x (:,i+1) = x(:,i)+ (h/6)*(k1+2*k2+2*k3+k4);
 end

 %backward Runge?Kutta 4th with 3 input algorithm for adjoint
 for i = 1:N
 j = N+2-i;
 k1(1:7,1) = adjoint(t(j),lambda(:,j), x(:,j),u(:,j),parameters);
 k2(1:7,1) = adjoint(t(j)-h/2,lambda(:,j)-h*k1/2,(x(:,j)+x(:,j-1))/2 ,(u (:,j)+u(:,j-1))/2,parameters);
 k3(1:7,1) = adjoint(t(j)-h/2,lambda(:,j)-h*k2/2,(x(:,j)+x(:,j-1))/2 ,(u (:,j)+u(:,j-1))/2,parameters);
 k4(1:7,1) = adjoint(t(j)-h/2,lambda(:,j)-h*k3/2,x(:,j-1),u(:,j-1),parameters);
 lambda(:,j-1) = lambda(:,j) - (h/6)*(k1+2*k2+2*k3+k4);

 end
 

 %Find controls
 for i = 1:N+1
 u(1, i ) = max(0,min(1,(beta*lambda(1,i)*(-x(1,i))*x(7,i)*phi +beta*lambda(2,i)*x(1,i)*x(7,i)*phi ...
 +(lambda(6,i)-lambda(5,i))*x(5,i)*phi*(theta2*x(2,i) +theta1*x(3,i)))/(2*B1*(x(1,i)+x(2,i)+x(3,i)+x(4,i))))) ;
 u(2, i ) = max(0,min(1,(lambda(3,i)-lambda(4,i))*r0*x(3,i)/(2*B2)));
 end

 %updates Control
 c = 0.8;
 u = (1-c)*u+c*oldu;
 % for i=1:N+1
 % if (u(1, i ) > oldu(1,i))
 % u(1, i ) = (1 ? c) + oldu(1,i)*c;
 % else
 % u(1, i ) = oldu(1,i)*c;
 % end
 % if (u(2, i ) > oldu(2,i))
 % u(2, i ) = (1 ? c) + oldu(2,i)*c;
 % else
 % u(2, i ) = oldu(2,i)*c;
 % end
 % end

 %error bewteen old and new

 tempu = min(delta*sum(abs(u),2)-sum(abs(oldu-u),2));
 tempx = min(delta*sum(abs(x),2)-sum(abs(oldx-x),2));
 templambda = min(delta*sum(abs(lambda),2)-sum(abs(oldlambda-lambda),2));
 test = min(tempu,min(tempx,templambda));

 %The cost
 %trapz(A1*oldx(2,:)+A2*oldx(3,:)+B1*oldu(1,:).^2+B2*oldu(2,:).^2)...
 % ?0.1*oldx(1,N+1)? 0.1*oldx(4,N+1)
 %trapz(A1*x(2,:)+A2*x(3,:)+B1*u(1,:).^2+B2*u(2,:).^2) ...
 % ?0.1*x(1,N+1)?0.1*x(4,N+1)
 k=k+1;
 end
 plot(t,x(3,:),t,x(7,:))
title (' infectious ')
 legend('host','vector')

state.m
function dxdt = state(t,x,u,parameters)
parameters = num2cell(parameters);
[epsilon,gamma,beta,theta1,theta2,mu,Lambda,rho,delta0,r, r0, delta1, d, psi, zeta1, zeta2,alpha,phi] = deal(parameters{:};

Invalid expression. When calling a function or indexing a variable, use parentheses. Otherwise, check for mismatched delimiters.

 dxdt = zeros(7,1);
 dxdt(1) = rho + Lambda*(x(1)+x(2)+x(4)) + Lambda*(1-zeta1)*x(3)+psi*x(4) -beta* phi*x(7)*x(1)*(1-u(1))/(x(1)+x(2)+x(3)+x(4))-mu*x(1);
 dxdt(2) = beta*phi*x(7)*x(1)*(1-u(1))/(x(1)+x(2)+x(3)+x(4))+ Lambda*zeta1*x(3)-d*x(2)-mu*x(2);
 dxdt(3) = d*x(2) - (r + alpha + mu+r0*u(2))*x(3);
 dxdt(4) = (r+r0*u(2))* x(3)- (psi + mu)*x(4);
 dxdt(5) = delta0 + delta1*(x(5)+x(6)) + delta1*(1-zeta2)*x(7) - phi*theta1*x(3)*x(5)*(1-u(1))/(x(1)+x(2)+x(3)+x(4)) ...
 - phi*theta2*x(2)*x(5)*(1-u(1))/(x(1)+x(2)+x(3)+x(4)) - gamma*x(5);
 dxdt(6) = phi*theta1*x(3)*x(5)*(1-u(1))/(x(1)+x(2)+x(3)+x(4)) + phi*theta2*x(2)*x(5)*(1-u(1))/(x(1)+x(2)+x(3)+x(4))...
 + delta1*zeta2*x(7)-epsilon*x(6)-gamma*x(6);
 dxdt(7) = epsilon*x(6)-gamma*x(7);
adjoint.m
function dlambdadt = adjoint(t,lambda,x,u,parameters)
 parameters = num2cell(parameters);
 [epsilon,gamma,beta,theta1,theta2,mu,Lambda,r, r0, delta1, d, psi, zeta1, zeta2,alpha,phi,A1,A2] = deal(parameters{:});

 dlambdadt = zeros(7,1);
 dlambdadt(1) = -(lambda(1)*(Lambda-mu+beta*(1-u(1))*x(1)*x(7)* phi/(x(1)+x(2)+x(3)+x(4))^2-beta*(1-u(1))*x(7)*phi/(x(1)+x(2)+x(3)+x(4))) +lambda(2)*(beta*(1-u(1))*x(7)*phi/(x(1)+x(2)+x(3)+x(4))-beta*(1-u(1))*x(1)*x(7)*phi/(x(1)+x(2)+x(3)+x(4))^2 +lambda(5)*(theta2*(1-u(1))*x(2)*x(5)*phi/(x(1)+x(2)+x(3)+x(4))^2+theta1*(1-u(1))*x(3)*x(5)*phi/(x(1)+x(2)+x(3)+x(4))^2) +lambda(6)*(-theta2*(1-u(1))*x(2)*x(5)*phi/(x(1)+x(2)+x(3)+x(4))^2 -theta1*(1-u(1))*x(3)*x(5)*phi/(x(1)+x(2)+x(3)+x(4))^2);
 dlambdadt(2) = -(A1+d*lambda(3) +lambda(2)*(-d-mu-beta*(1-u(1))*x(1)*x(7)*phi/(x(1)+x(2)+x(3)+x(4))^2) +lambda(1)*(Lambda+beta*(1-u(1))*x(1)*x(7)*phi/(x(1)+x(2)+x(3)+x(4))^2) +lambda(5)*(theta2*(1-u(1))*x(2)*x(5)*phi/(x(1)+x(2)+x(3)+x(4))^2+theta1*(1-u(1))*x(3)*x(5)*phi/(x(1)+x(2)+x(3)+x(4))^2-theta2*(1-u(1))*x(5)*phi/(x(1)+x(2)+x(3)+x(4))) +lambda(6)*(-theta2*(1-u(1))*x(2)*x(5)*phi/(x(1)+x(2)+x(3)+x(4))^2 -theta1*(1-u(1))*x(3)*x(5)*phi/(x(1)+x(2)+x(3)+x(4))^2+theta2*(1-u(1))*x(5)*phi/(x(1)+x(2)+x(3)+x(4))));
 dlambdadt(3) = -(A2+lambda(3)*(-alpha-mu-r0*u(2)-r)+lambda(4)*(r0*u(2)+r) +lambda(1)*((1-zeta1)*Lambda+beta*(1-u(1))*x(1)*x(7)*phi/(x(1)+x(2)+x(3)+x(4))^2) +lambda(2)*(zeta1*Lambda-beta*(1-u(1))*x(1)*x(7)*phi/(x(1)+x(2)+x(3)+x(4))^2) +lambda(5)*(theta2*(1-u(1))*x(2)*x(5)*phi/(x(1)+x(2)+x(3)+x(4))^2 +theta1*(1-u(1))*x(3)*x(5)*phi/(x(1)+x(2)+x(3)+x(4))^2 -theta1*(1-u(1))*x(5)*phi/(x(1)+x(2)+x(3)+x(4))) +lambda(6)*(-theta2*(1-u(1))*x(2)*x(5)*phi/(x(1)+x(2)+x(3)+ x(4))^2 -theta1*(1-u(1))*x(3)*x(5)*phi/(x(1)+x(2)+x(3)+x(4))^2+theta1*(1-u(1))*x(5)*phi/(x(1)+x(2)+x(3)+x(4)));
 dlambdadt(4) = -(lambda(4)*(-mu-psi)+ lambda(1)*(Lambda+beta*(1-u(1))*x(1)*x(7)*phi/(x(1)+x(2)+x(3)+x(4))^2+psi) -lambda(2)*(beta*(1-u(1))*x(1)*x(7)*phi/(x(1)+x(2)+x(3)+x(4))^2) +lambda(5)*(theta2*(1-u(1))*x(2)*x(5)*phi/(x(1)+x(2)+x(3)+x(4))^2 +theta1*(1-u(1))*x(3)*x(5)*phi/(x(1)+x(2)+x(3)+x(4))^2) +lambda(6)*(-theta2*(1-u(1))*x(2)*x(5)*phi/(x(1)+x(2)+x(3)+x(4))^2) -theta1*(1-u(1))*x(3)*x(5)*phi/(x(1)+x(2)+x(3)+x(4))^2);
 dlambdadt(5) = -(lambda(5)*(-gamma+delta1-theta2*(1-u(1))*x(2)*phi/(x(1)+x(2)+x(3)+x(4))-theta1*(1-u(1))*x(3)*phi/(x(1)+x(2)+x(3)+x(4))) +lambda(6)*(theta2*(1-u(1))*x(2)*phi/(x(1)+x(2)+x(3)+x(4))+theta1*(1-u(1))*x(3)*phi/(x(1)+x(2)+x(3)+x(4))));
 dlambdadt(6) = -(lambda(6)*(-gamma-epsilon)+delta1*lambda(5)+lambda(7)*epsilon);
 dlambdadt(7) = -(-gamma*lambda(7) +delta1*(1-zeta2)*lambda(5)+ delta1*zeta2*lambda(6) -beta*lambda(1)*(1-u(1))*x(1)*phi/(x(1)+x(2)+x(3)+x(4)) +beta*lambda(2)*(1-u(1))*x(1)*phi/(x(1)+x(2)+x(3)+x(4)));

Answer 1

[Anudeep Kumar]

Hey Sudipa,

I ran your code on my end and was able to reproduce the errors you encountered. Based on the R2024b definition of the functions you used; I identified two main issues:

Firstly, for line 3 in “state.m” -and “adjoint.m” the use of the MATLAB function “deal()” requires the number of input arguments to be equal to the number of output arguments. You were passing all 22 arguments in your variable “parameters” which was leading to an error.

Assuming you want to pass the first 18 arguments the following change in state.m:

>>[epsilon,gamma,beta,theta1,theta2,mu,Lambda,rho,delta0,r, r0, delta1, d, psi, zeta1, zeta2,alpha,phi] = deal(parameters{1:18}); 

And the below in adjoint.m:

>>[epsilon,gamma,beta,theta1,theta2,mu,Lambda,r, r0, delta1, d, psi, zeta1, zeta2,alpha,phi,A1,A2] = deal(parameters{1:18}); 

This should resolve the error.

If you need to pass specific arguments, you can use cell indexing or create a new variable accordingly.

Here's a link to guide you on accessing data in a cell array:

https://www.mathworks.com/help/matlab/matlab_prog/access-data-in-a-cell-array.html

Also, here’s the documentation for the deal() function:

https://www.mathworks.com/help/matlab/ref/deal.html

Secondly, in the definition of “adjoint.m”, there was an imbalance in the number of opening and closing parentheses. Assuming the formula was missing parentheses at the end, the following changes resolved the issue:

At line 6:

>> dlambdadt(1) = -(lambda(1)*(Lambda-mu+beta*(1-u(1))*x(1)*x(7)* ...
    phi/(x(1)+x(2)+x(3)+x(4))^2-beta*(1-u(1))*x(7)*phi/(x(1)+x(2)+x(3)+x(4)))...
    +lambda(2)*(beta*(1-u(1))*x(7)*phi/(x(1)+x(2)+x(3)+x(4))...
    -beta*(1-u(1))*x(1)*x(7)*phi/(x(1)+x(2)+x(3)+x(4))^2 ...
    +lambda(5)*(theta2*(1-u(1))*x(2)*x(5)*phi/(x(1)+x(2)+x(3)+x(4))^2 ...
    +theta1*(1-u(1))*x(3)*x(5)*phi/(x(1)+x(2)+x(3)+x(4))^2) ...
    +lambda(6)*(-theta2*(1-u(1))*x(2)*x(5)*phi/(x(1)+x(2)+x(3)+x(4))^2 ...
    -theta1*(1-u(1))*x(3)*x(5)*phi/(x(1)+x(2)+x(3)+x(4))^2))); 

And at line 8 :

>>dlambdadt(3) = -(A2+lambda(3)*(-alpha-mu-r0*u(2)-r)+lambda(4)*(r0*u(2)+r)...
    +lambda(1)*((1-zeta1)*Lambda+beta*(1-u(1))*x(1)*x(7)*phi/(x(1)+x(2)+x(3)+x(4))^2)...
    +lambda(2)*(zeta1*Lambda-beta*(1-u(1))*x(1)*x(7)*phi/(x(1)+x(2)+x(3)+x(4))^2)...
    +lambda(5)*(theta2*(1-u(1))*x(2)*x(5)*phi/(x(1)+x(2)+x(3)+x(4))^2 ...
    +theta1*(1-u(1))*x(3)*x(5)*phi/(x(1)+x(2)+x(3)+x(4))^2 ...
    -theta1*(1-u(1))*x(5)*phi/(x(1)+x(2)+x(3)+x(4))) +lambda(6)*(-theta2*(1-u(1))*x(2)*x(5)*phi/(x(1)+x(2)+x(3)+ x(4))^2 ...
    -theta1*(1-u(1))*x(3)*x(5)*phi/(x(1)+x(2)+x(3)+x(4))^2 ...
    +theta1*(1-u(1))*x(5)*phi/(x(1)+x(2)+x(3)+x(4)))); 

If the formula is different, please review and ensure that all parentheses are correctly balanced.

I hope this helped resolve your errors!