Heat conduction through 2D surface using Finite Difference Equation

Hello,

I struggle with Matlab and need help on a Numerical Analysis project. The assignment requires a 2D surface be divided into different sizes of equal increments in each direction, I'm asked to find temperature at each node/intersection.

The code I've developed so far for the first part is:

inodes=input('\nInput the Number of Nodes i Direction:');
jnodes=input('\nInput the Number of Nodes j Direction:');

%Rectangular Flat Plate
w=2.5;
h=2.0;
%Discretize the Plate Area
deli=w/(inodes-1);
delj=h/(jnodes-1);
if deli~=delj
  disp('inappropriate node dimmensions');
  return
else
  disp('Assume Delta i = Delta j');
  dx=deli;
end 

xrange=0:deli:w;
yrange=0:delj:h;

% generate appropriate x,y matrices
[x y] = meshgrid(xrange,yrange);
% boundary conditions

for t=i*j
    t(0,:)=0;
    t(2.5,:)=0;
    t(:,0)=0;
    t(:,2.0)=0;
    for t=~0
         t(x,y)=(t(x+1,y)+t(x-1,y)+t(x,y+1)+t(x,y-1))/4;
    end
end

I'm having trouble assigning the temperature reference to each node. Once this is complete I will use a matrix method with corresponding boundary conditions.

I'm also required to calculate an exact solution for comparison. My exact solution code so far is as follows.

%Plate dimensions

%Exact Equation
for y=0:2
    w=2.5;
    h=2.0;
    Qdot=400.0;  %J/(cm^3*s)
    k=.4;  %J/(cm*s*C)
    for x=0:2.5
    for n=0:201
     Ta=(((-1^n).*cos((2.*n+1).*(x-w/2).*(pi()./w)).*cosh((2.*n+1).*(y-h/2).*pi()./w)))/(((2.*n+1)^3).*cosh((2.*n+1).*pi().*h/(2.*w)));
     T1(:,n+1)=Ta
    end
    Tc=(Qdot/k)*((w.^2)/4)-(x-(w/2)^2);
    Tb=-((4*w^2*(Qdot/k))/pi()^3)
    Tval(x)=Tc+sum(T1(x,:))*Tb;
    end
    T(y+1,:)=Tval
end

Currently this code doesn't work, causing a matlab error that says "the cell can not be evaluated be cause it contains an incomplete statement." I'm unable to find the error. Also, I'm not sure that my code will accurately find a temperature at each node.

Ideally, I will construct 3D plots for each code for comparison.

Answer 1

[Walter Roberson]

Have another look at this loop:

for t=i*j
    t(0,:)=0;
    t(2.5,:)=0;
    t(:,0)=0;
    t(:,2.0)=0;
    for t=~0
         t(x,y)=(t(x+1,y)+t(x-1,y)+t(x,y+1)+t(x,y-1))/4;
    end
end

You have not assigned anything to i or j, so both values will have their default meaning of sqrt(-1) . sqrt(-1)*sqrt(-1) is -1, so your outer loop is "for t=-1" which means that the loop is to be executed once with the value of t set to the scalar value -1 . Inside the loop, you attempt to access t(0,:) . As t is a scalar, it's last dimension is 1, so that will be the same as trying to access t(0,1) which will fail because arrays cannot be subscripted at 0.

If it did not fail there, it would fail at the next line, t(2.5,:) as it is not allowed to access arrays at fractional indices.

If you managed to get through the next several lines, you would hit your "for t=~0" . ~0 is the scalar containing the logical value "true", so you would be asking for the loop to be executed once with t temporarily assigned true . And then you would try to access that 1x1 scalar at all kinds of non-integral points defined by x and y...

Arrays must be indexed with non-negative integers. The array position might logically correspond to some x value such as 2.5 but you need to maintain that list of correspondences independently. If you have a point in logical coordinates that you want to convert to array coordinates, histc() or mod() or interp1() can be useful.

[Matt Fig]

Walter, that line cracks me up (quoting still not fixed, grrrr):

Arrays must be indeed with non-negative integers.

Indeed, some arrays must be with non-negative integers! Some must be with negative integers or even zeros! (Of course I know it is a typo!)

[Matt Fig]

@Travis,

In case Walter doesn't fix that line, he meant to say that arrays in MATLAB must be indexed with integers greater than 0.

[Walter Roberson]

Blush ;-)

[Travis]

Thanks for the quick response Walter and Matt. I'm attempting to correct my mistakes but as I said before, I struggle with Matlab. Perhaps I should have said I struggle with code writing. I'm not sure how to correct my array even with the knowledge of Matlabs index requirement... I changed my i*j reference to be inodes*jnodes, but my nested (i think that's the name on the loop within a loop) for loop is still hung up on cases when t>0.

[Walter Roberson]

for idxy = 2:size(t,2)-1

for idxx = 2:size(t,1)-1

if t(x,y) ~= 0

t(x,y)=(t(x+1,y)+t(x-1,y)+t(x,y+1)+t(x,y-1))/4;

end

end

end

However, I have to question this. No matter whether you sweep rows or columns, no matter whether you sweep top to bottom or bottom to top, this code will be accessing t values that have already been updated along with t values that have not yet been updated. You should be consistent, such as with:

newt = t;

for idxy = 2:size(t,2)-1

for idxx = 2:size(t,1)-1

if t(x,y) ~= 0

newt(x,y)=(t(x+1,y)+t(x-1,y)+t(x,y+1)+t(x,y-1))/4;

end

end

end

t = newt;

Note: The code can be vectorized to at least some degree, but you need to get the code working first before you think about vectorizing it.

[Travis]

Thanks Walter, I'm going to work some more on it tonight and I will post back about my progress and struggle.

[Walter Roberson]

Correction to my earlier comment:

newt = t;

for idxy = 2:size(t,2)-1

for idxx = 2:size(t,1)-1

if t(idxx,idxy) ~= 0

newt(idxx,idxy)=(t(idxx+1,idxy)+t(idxx-1,idxy)+t(idxx,idxy+1)+t(idxx,idxy-1))/4;

end

end

end

t = newt;

Answer 2

[PAPA MAMA]

This is my project topic: Consider the two dimensional heat conduction equation, δ2φ/δx2 + δ2φ/δy2 = δφ/δt 0≤ x,y ≤2; t>0 subject to the boundary condition φ(x,y,t) =0, on the boundary, for t>0 and the initial conditions φ(x,y,0)= cos(

[Walter Roberson]

http://www.mathworks.com/matlabcentral/answers/6200-tutorial-how-to-ask-a-question-on-answers-and-get-a-fast-answer

[PAPA MAMA]

only what am asking its about the challeges i was facing to run my home work in matlab

T(i,j,m+1)=(1-4*r)*T(i,j,m)*(T(i+1,j,m)+T(i-1,j,m)+T(i,j+1,m)+T(i,j-1,m));

this is major area of comment of hypothesis or bracket either

thanls

[Walter Roberson]

newt = T(:,:.m);
for i = 2:size(T,2)-1
  for j = 2:size(T,1)-1
    newt(idxx,idxy) = (1-4*r)*T(i,j,m)*(T(i+1,j,m)+T(i-1,j,m)+T(i,j+1,m)+T(i,j-1,m));
  end
end
T(:,:,m+1) = newt;