No error but I am getting Complex double!

clear ALL
close ALL 
clc 

Area=0.1963;
L=4;
porosity=0.18;
K=61.5;
injrate=0.01;
oilvisc=4.2;
watervisc=0.85;
Boi=1.21;
Swr=0.28;
Sorw=0.35;
KrwMax=0.31;
KroMax=0.65;
Nw=2.25;
No=2.45;
CO=0.000001;
CW=0.000004;
CR=0.000003;

POi=1500;
PWF=1460;
dt=0.2;
Vr=Area*0.4;

for i=1:10
     Swi(1,i)=Swr;
     PO(1,i)=POi;
     PCOW(1,i)=5;
     
     Krw(1,i)=KrwMax*((Swi(1,i)-Swr)/(1-Sorw-Swr))^Nw;  
     Kro(1,i)=KroMax*((1-Swi(1,i)-Sorw)/(1-Sorw-Swr))^No;
     Omob(1,i)= Kro(1,i)/oilvisc;
     Wmob(1,i)=Krw(1,i)/watervisc;
    
end
for n=1:15
    
    for i=1:10
         Ctot(n,:)=CR+Swi(n,:)*CW+(1-Swi(n,:))*CO;

         if i==1
             
             Ttot=(2.6366e-4)*K*(Wmob(n,i)+Omob(n,i))*(Area/0.4);
             Tw=(2.6366e-4)*K*(Wmob(n,i))*(Area/0.4);
             
             A(i,i)=-(Ttot+(Vr*porosity*Ctot(n,i))/(dt)); %E
             A(i,i+1)=Ttot; %F
             R(i,1)= (-(Vr*porosity*Ctot(n,i)*PO(n,i))/(dt))+Tw*(PCOW(n,i+1)-PCOW(n,i))-injrate; %R
     
         elseif i==10
             WI(n,i)=(2.6366e-4)*K*(Wmob(n,i)+Omob(n,i))*(Area/(0.4/2));
             TtotL=Ttot;
             TwL=Tw;
             Ttot=0;
             Tw=0;
           
             %Ttot=2.6366*10^-4*K*(Wmob(n,i)+Omob(n,i))*(A/0.4);
             %Tw=2.6366*10^-4*K*(Wmob(n,i))*(A/0.4);
   
             
              A(i,i-1)=TtotL; %D
              A(i,i)=-(TtotL+((Vr*porosity*Ctot(n,i))/(dt))+WI(n,i)); %E
              R(i,1)=(-(Vr*porosity*Ctot(n,i)*PO(n,i))/(dt))-TwL*(PCOW(n,i)-PCOW(n,i-1))-WI(n,i)*PWF; %R

         else 
               TtotL=Ttot;
               TwL=Tw;
               Ttot=(2.6366e-4)*K*(Wmob(n,i)+Omob(n,i))*(Area/0.4);
               Tw=(2.6366e-4)*K*(Wmob(n,i))*(Area/0.4);
               
               A(i,i)=-(TtotL+Ttot+(Vr*porosity*Ctot(n,i))/dt); %E
               A(i,i+1)=Ttot; %F
               A(i,i-1)=TtotL; %D
               R(i,1)= (-(Vr*porosity*Ctot(n,i)*PO(n,i))/(dt))+Tw*(PCOW(n,i+1)-PCOW(n,i))-TwL*(PCOW(n,i)-PCOW(n,i-1)); %R

        
         end   
      
    end
    
   PO(n+1,:)=A\R;
   
    for i=1:10
       if i==1
          q_t(n,i)=injrate;
          Tw=(2.6366e-4)*K*(Wmob(n,i))*(Area/0.4);
          Swi(n+1,i)= (1+(CR+CW)*(PO(n+1,i)-PO(n,i)))*Swi(n,i)+...
              (dt/(Vr*porosity))*(Tw*(PO(n+1,i+1)-PO(n+1,i))-...
              Tw*(PCOW(n,i+1)-PCOW(n,i))+q_t(n,i));
       elseif i==10
            TwL=Tw;
            Tw=0;
            WI(n,i)=(2.6366e-4)*K*(Wmob(n,i)+Omob(n,i))*(Area/(0.4/2));
            qt(n,i)= -WI(n,i)*(PO(n+1,i)-PWF);
            Swi(n+1,i)=(1+(CR+CW)*(PO(n+1,i)-PO(n,i)))*Swi(n,i)+dt/(Vr*porosity)*...
                (-TwL*(PO(n+1,i)-PO(n+1,i-1))+TwL*(PCOW(n,i)-PCOW(n,i-1))+qt(n,i));
       else 
         
           TwL=Tw;
           Tw=2.6366e-4*K*(Wmob(n,i))*(Area/0.4);
           Swi(n+1,i)=(1+(CR+CW)*(PO(n+1,i)-PO(n,i)))*Swi(n,i)+ ...
               (0.2/(Vr*porosity))*Tw*(PO(n+1,i+1)-PO(n+1,i))-...
               TwL*(PO(n+1,i)-PO(n+1,i-1))-Tw*(PCOW(n,i+1)-PCOW(n,i))+TwL*(PCOW(n,i)-PCOW(n,i-1));
       end 
       if Swi(n+1,i)<Swr
           Swi(n+1,i)=Swr;
       end 
    end 
    
    for i=1:10 
        
        if   Swi(n+1,i)<=Swr
             PCOW(n+1,i)=5;
        elseif   Swi(n+1,i)>=(1-Sorw)
             PCOW(n+1,i)=-5;
        elseif Swi(n+1,i)>Swr && Swi(n+1,i)<0.5
             PCOW(n+1,i)= min(-0.8*log((Swi(n+1,i)-Swr)/(0.5-Swr)),5);
        elseif Swi(n+1,i)>0.5 && Swi(n+1)<(1-Sorw)
             PCOW(n+1,i)=max(0.8*log((1-Swi(n+1,i)-Sorw)/(1-0.5-Sorw)),-5); 
        end 
    end
    for i=1:10
        
         Krw(n+1,i)=KrwMax*((Swi(n+1,i)-Swr)/(1-Sorw-Swr))^Nw;  
         Kro(n+1,i)=KroMax*((1-Swi(n+1,i)-Sorw)/(1-Sorw-Swr))^No;
         Omob(n+1,i)= Krw(n+1,i)/watervisc;
         Wmob(n+1,i)=Kro(n+1,i)/oilvisc;
    end 
end 

Answer 1

[Mark Sherstan]

You should not use i or j as the names of loop variables as these are both names for the inbuilt imaginary unit. Try changing that and see if the problem persisits.

[Sarah A Alruwayi]

My friend has exactly the same code and got it right and she named the loops i as well.So, I do not think this is the problem

Answer 2

[Walter Roberson]

When n = 4, i = 1, then 1-Swi(n+1),i)-Sorw becomes negative, so (1-Swi(n+1,i)-Sorw)/(1-Sorw-Swr) becomes negative. You are raising a negative value to a non-integer. A^B with scalar A, B is defined as exp(B*LOG(A)) so when A<0 the log is complex and unless B is an integer you get a complex result from the ^ operation. Everything else gets tainted from there.