Matlab help!!!

[pseudonym]

I hope some of you know Matlab, because this is starting to drive me nuts. The project is for my heat transfer class and the job is to do a program that calculates the change in temperature of a piece of metal. Now I believe that my equations are right, so thats not the problem, but when I run the program the second iteration causes many of the elements in the matrix to go from 513 (inital temp) to zero. So I end up with some funky matrix that looks like this (in much shortened form, in reality its 30x82 from this program and in the end it will be 460000x82. The first row is all 513 because that is inital temp, it should increase after the first iteration)

513 513 513 513 513 513 513
522 519 0 0 0 0 0

And so on. Every once in awhile the point in the matrix works but backwards and temperature decreases to a very low number rather than the increase I'm supposed to see. Anyways, I'm starting to think the problem is in the code because I can manually load the matrix and run the equations manually on the nodes where it zeroes out and it works fine. Can someone see a problem with this program? I believe its in the elseif statements and it must be a usage problem, but I'm not sure why it zeroes out nodes.


clear all
clc

format long

k=186;
rhoc=2907180;
deltx=.0009091;
delty=.0005;
deltt=.0008;
Tinfpre=893;
sigma=.0000000567;
epsilon=.8;

%Conduction(sidenodes)=(((2*deltt*(1.3*(Tinfpre-Tpre(m-1,n))^.33)+4*epsilon*sigma*((Tinfpre+Tpre(m-1,n))/2)^3))/(rhoc*deltx))*(Tinfpre-Tpre(m-1,n))
%Conduction(topnodes)= (((2*deltt*(1.3*(Tinfpre-Tpre(m-1,n))^.33)+4*epsilon*sigma*((Tinfpre+Tpre(m-1,n))/2)^3))/(rhoc*delty))*(Tinfpre-Tpre(m-1,n))


alpha=k/(rhoc);
taux=(deltt*alpha)/(deltx^2);
tauy=(deltt*alpha)/(delty^2);

for n=1:82
Tpre(1,n)=513;
end

for m=2:30
for n=1:82
if n==1
Tpre(m,n)=2*tauy*(Tpre(m-1,n+1)-Tpre(m-1,n))+2*taux*(Tpre(m-1,n+6)-Tpre(m-1,n))+Tpre(m-1,n);

elseif n==2,3,4,5
Tpre(m,n)=2*tauy*(Tpre(m-1,n-1)-Tpre(m-1,n))+2*tauy*(Tpre(m-1,n+1)-Tpre(m-1,n))+2*taux*(Tpre(m-1,n+6)-Tpre(m-1,n))+Tpre(m-1,n);

elseif n==6
Tpre(m,n)=(((2*deltt*(1.3*(Tinfpre-Tpre(m-1,n))^.33)+4*epsilon*sigma*((Tinfpre+Tpre(m-1,n))/2)^3))/(rhoc*delty))*(Tinfpre-Tpre(m-1,n))+2*tauy*(Tpre(m-1,n-1)-Tpre(m-1,n))+2*taux*(Tpre(m-1,n+6)-Tpre(m-1,n))+Tpre(m-1,n);

elseif n==12,18,24,30,36,42,48,54,60,66
Tpre(m,n)=(((2*deltt*(1.3*(Tinfpre-Tpre(m-1,n))^.33)+4*epsilon*sigma*((Tinfpre+Tpre(m-1,n))/2)^3))/(rhoc*delty))*(Tinfpre-Tpre(m-1,n))+2*taux*(Tpre(m-1,n-6)-Tpre(m-1,n))+2*taux*(Tpre(m-1,n+6)-Tpre(m-1,n))+2*tauy*(Tpre(m-1,n-1)-Tpre(m-1,n))+Tpre(m-1,n);

elseif n==7,13,19,25,31,37,43,49,55,61
Tpre(m,n)=2*taux*(Tpre(m-1,n-6)-Tpre(m-1,n))+2*taux*(Tpre(m-1,n+6)-Tpre(m-1,n))+2*tauy*(Tpre(m-1,n+1)-Tpre(m-1,n))+Tpre(m-1,n);

elseif n==8,9,10,11,14,15,16,17,20,21,22,23,26,27,28,29,32,33,34,35,38,39,40,41,44,45,46,47,50,51,52,53,56,57,58,59,62,63,64,65
Tpre(m,n)=(2*taux*(Tpre(m-1,n+6)-Tpre(m-1,n))+2*taux*(Tpre(m,n-6)-Tpre(m-1,n))+2*tauy*(Tpre(m-1,n+1)-Tpre(m-1,n))+2*tauy*(Tpre(m,n-1)-Tpre(m-1,n)))+Tpre(m-1,n);

elseif n==67
Tpre(m,n)=2*taux*(Tpre(m-1,n-6)-Tpre(m-1,n))+2*tauy*(Tpre(m-1,n+1)-Tpre(m-1,n))+Tpre(m-1,n);

elseif n==68,69,70,71
Tpre(m,n)=2*tauy*(Tpre(m-1,n-1)-Tpre(m-1,n))+2*tauy*(Tpre(m-1,n+1)-Tpre(m-1,n))+2*taux*(Tpre(m-1,n-6)-Tpre(m-1,n))+Tpre(m-1,n);

elseif n==72
Tpre(m,n)=(((2*deltt*(1.3*(Tinfpre-Tpre(m-1,n))^.33)+4*epsilon*sigma*((Tinfpre+Tpre(m-1,n))/2)^3))/(rhoc*deltx))*(Tinfpre-Tpre(m-1,n))+2*taux*(Tpre(m-1,n-6)-Tpre(m-1,n))+2*tauy*(Tpre(m-1,n-1)-Tpre(m-1,n))+2*tauy*(Tpre(m-1,n+1)-Tpre(m-1,n))+Tpre(m-1,n);

elseif n==73,74,75,76,77,78,79,80,81
Tpre(m,n)=(((2*deltt*(1.3*(Tinfpre-Tpre(m-1,n))^.33)+4*epsilon*sigma*((Tinfpre+Tpre(m-1,n))/2)^3))/(rhoc*deltx))*(Tinfpre-Tpre(m-1,n))+2*tauy*(Tpre(m-1,n+1)-Tpre(m-1,n))+2*tauy*(Tpre(m-1,n-1)-Tpre(m-1,n))+Tpre(m-1,n);

elseif n==82
Tpre(m,n)=(((2*deltt*(1.3*(Tinfpre-Tpre(m-1,n))^.33)+4*epsilon*sigma*((Tinfpre+Tpre(m-1,n))/2)^3))/(rhoc*deltx))*(Tinfpre-Tpre(m-1,n))+2*tauy*(Tpre(m-1,n-1)-Tpre(m-1,n))+Tpre(m-1,n);
end
end
end
end

x=1:30;
plot(x,Tpre);
axis([0 30 -100 800]);

Also, heres a pic of what the nodes do. Notice some of them actually increase like they should, but others just drop. The culprits are under the one that drops instantly to 0. Why the zero out I have no clue. Any help is appreciated!!

[TheBaron]

I can try to work through this with you, but you'll have to show me what equations youre using. since you're not using matrix multiplication at all, that can't be the issue, so one or more of your equations must be wrong

put a keyboard or a breakpoint before you enter the loop and MAKE SURE all your constants are right.

i'll look at it a little more

[TheBaron]

also as a tip, use smart indent (or do it manually)
ctrl-a,ctrl-i
makes everything easier to read (since you had an extra 'end')

[pseudonym]

All the constants are right, and I have gone through all of the equations and I don't think there are any problems with it, but I could be wrong. No matter what the equations should never zero out because they should be adding the previous value to them, which is 513.

I'm really starting to think that somehow the way i've written the else statements has caused it.

[pseudonym]

also as a tip, use smart indent (or do it manually)
ctrl-a,ctrl-i
makes everything easier to read (since you had an extra 'end')

Smart indent is on, thats just the way it ended up when it copied it into the post window. The extra end is probably from when I tried removing the elseifs and went to else if and had to add all the ends.

[TheBaron]

alright I dont like the way you have this comma seperated list in your if statements. i dont know if matlab can do this (i have never done it this way) but it seems to be printing out a bunch of garbage when its trying to evaluate those statements

[TheBaron]

yeah changing all the commas to |n== fixed it

k=186;
rhoc=2907180;
deltx=.0009091;
delty=.0005;
deltt=.0008;
Tinfpre=893;
sigma=.0000000567;
epsilon=.8;

%Conduction(sidenodes)=(((2*deltt*(1.3*(Tinfpre-Tpre(m-1,n))^.33)+4*epsilon*sigma*((Tinfpre+Tpre(m-1,n))/2)^3))/(rhoc*deltx))*(Tinfpre-Tpre(m-1,n))
%Conduction(topnodes)=(((2*deltt*(1.3*(Tinfpre-Tpre(m-1,n))^.33)+4*epsilon*sigma*((Tinfpre+Tpre(m-1,n))/2)^3))/(rhoc*delty))*(Tinfpre-Tpre(m-1,n))

alpha=k/(rhoc);
taux=(deltt*alpha)/(deltx^2);
tauy=(deltt*alpha)/(delty^2);

for n=1:82
Tpre(1,n)=513;
end

for m=2:30
for n=1:82

if n==1
Tpre(m,n)=2*tauy*(Tpre(m-1,n+1)-Tpre(m-1,n))+2*taux*(Tpre(m-1,n+6)-Tpre(m-1,n))+Tpre(m-1,n);

elseif n==2|n==3|n==4|n==5
Tpre(m,n)=2*tauy*(Tpre(m-1,n-1)-Tpre(m-1,n))+2*tauy*(Tpre(m-1,n+1)-Tpre(m-1,n))+2*taux*(Tpre(m-1,n+6)-Tpre(m-1,n))+Tpre(m-1,n);
Tpre(m,n)


elseif n==6
Tpre(m,n)=(((2*deltt*(1.3*(Tinfpre-Tpre(m-1,n))^.33)+4*epsilon*sigma*((Tinfpre+Tpre(m-1,n))/2)^3))/(rhoc*delty))*(Tinfpre-Tpre(m-1,n))+2*tauy*(Tpre(m-1,n-1)-Tpre(m-1,n))+2*taux*(Tpre(m-1,n+6)-Tpre(m-1,n))+Tpre(m-1,n);


elseif n==12|n==18|n==24|n==30|n==36|n==42|n==48|n==54|n==60|n==66
Tpre(m,n)=(((2*deltt*(1.3*(Tinfpre-Tpre(m-1,n))^.33)+4*epsilon*sigma*((Tinfpre+Tpre(m-1,n))/2)^3))/(rhoc*delty))*(Tinfpre-Tpre(m-1,n))+2*taux*(Tpre(m-1,n-6)-Tpre(m-1,n))+2*taux*(Tpre(m-1,n+6)-Tpre(m-1,n))+2*tauy*(Tpre(m-1,n-1)-Tpre(m-1,n))+Tpre(m-1,n);


elseif n==7|n==13|n==19|n==25|n==31|n==37|n==43|n==49|n==55|n==61
Tpre(m,n)=2*taux*(Tpre(m-1,n-6)-Tpre(m-1,n))+2*taux*(Tpre(m-1,n+6)-Tpre(m-1,n))+2*tauy*(Tpre(m-1,n+1)-Tpre(m-1,n))+Tpre(m-1,n);


elseif n==8|n==9|n==10|n==11|n==14|n==15|n==16|n==17|n==20|n==21|n==22|n==23|n==26|n==27|n==28|n==29|n==32|n==33|n==34|n==35|n==38|n==39|n==40|n==41|n==44|n==45|n==46|n==47|n==50|n==51|n==52|n==53|n==56|n==57|n==58|n==59|n==62|n==63|n==64|n==65
Tpre(m,n)=(2*taux*(Tpre(m-1,n+6)-Tpre(m-1,n))+2*taux*(Tpre(m,n-6)-Tpre(m-1,n))+2*tauy*(Tpre(m-1,n+1)-Tpre(m-1,n))+2*tauy*(Tpre(m,n-1)-Tpre(m-1,n)))+Tpre(m-1,n);


elseif n==67
Tpre(m,n)=2*taux*(Tpre(m-1,n-6)-Tpre(m-1,n))+2*tauy*(Tpre(m-1,n+1)-Tpre(m-1,n))+Tpre(m-1,n);


elseif n==68|n==69|n==70|n==71
Tpre(m,n)=2*tauy*(Tpre(m-1,n-1)-Tpre(m-1,n))+2*tauy*(Tpre(m-1,n+1)-Tpre(m-1,n))+2*taux*(Tpre(m-1,n-6)-Tpre(m-1,n))+Tpre(m-1,n);


elseif n==72
Tpre(m,n)=(((2*deltt*(1.3*(Tinfpre-Tpre(m-1,n))^.33)+4*epsilon*sigma*((Tinfpre+Tpre(m-1,n))/2)^3))/(rhoc*deltx))*(Tinfpre-Tpre(m-1,n))+2*taux*(Tpre(m-1,n-6)-Tpre(m-1,n))+2*tauy*(Tpre(m-1,n-1)-Tpre(m-1,n))+2*tauy*(Tpre(m-1,n+1)-Tpre(m-1,n))+Tpre(m-1,n);


elseif n==73|n==74|n==75|n==76|n==77|n==78|n==79|n==80|n==81
Tpre(m,n)=(((2*deltt*(1.3*(Tinfpre-Tpre(m-1,n))^.33)+4*epsilon*sigma*((Tinfpre+Tpre(m-1,n))/2)^3))/(rhoc*deltx))*(Tinfpre-Tpre(m-1,n))+2*tauy*(Tpre(m-1,n+1)-Tpre(m-1,n))+2*tauy*(Tpre(m-1,n-1)-Tpre(m-1,n))+Tpre(m-1,n);


elseif n==82
Tpre(m,n)=(((2*deltt*(1.3*(Tinfpre-Tpre(m-1,n))^.33)+4*epsilon*sigma*((Tinfpre+Tpre(m-1,n))/2)^3))/(rhoc*deltx))*(Tinfpre-Tpre(m-1,n))+2*tauy*(Tpre(m-1,n-1)-Tpre(m-1,n))+Tpre(m-1,n);
end

end
end

x=1:30;
plot(x,Tpre);
axis([0 30 -100 800]);

[pseudonym]

alright I dont like the way you have this comma seperated list in your if statements. i dont know if matlab can do this (i have never done it this way) but it seems to be printing out a bunch of garbage when its trying to evaluate those statements

Agreed, that sounds about right. The way I have it written should mean NOTHING gets printed till the graph, but it still does that, Matlab must not be able to handle that. The alternative is to write 82 if statements, and I obviously don't want to do that.

Heres an interesting little thing.

» Tpre(2,3)

ans =

0

» 2*tauy*(Tpre(1,2)-Tpre(1,3))+2*tauy*(Tpre(1,4)-Tpre(1,3))+2*taux*(Tpre(1,9)-Tpre(1,3))+Tpre(1,3)

ans =

513


When I type this equation in I get the answer 513 (Correct), when the prog runs the same operation it gets zero. I think you are right about the comma seperated ifs.

[pseudonym]

Dear god!!!!! Thank you so very much Baron!!!! Looks like the temp is increasing too fast, but that means the equations are a bit off, and I suspected they were off value wise but that they should function correctly. Thanks very much, I owe you!!

[mmonnin]

I hate MatLab. That is all.

[TheBaron]

thats because Matlab rawks you

[mmonnin]

Its an awesome tool, but in the class I am using it with now we are just said use Matlab to do this, this and this. We dont know the subject or Matlab and it sucks.