[Prahallada]

Hi guys...Im trying to formulate a mathematical model, which was already proposed in a paper using matlab for my B.Tech final year project work.I finished most of the coding part and when i try to run the code i get an error,which i'm unaware of.

clear
%chan = ddeinit('excel', 'Sheet1');
R=.02%input('enter the value of  inside radius :');
Ro=.03;
l=0.5%input('enter the value of length of heat source:');
%input('enter the length of pipe:');
l1=0.1%input('enter lengh between two heat source:');
Te=20%input('enter the value of free stream temperature:');
q=1000%input('enter the value of heat flux:');
U=5%input('enter free stream velocty:');
k=12%input('enter the value of thermal conductivity ofi inside fluid:');
kw=12;
alpha=.01%input('enter the value of thermal difusivity of inside fluid:');
alpha1=4*alpha;
a=input('enter no of grid in radial directon must be odd no:');
b=input('enter no of grid in axial direction :');
n=input('enter no of heat sources on periphery:');
er=input('enter value of max possible percentage error:');
s=input('enter relaxatation factor around 1.05:')
Pe=U*2*R/alpha
deln1=(Ro/R)/(a-1);
L=n*(l+l1);
c=L/(R*(b-1));
K=kw/k;
A1=alpha/alpha1;
h=0;
delt(1,1)=0;
Nu=0;
g=((a-1)*(R/Ro)+1)
warning off all
  for i=1:a
             for j=1:b
                
                % T(i,j)=Te;
                 %t(i,j)=(T(i,j)-Te)/(q*R/k);
                 t(i,j)=0.1;
             end
   end
 format long g
 erfmx=10;
 erf=1;
 erf1=0;
 itn=1;
while(er<erfmx)
     erfmx=0;
    for i=1:a
     for j=1:b
          if j==1
          t(i,1)=0;
          end
           tn(i,b)=0;
        if j>1 && j<b
            if i==1
                X=Pe*c + 1;
                Y=Pe*c + 2;
                tn(i,j)=(t(i,j+1)+X*t(i,j-1))/Y;
            end
            if i>1 && i<((a-1)*(R/Ro)+1)
                  A=(Pe*(1-(((i-1)*Ro)/((a-1)*R))^2))/c + 2*((a-1)*(R/Ro))^2 + 2/c^2 + (((a-1)*(R/Ro))^2)/(i-1);
                  B=((a-1)*(R/Ro))^2  - (((a-1)*(R/Ro))^2)/(2*(i-1));
                  C=((a-1)*(R/Ro))^2 + (((a-1)*(R/Ro))^2)/(2*(i-1));
                  D=(Pe*(1-(((i-1)*Ro)/((a-1)*R))^2))/c + 1/c^2;
                  E=1/c^2;
             tn(i,j)=(B*t(i-1,j)+C*t(i+1,j)+D*t(i,j-1)+E*t(i,j+1))/A;
            end
            if i==((a-1)*(R/Ro)+1)
                A=(Pe*(1-(((i-1)*Ro)/((a-1)*R))^2))/c + 2*((a-1)*(R/Ro))^2 + 2/c^2 + (((a-1)*(R/Ro))^2)/(i-1);
                B=((a-1)*(R/Ro))^2  - (((a-1)*(R/Ro))^2)/(2*(i-1));
                C=((a-1)*(R/Ro))^2 + (((a-1)*(R/Ro))^2)/(2*(i-1));
                D=(Pe*(1-(((i-1)*Ro)/((a-1)*R))^2))/c + 1/c^2;
                E=1/c^2;
                x=(L/(b-1))*(j-1);
                o=3;
                if   x>(o-1)*(l+l1) && x<=(o*l+(o-1)*l1)
                 tn(i,j)=(B*t(i-1,j)+C*(t(i-1,j)+2/(a-1))+D*t(i,j-1)+E*t(i,j+1))/A; 
                end
                 if  x<=(o-1)*(l+l1)  &&  x>(o*l+(o-1)*l1) 
                   tn(i,j)=((B+C)*t(i-1,j)+D*t(i,j-1)+E*t(i,j+1))/A;
                 end  
            end 
              if i>((a-1)*(R/Ro)+1) && i<a
                  A=(Pe*(1-(((i-1)*Ro)/((a-1)*R))^2))/c + 2*((a-1)*(R/Ro))^2 + 2/c^2 + (((a-1)*(R/Ro))^2)/(i-1);
                  B=((a-1)*(R/Ro))^2  - (((a-1)*(R/Ro))^2)/(2*(i-1));
                  C=((a-1)*(R/Ro))^2 + (((a-1)*(R/Ro))^2)/(2*(i-1));
                  D=(Pe*(1-(((i-1)*Ro)/((a-1)*R))^2))/c + 1/c^2;
                  E=1/c^2;
                  tn(i,j)=(B*t(i-1,j)+C*t(i+1,j)+E*(t(i,j-1)+t(i,j+1)))/(2*((a-1)*(R/Ro))^2 + 2/c^2);
              end
                  if i==a
                    A=(Pe*(1-(((i-1)*Ro)/((a-1)*R))^2))/c + 2*((a-1)*(R/Ro))^2 + 2/c^2 + (((a-1)*(R/Ro))^2)/(i-1);
                    B=((a-1)*(R/Ro))^2  - (((a-1)*(R/Ro))^2)/(2*(i-1));
                    C=((a-1)*(R/Ro))^2 + (((a-1)*(R/Ro))^2)/(2*(i-1));
                    D=(Pe*(1-(((i-1)*Ro)/((a-1)*R))^2))/c + 1/c^2;
                    E=1/c^2;
                    x=(L/(b-1))*(j-1);
                    o=3;
                      if   x>(o-1)*(l+l1) && x<=(o*l+(o-1)*l1)
                          tn(i,j)=(B*t(i-1,j)+C*(Ro/(R*(a-1))+t(i-1,j))+E*(t(i,j-1)+t(i,j+1)))/(2*((a-1)*(R/Ro))^2 + 2/c^2);
                      end
                      if  x<=(o-1)*(l+l1)  &&  x>(o*l+(o-1)*l1) 
                      tn(i,j)=((B+C)*t(i-1,j)+E*(t(i,j-1)+t(i,j+1)))/(2*((a-1)*(R/Ro))^2 + 2/c^2);
                      end  
                  end
        end
        if j==b
             if i==1
             tn(i,j)=t(1,j-1);
             end
             if i>1 && i<((a-1)*(R/Ro)+1)
                  A=(Pe*(1-(((i-1)*Ro)/((a-1)*R))^2))/c + 2*((a-1)*(R/Ro))^2 + 2/c^2 + (((a-1)*(R/Ro))^2)/(i-1);
                  B=((a-1)*(R/Ro))^2  - (((a-1)*(R/Ro))^2)/(2*(i-1));
                  C=((a-1)*(R/Ro))^2 + (((a-1)*(R/Ro))^2)/(2*(i-1));
                  D=(Pe*(1-(((i-1)*Ro)/((a-1)*R))^2))/c + 1/c^2;
                  E=1/c^2;
                  tn(i,j)=(B*t(i-1,j)+(D+E)*t(i,j-1)+ C*t(i+1,j))/A;
             end
             if i==((a-1)*(R/Ro)+1)
                  A=(Pe*(1-(((i-1)*Ro)/((a-1)*R))^2))/c + 2*((a-1)*(R/Ro))^2 + 2/c^2 + (((a-1)*(R/Ro))^2)/(i-1);
                  B=((a-1)*(R/Ro))^2  - (((a-1)*(R/Ro))^2)/(2*(i-1));
                  C=((a-1)*(R/Ro))^2 + (((a-1)*(R/Ro))^2)/(2*(i-1));
                  D=(Pe*(1-(((i-1)*Ro)/((a-1)*R))^2))/c + 1/c^2;
                  E=1/c^2;
                  tn(i,j)=((B+C)*t(i-1,j)+(D+E)*t(i,j-1))/A;
             end
              if i>((a-1)*(R/Ro)+1) && i<a
                  A=(Pe*(1-(((i-1)*Ro)/((a-1)*R))^2))/c + 2*((a-1)*(R/Ro))^2 + 2/c^2 + (((a-1)*(R/Ro))^2)/(i-1);
                  B=((a-1)*(R/Ro))^2  - (((a-1)*(R/Ro))^2)/(2*(i-1));
                  C=((a-1)*(R/Ro))^2 + (((a-1)*(R/Ro))^2)/(2*(i-1));
                  D=(Pe*(1-(((i-1)*Ro)/((a-1)*R))^2))/c + 1/c^2;
                  E=1/c^2;
                  tn(i,j)=(B*t(i-1,j)+C*t(i+1,j)+2*E*t(i,j-1))/(2*((a-1)*(R/Ro))^2 + 2/c^2);
              end
                  if i==a
                    A=(Pe*(1-(((i-1)*Ro)/((a-1)*R))^2))/c + 2*((a-1)*(R/Ro))^2 + 2/c^2 + (((a-1)*(R/Ro))^2)/(i-1);
                    B=((a-1)*(R/Ro))^2  - (((a-1)*(R/Ro))^2)/(2*(i-1));
                    C=((a-1)*(R/Ro))^2 + (((a-1)*(R/Ro))^2)/(2*(i-1));
                    D=(Pe*(1-(((i-1)*Ro)/((a-1)*R))^2))/c + 1/c^2;
                    E=1/c^2;
                    tn(i,j)=((B+C)*t(i-1,j)+2*E*t(i,j-1))/(2*((a-1)*(R/Ro))^2 + 2/c^2);
                   
                  end
        end
        erf=abs(((tn(i,j)-t(i,j))/(tn(i,j)+t(i,j)))*200);
                t(i,j)=(1-s)*t(i,j)+s*tn(i,j);
                     if erf>erfmx
                        erfmx=erf;
                     end   
      itn=itn+1;   
     end
    end
        
        itn=itn+1;
end
itn;
t';
erfmx;
t;
T= t*(q*R/k)+Te;
 itn 
 h=(Ro/R)/(a-1);
 v(1)=0;
f(1)=0;
s=f(1);
for j=1:b
for i=2:2:g-2
    v(i)=(i-1)*(Ro/R)/(a-1);
    f(i)=(v(i)-v(i)^3)*t(i,j);
    s=s+4*f(i);
    v(i+1)=(i)*(Ro/R)/(a-1);
    f(i+1)=(v(i+1)-v(i+1)^3)*t(i+1,j);
    s=s+2*f(i+1);
end
i=g-1;
v(g-1)=(g-2)*h;
p=(g-2)*h;
s=s+4*(v(g-1)-(v(g-1)^3))*t(i,j);
s=s;
sum=1*s/(3*(g-1));
i=g;
tb(i,j)=4*sum;
delt(i,j)=t(g,j)-tb(i,j);
Nu(i,j)=0;
x=(L/(b-1))*(j-1);
             o=3;
    if x>(o-1)*(l+l1)&& x<= (o*l+(o-1)*l1)
        Q(i,j)=1;
        Nu(i,j)=2*(Ro/R)*Q(i,j)/delt(i,j);
    end
end
Nu;
tb;
delt;

fid=fopen('D:\matlabresult\hs1.txt','wt')
          fprintf(fid,'radius    nonrad   inc   axis      ndaxis      incx     temp       nondt     Nu     delt     tb\n');
for i=1:a
       n1(i)=(i-1)/(a-1);
       r(i)=n1(i)*R;
   
 for j=1:b
        x(j)= (L/(b-1))*(j-1);
          xi(j)=x(j)/R;
         rt=[r(i);n1(i);i;x(j);xi(j);j;T(i,j);t(i,j);Nu(i,j);delt(i,j);tb(i,j)];
        % csvwrite('rk.txt',rt); 
          %[c,h]=contour(x,r,T);
          %set(h,'ShowText','on','TextStep',get(h,'LevelStep')*2)
          %colormap cool
          %o=((i-1)*b)+j;
          %xlabel('x')
          %ylabel('r')
         %clabel(cs)
         
         fprintf(fid,'%1.4f   %2.4f  %3d    %1.4f     %2.4f     %3d       %6.4f    %1.4f    %1.2f    %1.2f    %1.2f\n',rt);
       
      
        
  end
 end
fclose(fid);

I'm getting an error in the line which is bolded and coloured red in the code.The error states "Subscript indices must either be real positive integers or logicals".

When i use a value of 0.02 for Ro the code is working fine without any errors.The error starts popping only for Ro=0.03.
Can some one please help me out in debugging this code.I would be thankful if some one can do so

Thanking You
Prahallada J

[f-w]

v(g-1)=(g-2)*h

In the line above, (g-1) must be a real, positive integer (it can't be 1.222 or 4.11).

Your formulations for g, R, Ro, a:
a= an integer
R=.02
Ro=.03;
g=((a-1)*(R/Ro)+1)

So, when Ro=0.02, "g" works out to be an integer, otherwise, it's not. So, make it a multiple of R.

[Prahallada]

Thanks a lot dude.....