| Prahallada |
August 14, 2009 12:05 |
Need Help in Debugging Matlab Code
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.
Code:
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 |
