So i have this code and it gives me the error "Unable to perform assignment because the left and right sides have a different number of elements".
In this part of the code, %update the new approximation
LSc(j-1)=LSc(j)-h6*(auxLSc1+2*(auxLSc2+auxLSc3)+auxLSc4); %here is the error!!!
LSa(j-1)=LSa(j)-h6*(auxLSa1+2*(auxLSa2+auxLSa3)+auxLSa4); %here is the error!!
LIc(j-1)=LIc(j)-h6*(auxLIc1+2*(auxLIc2+auxLIc3)+auxLIc4); %here is the error!
LIa(j-1)=LIa(j)-h6*(auxLIa1+2*(auxLIa2+auxLIa3)+auxLIa4); %here is the error!
I enclose the entire code for more clarity. I understand what the error means but i can't see the problem in my code! sorry i'm rather new to programming.
function dy=odeControlage(T)
%parameters of the model
g=0.2; A=0.00018265; bcc=1; bca=0.2; bac=0.2; baa=1; g=0.2; C=0.001; K=1000; B=1; %balancing parameter
%parameters of the Runge-Kutta
test=-1; deltaError=0.001; M=100;
t=linspace(0,T,M+1);
h=T/M; h2=h/2; h6=h/6;
Sa=zeros(1,M+1); Sa=zeros(1,M+1); Ic=zeros(1,M+1); Ia=zeros(1,M+1);
Rc=zeros(1,M+1); Ra=zeros(1,M+1);
%initial conditions of the model
Sc(1)=0.199; Sa(1)=0.7; Ic(1)=0.001; Ia(1)=0.1 ; Ra(1)=0; Rc=0;
% vectors for system restrictions and control
LSc=zeros(1,M+1); LSa=zeros(1,M+1); LIc=zeros(1,M+1);LIa=zeros(1,M+1); U=zeros(1,M+1);
%iterations of the method
while (test<0)
oldSc=Sc; oldSa=Sa; oldIc=Ic; oldIa=Ia; oldRc=Rc; oldRa=Ra;
oldLSc=LSc; oldLSa=LSa; oldLIc=LIc; oldLIa=LIa;
oldU=U;
%forward Runge-Kutta Method
for i=1:M
%differential equations of the state system
%1st Runge-Kutta parameter
auxSc1=-bcc*U(i)*Sc(i)*Ic(i)-bca*U(i)*Sc(i)*Ia(i)-A*Sc(i);
auxSa1=-bac*U(i)*Sa(i)*Ic(i)-baa*U(i)*Sa(i)*Ia(i)+A*Sc(i);
auxIc1=bcc*U(i)*Sc(i)*Ic(i)+bca*U(i)*Sc(i)*Ia(i)-g*Ic(i)-A*Ic(i);
auxIa1=bac*U(i)*Sa(i)*Ic(i)+baa*U(i)*Sa(i)*Ia(i)-g*Ia(i)+A*Ic(i);
auxRc1=g*Ic(i)-A*Rc(i);
auxRa1=g*Ia(i)+A*Rc(i);
%2nd order Runge-Kutta
auxU=0.5* (U(i)+U(i+1));
auxSc=Sc(i)+ h2*auxSc1;
auxSa=Sa(i)+ h2*auxSa1;
auxIc=Ic(i)+ h2*auxIc1;
auxIa=Ia(i)+ h2*auxIa1;
auxRc=Rc(i)+ h2*auxRc1;
auxRa=Ra(i)+ h2*auxRa1;
%implement the update
auxSc2=-bcc*auxU*auxSc*auxIc -bca*auxU*auxSc*auxIa -A*auxSc;
auxSa2=-bac*auxU*auxSa*auxIa -baa*auxU*auxSa*auxIa +A*auxSc;
auxIc2=bcc*auxU*auxSc+bca*auxU*auxSc*auxIa - g*auxIc -A*auxIc;
auxIa2=bac*auxU*auxSa*auxIc + baa*auxU*auxSa*auxIa -g*auxIa + A*auxIc;
auxRc2=g*auxIc - A*auxRc;
auxRa2=g*auxIa + A*auxRc;
%3rd order Runge-Kutta
auxSc=Sc(i)+ h2*auxSc2;
auxSa=Sa(i)+ h2*auxSa2;
auxIc=Ic(i)+ h2*auxIc2;
auxIa=Ia(i)+ h2*auxIa2;
auxRc=Rc(i)+ h2*auxRc2;
auxRa=Ra(i)+ h2*auxRa2;
%update
auxSc3=-bcc*auxU*auxSc*auxIc -bca*auxU*auxSc*auxIa - A*auxSc;
auxSa3=-bac*auxU*auxSa*auxIa -baa*auxU*auxSa*auxIa +A*auxSc;
auxIc3=bcc*auxU*auxSc +bca*auxU*auxSc*auxIa -g*auxIc -A*auxIc;
auxIa3=bac*auxU*auxSa*auxIc + baa*auxU*auxSa*auxIa - g*auxIa + A*auxIc;
auxRc3=g*auxIc - A*auxRc;
auxRa3=g*auxIa + A*auxRc;
%4th order Runge-Kutta
auxSc=Sc(i)+ h*auxSc3;
auxSa=Sa(i)+ h*auxSa3;
auxIc=Ic(i)+ h*auxIc3;
auxIa=Ia(i)+ h*auxIa3;
auxRc=Rc(i)+ h*auxRc3;
auxRa=Ra(i)+ h*auxRa3;
%update
auxSc4=-bcc*U(i+1)*auxSc*auxIc -bca*U(i+1)*auxSc*auxIa- A*auxSc;
auxSa4=-bac*U(i+1)*auxSa*auxIc -baa*U(i+1)*auxSa*auxIa +A*auxSc;
auxIc4=bcc*U(i+1)*auxSc*auxIc +bca*U(i+1)*auxSc*auxIa-g*auxIc -A*auxIc;
auxIa4=bac*U(i+1)*auxSa*auxIc + baa*U(i+1)*auxSa*auxIa -g*auxIa +A*auxIc;
auxRc4=g*auxIc - A* auxRc;
auxRa4=g*auxIa + A*auxRc;
%Runge-Kutta New approximation
Sc(i+1)=Sc(i)+h6*(auxSc1+ 2*(auxSc2+auxSc3)+auxSc4);
Sa(i+1)=Sa(i)+h6*(auxSa1+ 2*(auxSa2+auxSa3)+auxSa4);
Ic(i+1)=Ic(i)+h6*(auxIc1+ 2*(auxIc2+auxIc3)+auxIc4);
Ia(i+1)=Ia(i)+h6*(auxIa1+ 2*(auxIa2+auxIa3)+auxIa4);
Rc(i+1)=Rc(i)+h6*(auxRc1+ 2*(auxRc2+auxRc3)+auxRc4);
Ra(i+1)=Ra(i)+h6*(auxRa1+ 2*(auxRa2+auxRa3)+auxRa4);
end
%Backward Runge-Kutta
%adjoint system
for i=1:M
j=M+2-i;
%update
auxLSc1=LSc(j)*(bcc*U(j)*Ic(j)+ bca*Ia(j)+A)-LIc(j)*(bcc*U(j)*Ic(j)+bca*U(j)*Ia(j))-LSa(j)*A;
auxLSa1=LSc(j)*(bac*U(j)*Ic(j)+baa*U(j)*Ia(j))-LIa(j)*(bac*U(j)*Ic(j)+baa*U(j)*Ia(j)-g);
auxx=B*K*exp(K*(C-(Ic(j)+Ia(j)))); %IS IT COMPLICATED MORE LIKE THIS?
auxLIc1=auxx + LSc(j)*bcc*U(j)*Sc(j)- LIc(j)*(bcc*U(j)*Sc(j)+bca*U(j)*Sc(j)-g -A)+ LSa(j)*bac*U(j)*Sa(j)+LIa*(-bac*U(j)*Sa(j)-A);
auxLIa1=auxx + LSc(j)*U(j)*bca*Sc(j)- LIc(j)*U(j)*bca*Sc(j)+ LSa(j)*baa*U(j)*Sa(j)-LIa(j)*(baa*Sa(j)-g);
%2nd order Runge-Kutta
auxU=0.5*(U(j)+U(j-1));
auxSc=0.5*(Sc(j)-Sc(j-1));
auxSa=0.5*(Sa(j)-Sa(j-1));
auxIc=0.5*(Ic(j)+Ic(j-1));
auxIa=0.5*(Ia(j)-Ia(j-1));
auxRc=0.5*(Rc(j)-Rc(j-1));
auxLSc=LSc(j)-h2*auxLSc1;
auxLSa=LSa(j)-h2*auxLSa1;
auxLIc=LIc(j)-h2*auxLIc1;
auxLIa=LIa(j)-h2*auxLIa1;
auxx=B*K*exp(K*(C-(Ic(j)+Ia(j))));
%update
auxLSc2=(bcc*auxU*auxIc+bca*auxU*auxIa +A)*auxLSc-auxLIc*(bcc*auxU*auxIc +bca*auxU*auxIa)-auxLSa*A;
auxLSa2=(bac*auxU*auxIc + baa*auxU*auxIa)*auxLSa- (bac*auxU*auxIc+baa*auxU*auxIa-g)*auxLIa;
auxLIc2=auxx + auxLSc*bcc*auxU*auxSc - auxLIc*(bcc*auxU*auxSc+bca*auxU*auxSc-g-A)+auxLSa*bac*auxU*auxSa + auxLIa*(-bac*auxU*auxSa-A);
auxLIa2=auxx + auxLSc*bca*auxU*auxSc - auxLIc*bca*auxU*auxSc + auxLSa*baa*auxU*auxSa - auxLIa*(baa*auxU*auxSa-g);
%3rd orde Runge-Kutta
auxLSc=LSc(j)-h2*auxLSc2;
auxLSa=LSa(j)-h2*auxLSa2;
auxLIc=LIc(j)-h2*auxLIc2;
auxLIa=LIa(j)-h2*auxLIa2;
%update
auxLSc3=(bcc*auxU*auxIc+bca*auxU*auxIa +A)*auxLSc-auxLIc*(bcc*auxU*auxIc +bca*auxU*auxIa)-auxLSa*A;
auxLSa3=(bac*auxU*auxIc + baa*auxU*auxIa)*auxLSa- (bac*auxU*auxIc+baa*auxU*auxIa-g)*auxLIa;
auxx=B*K*exp(K*(C-(Ic(j)+Ia(j))));
auxLIc3=auxx + auxLSc*bcc*auxU*auxSc - auxLIc*(bcc*auxU*auxSc+bca*auxU*auxSc-g-A)+auxLSa*bac*auxU*auxSa + auxLIa*(-bac*auxU*auxSa-A);
auxLIa3=auxx + auxLSc*bca*auxU*auxSc - auxLIc*bca*auxU*auxSc + auxLSa*baa*auxU*auxSa - auxLIa*(baa*auxU*auxSa-g);
%4th order Runge-Kutta
auxU=U(j-1);
auxSc=Sc(j-1);
auxSa=Sa(j-1);
auxIc=Ic(j-1);
auxIa=Ia(j-1);
auxRc=Rc(j-1);
auxRa=Ra(j-1);
auxLSc=LSc(j)-h*auxLSc3;
auxLSa=LSa(j)-h*auxLSa3;
auxLIc=LIc(j)-h*auxLIc3;
auxLIa=LIa(j)-h*auxLIa3;
%update
auxLSc4=(bcc*auxU*auxIc+bca*auxU*auxIa +A)*auxLSc-auxLIc*(bcc*auxU*auxIc +bca*auxU*auxIa)-auxLSa*A;
auxLSa4=(bac*auxU*auxIc + baa*auxU*auxIa)*auxLSa- (bac*auxU*auxIc+baa*auxU*auxIa-g)*auxLIa;
auxx=B*K*exp(K*(C-(Ic(j)+Ia(j))));
auxLIc4=auxx + auxLSc*bcc*auxU*auxSc - auxLIc*(bcc*auxU*auxSc+bca*auxU*auxSc-g-A)+auxLSa*bac*auxU*auxSa + auxLIa*(-bac*auxU*auxSa-A);
auxLIa4=auxx + auxLSc*bca*auxU*auxSc - auxLIc*bca*auxU*auxSc + auxLSa*baa*auxU*auxSa - auxLIa*(baa*auxU*auxSa-g);
%update the new approximation
LSc(j-1)=LSc(j)-h6*(auxLSc1+2*(auxLSc2+auxLSc3)+auxLSc4); %here is the error!!!
LSa(j-1)=LSa(j)-h6*(auxLSa1+2*(auxLSa2+auxLSa3)+auxLSa4); %here is the error!!
LIc(j-1)=LIc(j)-h6*(auxLIc1+2*(auxLIc2+auxLIc3)+auxLIc4); %here is the error!
LIa(j-1)=LIa(j)-h6*(auxLIa1+2*(auxLIa2+auxLIa3)+auxLIa4); %here is the error!
end
%new control vector
for i=1:M+1
vAux(i)=0.5*(bcc+bca+bac+baa+bcc*LSc(i)*Sc(i)*Ic(i)+bca*Sc(i)*Ia(i)-LIc(i)*(bcc*Sc(i)*Ic(i)+bca*Sc(i)*Ia(i))+LSa(i)*(bac*Sa(i)*Ic(i)+baa*Sc(i)*Ia(i))-LIa(i)*(bac*Sa(i)*Ic(i)+baa*Sa(i)*Ia(i)));
auxU = min([max([0.0 vAux(i)]) 1]);
U(i) = 0.5 * (auxU + oldU(i));
end
%absolute error for convergence
temp1=deltaError*sum(abs(Sc))-sum(abs(oldSc-Sc));
temp2=deltaError*sum(abs(Sa))-sum(abs(oldSa-Sa));
temp3=deltaError*sum(abs(Ic)-sum(abs(oldIc-Ic)));
temp4=deltaError*sum(abs(Ia))-sum(abs(oldIa-Ia));
temp5=deltaError*sum(abs(U))-sum(abs(oldU-U));
temp6=deltaError*sum(abs(LSc))-sum(abs(oldLSc-LSc));
temp7=deltaError*sum(abs(LSa))-sum(abs(oldLSa-LSa));
temp8=deltaError*sum(abs(LIc))-sum(abs(oldLIc-LIc));
temp9=deltaError*sum(abs(Rc))-sum(abs(oldRc-Rc));
temp10=deltaError*sum(abs(Ra))-sum(abs(oldRa-Ra));
test = min(temp1,min(temp2,min(temp3,min(temp4,min(temp5,min(temp6,min(temp7,min(temp8,min(temp9,min(temp10))))))))));
end
dy(1,:) = t; dy(2,:) = Sc; dy(3,:) = Sa;
dy(4,:) = Ic; dy(5,:) = Rc; dy(6,:) = Ra; dy(7,:)=U;
figure(2);
plot(t,Rc,"m");
hold off;disp("Value of LAMBDA at FINAL TIME");
disp([LS(M+1) LE(M+1) LI(M+1)]);
end
