MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn moreOpportunities for recent engineering grads.

Apply Today**New to MATLAB?**

Asked by Charles
on 27 Mar 2012

1. The problem statement, all variables and given/known data

Having trouble with code as seen by the gaps left where it asks me to add things like the coefficient matrices. Any Numerical Conduction matlab wizzes out there?

A long square bar with cross-sectional dimensions of 30 mm x 30 mm has a specied temperature on each side, The temperatures are:

Tbottom = 100 C Ttop = 150 C Tleft = 250 C Tright = 300 C

Assuming isothermal surfaces, write a software program to solve the heat equation to determine the two-dimensional steady-state spatial temperature distribution within the bar. Your analysis should use a finite difference discretization of the heat equation in the bar to establish a system of equations:

2. Relevant equations

AT = C

Must use Gauss-Seidel Method to solve the system of equations

3. The attempt at a solution

clear all close all %Specify grid size Nx = 10; Ny = 10; %Specify boundary conditions Tbottom = 100 Ttop = 150 Tleft = 250 Tright = 300 % initialize coefficient matrix and constant vector with zeros A = zeros(Nx*Ny); C = zeros(Nx*Ny,1); % initial 'guess' for temperature distribution T(1:Nx*Ny,1) = 100; % Build coefficient matrix and constant vector % inner nodes for n = 2:(Ny-1) for m = 2:(Nx-1) i = (n-1)*Nx + m; A(i,i+Nx) = 1; A(i,i-Nx) = 1; A(i,i+1) = 1; A(i,i-1) = 1; A(i,i) = -4; end end % Edge nodes % bottom for m = 2:(Nx-1) %n = 1 i = m;

A(i,i+Nx) = 1; A(i,i+1) = 1; A(i,i-1) = 1; A(i,i) = -4; C(i) = -Tbottom;

end %top: for m = 2:(Nx-1) % n = Ny i = (Ny-1)*Nx + m;

A(i,i-Nx) = 1; A(i,i+1) = .5; A(i,i-1) = .5; A(i,i) = -5; C(i) = -Ttop;

end %left: for n=2:(Ny-1) %m = 1 i = (n-1)*Nx + 1; A(i,i+Nx) = .5; A(i,i+1) = 1; A(i,i-Nx) = .5; A(i,i) = -2; end %right: for n=2:(Ny-1) %m = Nx i = (n-1)*Nx + Nx; A(i,i+Nx) = 1; A(i,i-1) = 1; A(i,i-Nx) = 1; A(i,i) = -4; C(i) = -Tright;

DEFINE COEFFICIENT MATRIX AND CONSTANT VECTOR ELEMENTS HERE end % Corners %bottom left (i=1): i=1 A(i,Nx+i) = 1; A(i,2) = 1; A(i,1) = -4; C(i) = -(Tbottom + Tleft); %bottom right: i = Nx A(i,Nx+i) = 1; A(i,2) = 1; A(i,1) = -4; C(Nx) = -(Tbottom + Tright); %top left: i = (Ny-1)*Nx + 1; A(i,i+1) = .5 A(i,i) = %top right: i = Nx*Ny; DEFINE COEFFICIENT MATRIX AND CONSTANT VECTOR ELEMENTS HERE %Solve using Gauss-Seidel residual = 100; iterations = 0; while (residual > 0.0001) % The residual criterion is 0.0001 in this example 7% You can test different values iterations = iterations+1 %Transfer the previously computed temperatures to an array Told Told = T; %Update estimate of the temperature distribution INSERT GAUSS-SEIDEL ITERATION HERE %compute residual deltaT = abs(T - Told); residual = max(deltaT); end iterations % report the number of iterations that were executed %Now transform T into 2-D network so it can be plotted. delta_x = 0.03/(Nx+1) delta_y = 0.03/(Ny+1) for n=1:Ny for m=1:Nx i = (n-1)*Nx + m; T2d(m,n) = T(i); x(m) = m*delta_x; y(n) = n*delta_y; end end T2d surf(x,y,T2d) figure contour(x,y,T2d)

*No products are associated with this question.*

Answer by Charles
on 28 Mar 2012

clear all close all %Specify grid size Nx = 10; Ny = 10; %Specify boundary conditions Tbottom = 100 Ttop = 150 Tleft = 250 Tright = 300 % initialize coefficient matrix and constant vector with zeros A = zeros(Nx*Ny); C = zeros(Nx*Ny,1); % initial 'guess' for temperature distribution T(1:Nx*Ny,1) = 100; % Build coefficient matrix and constant vector % inner nodes for n = 2:(Ny-1) for m = 2:(Nx-1) i = (n-1)*Nx + m; A(i,i+Nx) = 1; A(i,i-Nx) = 1; A(i,i+1) = 1; A(i,i-1) = 1; A(i,i) = -4; end end % Edge nodes % bottom for m = 2:(Nx-1) %n = 1 i = m;

A(i,i+Nx) = 1; A(i,i+1) = 1; A(i,i-1) = 1; A(i,i) = -4; C(i) = -Tbottom;

end %top: for m = 2:(Nx-1) % n = Ny i = (Ny-1)*Nx + m;

A(i,i-Nx) = 1; A(i,i+1) = 1; A(i,i-1) = 1; A(i,i) = -4; C(i) = -Ttop;

end %left: for n=2:(Ny-1) %m = 1 i = (n-1)*Nx + 1; A(i,i+Nx) = 1; A(i,i+1) = 1; A(i,i-Nx) = 1; A(i,i) = -4; end %right: for n=2:(Ny-1) %m = Nx i = (n-1)*Nx + Nx; A(i,i+Nx) = 1; A(i,i-1) = 1; A(i,i-Nx) = 1; A(i,i) = -4; C(i) = -Tright;

end % Corners %bottom left (i=1): i=1; A(i,Nx+i) = 1; A(i,2) = 1; A(i,1) = -4; C(i) = -(Tbottom + Tleft); %bottom right: i = Nx; A(i,i+Nx) = 1; A(i,i-1) = 1; A(i,i) = -4; C(i) = -(Tbottom + Tright); %top left: i = (Ny-1)*Nx + 1; A(i,i+1) = 1; A(i,i) = -4; A(i,i-Nx) = 1; C(i) = -(Ttop + Tleft); %top right: i = Nx*Ny; A(i,i-1) = 1; A(i,i) = -4; A(i,i-Nx) = 1; C(i) = -(Tright + Ttop); %Solve using Gauss-Seidel residual = 100; iterations = 0; while (residual > 0.0001) % The residual criterion is 0.0001 in this example % You can test different values iterations = iterations+1; %Transfer the previously computed temperatures to an array Told Told = T; %Update estimate of the temperature distribution

for n=1:Ny for m=1:Nx i = (n-1)*Nx + m; Told(i) = T(i); end end % iterate through all of the equations for n=1:Ny for m=1:Nx i = (n-1)*Nx + m; %sum the terms based on updated temperatures sum1 = 0; for j=1:i-1 sum1 = sum1 + A(i,j)*T(j); end %sum the terms based on temperatures not yet updated sum2 = 0; for j=i+1:Nx*Ny sum2 = sum2 + A(i,j)*Told(j); end % update the temperature for the current node T(i) = (1/A(i,i)) * (C(i) - sum1 - sum2); end end residual = max(T(i) - Told(i)); end

%compute residual deltaT = abs(T - Told); residual = max(deltaT);

iterations; % report the number of iterations that were executed %Now transform T into 2-D network so it can be plotted. delta_x = 0.03/(Nx+1); delta_y = 0.03/(Ny+1); for n=1:Ny for m=1:Nx i = (n-1)*Nx + m; T2d(m,n) = T(i); x(m) = m*delta_x; y(n) = n*delta_y; end end T2d; surf(x,y,T2d) figure contour(x,y,T2d)

Answer by Ahmed Hakim
on 17 Nov 2012

Very nice Code

I would like to use SOR method for finding the optimum omega...can u help me?

thanks

## 7 Comments

## Walter Roberson

Direct link to this comment:http://www.mathworks.com/matlabcentral/answers/33651#comment_70538

Please be more specific about the difference between what you observe and what you expect.

## Charles

Direct link to this comment:http://www.mathworks.com/matlabcentral/answers/33651#comment_70547

I am not sure what to expect, as this is my first time working with matlab. I just edited the code adding in what I think is the correct coefficient matrices, do the values look right?

## Walter Roberson

Direct link to this comment:http://www.mathworks.com/matlabcentral/answers/33651#comment_70548

I am not familiar with heat transfer equations, so I have no idea if the coefficients "look right".

I _am_ familiar with many MATLAB error messages, and with tracing back _specific_ problems to causes (given inputs), but there isn't much I can help with unless you can clarify

Having trouble with code as seen by the gaps left where it asks me to add things like the coefficient matrices"

## Geoff

Direct link to this comment:http://www.mathworks.com/matlabcentral/answers/33651#comment_70549

Most people here are probably not involved in your course, so we don't necessarily know what the correct values are. We are here to help with MatLab questions or issues, not to verify your results. I expect your results can be used to reproduce the initial conditions. That is verification. May I also suggest (if you doubt the correctness of your code) that you set breakpoints at relevant places (F12), run it, and check your intermediate values for sanity. Use F10 to step through code one line at a time, F5 to resume, Shift-F5 to cancel....

## Charles

Direct link to this comment:http://www.mathworks.com/matlabcentral/answers/33651#comment_70552

hold on i think im onto something. be right back

## Charles

Direct link to this comment:http://www.mathworks.com/matlabcentral/answers/33651#comment_70566

okay I cranked it out below, its working, but the distributions of the heat looks kind of odd. It is so overwhelming hotter at the top. Idk if anyone knows how to check for wrong coefficient matrice values or something, please give me a tip :)

## Image Analyst

Direct link to this comment:http://www.mathworks.com/matlabcentral/answers/33651#comment_70570

You mean like how to use the debugger? Like how to stop at breakpoints and how you can "check for wrong coefficient matrice values or something"? You can click in the margin to set a breakpoint and, after you've stopped at it, click in a variable and type control-D to check it's value in the variable editor, or look in the workspace panel.