Stuck at a For Loop

Hi, good day! I am having a problem with my coding.

I want to develop a code for Neumann boundary condition, where one edge of a heat plate is insulated, using Liebmann (or Jacobi?) Method.

However, my final answer returns the original matrix T (representing the temperatures), which means my iteration failed, but I don't know why...

Appreciate much for any guidance! Thanks!

% Overrelaxation, lambda
lam = input('Please enter the overrelaxation value, lambda:   ');
e_s = input('Please enter the stopping criterion in %:   ');
err = e_s/100;
disp('The unit of k value is suggested to be in cal/s.cm.Celcius')
disp('As long as the unit for dimensions is consistent')
k = input('Please enter the k value of the plate material:   ');
disp('The dimensions are suggested to be in cm')
l = input('Please enter the length of the heated plate:   ');
w = input('Please enter the width of the heated plate:   ');
pt = input('Please enter the desired number of points per row within the plate:   ');

% Calculating the temperatures
delx = w/(pt + 1);
dely = l/(pt + 1);
T = zeros(pt + 2);

T(:,1) = 75;
T(:, (pt + 2)) = 50;
T(pt + 2, (2:pt + 1)) = 100;

disp('T_up + T_down + T_left + T_right = 4*T_current')
disp('At insulated edge, T_up is not available, use fictitious point')
disp('Hence the equation becomes:')
disp('2*T_down + 2*dely*dTdy + T_left + T_right = 4*T_current')
disp('Since dTdy = 0 at insulated edge, the equation becomes:')
disp('T_current = (2*T_down + T_left + T_right)/4')

A = T;
e1 = ones((pt + 1), pt);

while max(max(e1)) > err
    % At insulated edge
    for r = 1
        for c = 2:(pt + 1)
            A(r,c) = (T(r,(c-1)) + T(r, (c+1)) + (2*T((r+1),c)))/4;
        end
    end
    
    % At other points
    for r = 2:(pt + 1)
        for c = 2:(pt + 1)
            A(r,c) = (T(r,(c-1)) + T(r,(c+1)) + T((r-1),c) + T((r+1),c))/4;
        end
    end
    
    T(r,c) = (lam*A(r,c)) + ((1-lam)*T(r,c));
    e1(r,(c-1)) = abs((T(r,c) - A(r,c))/T(r,c));
    
end

disp('The temperature distribution on the plate is:')
T_plate = flipud(T);
disp(T_plate)

[Walter Roberson]

What input values are you using, so we can test the code?

[LEW JUN KIT]

Hi, Mr. Roberson! These are my inputs.

lam = 1.5
e_s = 1
k = 0.49
l = w = 40
pt = 3

Answer 1

[per isakson]

Notes from my debugging session:

• added your assignments of input values to the top of the script
• commented out the input() statements
• ran the script
• the while-loop looped forever
• the following two line seems to be in the wrong place

    T(r,c) = (lam*A(r,c)) + ((1-lam)*T(r,c));
    e1(r,(c-1)) = abs((T(r,c) - A(r,c))/T(r,c));

• using a loop counter, in this case r and c, outside the loop is allowed in Matlab, but it smells
• the values of r and c, don't change. ( pt doesn't change value inside the while-loop.) At these two lines r and c, holds the end-values of the previous for-loops, respectively.
• thus only the lower right values of T and e are changed, which explains why the while condition never takes the value false

Your turn

P.S. See: Debug a MATLAB Program and Examine Values While Debugging