I am having trouble trying to maintain the loop

Nt=161224;
dt=6e-12;
T=9.6734e-7;
Nz=8061;
Q=2.3131e-6;
dz=0.0012;
N=8061;
n=1.45;
gammaE=0.8967;
Z0=376.7343;
A=8e-11;
z=0:dz:(dz*N-dz);
c=3e8;
G=255.747;
Gth=21;

So basically I am trying ti solve this coupled PDE's , here I have defined ELt as 8061x8061 matrix where time increments is along rows and space along clumns .The program solves for each time iteration and I am trying to save all the values in the matrix.The problem arises when index=8061 and for 8062 it should go to the first row bit I am unable to figure it out. The program just says index array is exceeded. Can anyone help with this loop...

q= 1.7125e07;
eps0=8.854e-12;
Omega=1.207e11;
sigma=826.2189;
GAMMAb=3.5904e8;
Nt1=round(Nt/20);
%t = Nt*dt;
t = (-T/2/dt:1:T/2/dt)*dt;
%fs=1/dbt;
fsine = 1e9;
vsine = 1;
phi = vsine*sin(2*pi*fsine*t);



%% initial conditions

    
    
    EL1t=1.5607e7*exp(1i*phi);
    %EL1t=El_s(1)*ones(1,round(Nt));
    ELt=zeros(Nt1,Nz);
    ESt=zeros(Nt1,Nz);
    ELnew=zeros(1,Nz);
    ESnew=ELnew;
     
    %roh=sqrt(n*Q/c/GAMMAb)*(randn(1,Nz)+1i*randn(1,Nz));
    %roh=2/GAMMAb*sqrt(Q*GAMMAb/L)*(randn(1,Nz)*(1+1i));
    
    %roh1=2/GAMMAb*sqrt(Q*GAMMAb/L)*(randn(1,Nz)*(1+1i)); %roh value for previous time instant
    roh1 = wgn(1,Nz,Q/dt/dz,'complex','linear')*2/(GAMMAb);
    roh2=roh1;    %roh value for current time instant
    

%
%% Numerical Process

for index=1:Nt
  

    %ELt(1)=EL1t(index);
    ELnew(1)=EL1t(index);  %EL(0,t)=ELoexp(i*phi)
  
    Elz=(ELt(index,2:end)-ELt(index,1:end-1))/dz;    
    ELnew(2:end)=ELt(index,2:end)+dt*c/n*(1i*sigma*roh2(2:end).*ESt(index,2:end)-Elz);%final equation for EL
    Esz=(ESt(index,2:end)-ESt(index,1:end-1))/dz;
    croh=conj(roh2);
    ESnew(1:end-1)=ESt(index,1:end-1)+dt*c/n*(1i*sigma*croh(1:end-1).*ELt(index,1:end-1)+Esz);  %final equation for ES

    %f=sqrt(n*Q/dt^2/c)*(randn(1,Nz)*(1+1i));
    %f =sqrt((Q*GAMMAb)/L)*(randn(1,Nz)+1i*randn(1,Nz));
    f =wgn(1,Nz,Q/dt/dz,'complex','linear');
    
    
    %rohnew=roh+dt*p;    this and below equations are for second order
    %pnew=p+dt*((-GAMMAb+2*1i*Omega)*p+1i*Omega*GAMMAb*roh+eps0*gammaE*q^2*ELt.*conj(ESt)-2*1i*Omega*f); 
    %rohnew=roh+dt*(-GAMMAb/2*roh+1i*LAMBDA*ELt.*conj(ESt)+f);  the equation is for first order equation(ignoring GAMMAb) 
    rohnew=1/(1/dt+GAMMAb-2*1i*Omega)*(-roh1/dt+roh2*(2/dt+GAMMAb-2*1i*Omega+1i*GAMMAb*Omega*dt)+dt*(gammaE*eps0*q^2*ELt(index,:).*conj(ESt(index,:))-2*1i*Omega*f));
    ELt(index,:)=ELnew;
    ESt(index,:)=ESnew;
    roh1=roh2;
    roh2=rohnew;
    %p=pnew;
         
% % %     Plotting the longitudial ntensities for following the process

    
    
   if (mod(index,5000)==0)
        if G<Gth
            
            subplot(2,1,1);
            plot (z,2*n/Z0*A*abs(ELt).^2);
            
                hold on
                %I_L=2*n/Z0*A*abs(El_s).^2;
                legend('Dynamic','Steady State')
                hold off
            title(['Iteration Number:   ',num2str(index)])
            xlabel('z [m]'); ylabel('I_L [Watt]');
            subplot(2,1,2);
            plot (z,2*n/Z0*A*abs(ESt).^2);
               
                
                hold on
                
                legend('Dynamic','Steady State')
            ylabel('I_S [Watt]'); xlabel('z [m]');
        else

            subplot(2,1,1);  
            plot (z,2*n/Z0*A*abs(ELt).^2)
            hold on
            %plot(z,I_L)
            legend('Pump Dynamic','Pump Steady State')
            hold off
            title(['Iteration Number:   ',num2str(index),' of ' num2str(Nt)])
            xlim([-1 11]); %ylim([-2 inf]);
            xlabel('z [m]');    ylabel ('Optical Power')
                        
                 subplot(2,1,2);
                  plot (z,2*n/Z0*A*abs(ESt).^2)
                  hold on
                 % plot(z,I_S)
                  legend('Stokes Dynamic','Stokes Steady State')
                  title(['Iteration Number:   ',num2str(index),' of ' num2str(Nt)])
                  xlim([-1 11]); %ylim([-2 inf]);
                  xlabel('z [m]');    ylabel ('Optical Power')
                  
              %figure;
             %plot(z,2*n/Z0*A*abs(ELt).^2);

        end 
        hold off
        drawnow;
   
    end
    





end

Index in position 1 exceeds array bounds. Index must not exceed 8061.

 if G<Gth  
     subplot(2,1,1);
 end
 
title('Intensity distibution along the fiber');

[Stephen23]

"The program just says index array is exceeded. Can anyone help with this loop..."

Here are the relevant lines of code:

Nt=161224;
...
Nt1=round(Nt/20);
...
ELt=zeros(Nt1,Nz);
...
for index=1:Nt
    ...	
    Elz=(ELt(index,2:end)-ELt(index,1:end-1))/dz;

So you preallocate an array ELT with 8601 rows. But your loop iterates from 1 to 161224, thus throwing an error when you try to access row 8602 (which does not exist).

"The problem arises when index=8061 and for 8062 it should go to the first row bit I am unable to figure it out"

It is unclear exactly what you mean, but perhaps you could use a modulo operation, e.g.:

    idx = 1+mod(index-1,Nt1)
    Elz = (ELt(idx,2:end)-ELt(idx,1:end-1))/dz;

[Yogesh]

Hey thats a good idea!!!...

Thank you. I got it now!!

Answer 1

[VBBV]

@Yogesh when array sizes dont match , those errors occur. check the array dimensions in for loop

Nt=161224;
dt=6e-12;
T=9.6734e-7;
Nz=8061

Nz = 8061

Q=2.3131e-6;
dz=0.0012;
N=8061;
n=1.45;
gammaE=0.8967;
Z0=376.7343;
A=8e-11;
z=0:dz:(dz*N-dz);
c=3e8;
G=255.747;
Gth=21;
q= 1.7125e07;
eps0=8.854e-12;
Omega=1.207e11;
sigma=826.2189;
GAMMAb=3.5904e8;

%t = Nt*dt;
t = (-T/2/dt:1:T/2/dt)*dt

t = 1x161224

1.0e-06 * -0.4837 -0.4837 -0.4837 -0.4837 -0.4836 -0.4836 -0.4836 -0.4836 -0.4836 -0.4836 -0.4836 -0.4836 -0.4836 -0.4836 -0.4836 -0.4836 -0.4836 -0.4836 -0.4836 -0.4836 -0.4836 -0.4835 -0.4835 -0.4835 -0.4835 -0.4835 -0.4835 -0.4835 -0.4835 -0.4835

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

%fs=1/dbt;
fsine = 1e9;
vsine = 1;
phi = vsine*sin(2*pi*fsine*t);
%% initial conditions

Nt1=round(Nt/20)

Nt1 = 8061

EL1t=1.5607e7*exp(1i*phi);
%EL1t=El_s(1)*ones(1,round(Nt));

ELt=zeros(Nt1,Nz);  % the size of this array 8061 x  161224
ESt=zeros(Nt1,Nz);  % the size of this array 8061 x  161224

ELnew=zeros(1,Nt);  % here the size needs to be checked 

ESnew=ELnew;

% in your code
for index = 1:1000:Nt
  ELt(index,:)=ELnew;  % Here the value of index exceeds the 8061 ...iterates upto 161224
  ESt(index,:)=ESnew;
end

ELnew=zeros(Nt1,1);  

% change needed for array size match 
for index = 1:1000:Nt
  ELt(:,index)=ELnew;  % Here the value of index exceeds the 8061 
  ESt(:,index)=ESnew;
end