# Index exceeds the number of array elements (0)

13 views (last 30 days)
cRISTIAN Almeida Ortiz on 5 Dec 2019
Answered: the cyclist on 5 Dec 2019
clear
clc
i=100000;
N=zeros(i:1);
s=zeros(i:1);
alfa=zeros(i:1);
Vs=zeros(i:1);
k2=zeros(i:1);
Va=zeros(i:1);
Vr=zeros(i:1);
t=zeros(i:1);
%parametros
kB=1.36e-25; %constante dealfa (i:1); boltzmann atm L/K
T=291; % temperatura absoluta ambiente de quito K
u=1.18e-05; % viscosidad del suero de leche atm*s
E= 99.23; %constante dielectrica
a=500; %radio de micelio de caseina A
Yo=0.012; %potencial de la superficie Volts
kr=0.1; %es el reciproco de espesor de la doble capa A-1
H=1e-22; %constante de Hamaker para micelas J
ko2=4*kB*T/(3*u); % constante cinetica aparente
deltat=0.01;
N(1:1)=10000000;
s(1:1)=2;
t(1:1)=0;
Vs(1:1)=Vr(1:1)+Va(1:1);
Vr(1:1)=((E*a*Yo^2)/2)*log(1+ exp(-kr*a*(s(1:1)-2)));
Va(1:1)=(-H/6)*((2/(S(1:1)^2-4))+(2/S(1:1)^2) + log((S(1:1)^2-4)/s(1:1)^2));
for k=1:i
funcionN=-k2(k)*(N(k)^2); %ec1
N(k+1)= N(k)+ deltat*funcionN;
k2(k+1)= ko2/(2*alfa(k+1)); %ec2
funcionalfa=exp((Vs(k)/(kB*T))/s(k)^2); %ec3
alfa(k+1)=alfa(k)+deltat*funcionalfa;
Vs(k+1)=Vr(k+1)+Va(k+1); %ec4
Vr(k+1)=((E*a*Yo^2)/2)*log(1+ exp(-kr*a*(s(k+1)-2)));
Va(k+1)=(-H/6)*((2/(S(k+1)^2-4))+(2/S(k+1)^2) + log((S(k+1)^2-4)/s(k+1)^2));
t(k+1)=t(k)+deltat;
end
figure(1) % prepara para crear algo y y x
plot(t,N, 'LineWidth',1.5) % line width da el tamaño de linea
xlabel('tiempo en h') % la etiqueta del eje x
ylabel ('particulas independientes')
grid on % permite cuadricular el grafico
hold on %asi entiende que viene otra figura

the cyclist on 5 Dec 2019
Everywhere that you wrote
zeros(i:1)
I think you meant
zeros(i,1)