Ex 1: Steady state temperature profile in a cylindrical wall

The equation for steady state diffusion of heat through a cylindrical wall is given by

$\frac{d^2 T}{dr^2} + \frac{1}{r} \frac{dT}{dr} = 0, T(1) = 1, T'(2) = -1$

Analytic solution

The analytic solution to this differential equation is given by

$T(r) = 1 - 2\log(r)$

Discretised Equation and Boundary Conditions

Using Matlab indexing the boundary condition $T(1) = 1$ is discretised as

$$T_1 = 1$$

The other boundary condition $T'(2) = -1$ is slightly more complicated. Using the second order backwards difference the equation can be discretised with hr = 0.25 as

$$2T_3 - 8T_4 + 6T_5 = -1$$

The inner points are simply discretised with a central difference as

$16T_{i-1} - 32T_i + 16T_{i+1} + \frac{8}{i+3}(T_{i+1} - T_{i-1}) = 0$

Where r is given by $\frac{4}{i+3}$

Matrix Form

The matrix form is dependent on a changing $\frac{1}{r}$ term which requires some calculation while building the matrix. The final matrix problem is $\left[ \begin{array}{ccccc} 1 & 0 & 0 & 0 & 0 \\ 14.4 & -32 & 17.4 & 0 & 0 \\ 0 & 14.67 & -32 & 17.33 & 0 \\ 0 & 0 & 14.86 & -32 & 17.4 \\ 0 & 0 & 2 & -8 & 6 \end{array} \right] \left[ \begin{array}{c} T_1\\ T_2\\ T_3\\ T_4\\ T_5 \end{array} \right] = \left[ \begin{array}{c} 1\\ 0\\ 0\\ 0\\ -1 \end{array} \right]$

Matrix and Numerical Solution

Inversion of the matrix allows the T values to easily be obtained.

A = [1 0 0 0 0;
    14.4 -32 17.4 0 0;
    0 14.67 -32 17.33 0;
    0 0 14.86 -32 17.14;
    0 0 2 -8 6]
b = [1; 0; 0; 0; -1];

t = A\b        %Use matrix inversion to find the T values

plot(1:0.25:2,t,'-xr');     %plot the values over the proper range for r
hold on;

ezplot('1 - 2*log(x)',[1,2])      %plot the analytic solution over the same range
title('Solutions to heat diffusion in a cylindrical wall');
xlabel('r');
ylabel('temperature');
legend('Matrix method','analytic solution');
hold off;

% The matrix method is suprisingly accurate given the small number
% of panels used in the solution. This is probably at least partially due
% to using 2nd order finite difference approximations throughout.

A =

    1.0000         0         0         0         0
   14.4000  -32.0000   17.4000         0         0
         0   14.6700  -32.0000   17.3300         0
         0         0   14.8600  -32.0000   17.1400
         0         0    2.0000   -8.0000    6.0000


t =

    1.0000
    0.5466
    0.1776
   -0.1347
   -0.4054

Exercise 2

The equation for diffusion of temperature from a metal plate, taking into account the transfer of heat from each surface is given by

$\nabla^2 T(x,y) - \beta(T(x,y)-T_a)=0$

where in the following problem we have $\beta = 1m^{-2}$ and $T_a = 25^{\circ}C$ and the hot plate is 0.75m by 1m.

For interior points the second order central difference approximation is given by

$\frac{T_{i(j-1)} + T_{i(j+1)} + T_{(i-1)j} + T_{(i+1)j} -4T_{ij}}{0.25^2} - (T_{ij} - T_a) = 0$

with the boundary condition along the short edges of the form

$T_{i1} = 0$

and the boundary condition along the longer edges of the form

$\frac{T_{1j} - 2T_{2j} + T_{3j}}{h^2} = 0$

where the above equation specifies one of the edges as an insulator.

We can the proceed to sum the various components for our 20x20 matrix with rows in the form $T_{11}, T_{12},...,T_{15},T_{21},...,T_{55}

clear;

A = zeros(20,20);         %pre allocate the size of the matrix
b = zeros(20,1);            %size of the results vector

for i=1:4
    A(5*i-4,5*i-4) = 1;     %The left side is set to 0 degrees
    b(5*i-4) = 0;
    A(5*i,5*i) = 1;         %The right side is 100 degrees
    b(5*i) = 100;
end

for j=2:4
    A(j,j) = A(j,j) - 6;                %the top side is an insulator
    A(j,j+5) = A(j,j+5) + 8;
    A(j,j+10) = A(j,j+10) - 2;
    b(j) = 0;
    A(j+15,j+15) = A(j+15,j+15) + 6;        %the bottom side is taking in heat
    A(j+15,j+10) = A(j+15,j+10) - 8;
    A(j+15,j+5) = A(j+15,j+5) + 2;
    b(j+15) = 100;
end

for i=1:2
    for j = 2:4
        A(j+5*i,j+5*i) = -64;        %filling the interior points
        A(j+5*i,j+5*i+1) = 16;
        A(j+5*i,j+5*i-1) = 16;
        A(j+5*i,j+5*(i-1)) = 16;
        A(j+5*i,j+5*(i+1)) = 16;
    end
end

T = A\b;        %invert the matrix

diffT = T - 25;     %find the difference to the ambient temperature
T = T - 0.25^2 * 1 * diffT;

gridT = [T(1:5)'; T(6:10)'; T(11:15)'; T(16:20)']      %temperatures in a position grid
surf(gridT);

gridT =

    1.5625   30.5259   56.4405   77.4009   95.3125
    1.5625   33.0030   59.4893   79.8780   95.3125
    1.5625   40.4345   68.6357   87.3095   95.3125
    1.5625   58.5366   87.3095  105.4116   95.3125

2d Matrix method discussion

It is worth noting that the corners of the grid were set to be at the same temperature as the shorter sides which were held at constant temperatures. The above method is not strictly scalable but replacing some already substituted values for the size of the grid for a variable would allow for the matrix method to be applied on larger scales. Unfortunately, even for a 4x5 grid we require a 20x20 matrix. In comparison a 160x200 grid would require a 2000x2000 matrix which would seemingly be more difficult to invert. On the other hand it's a sparse matrix so let's see how matlab does.

clear;
tic
[A, b, gridT] = makingoutputshorter(40,50);     %Hiding the code and all output other than the surface graph.
toc

Elapsed time is 0.569434 seconds.

The above routine evaluates surprisingly quickly, however increasing the size of the matrix to (120,150) resulted in out of memory errors which is understandable. Unless Matlab was to compress arrays while they were in memory a 2000x2000 matrix would occupy 4MB. On the other hand, the required matrix for a (120,150) grid would occupy a GB of memory. While matlab may be able to efficiently invert this sparse matrix, on this computer Matlab isn't given enough memory for the matrix to be recorded.

Ultimately, the matrix method is a powerful, simple to implement technique for numerical approximation, however over increasing grid sizes the technique becomes quickly burdensome.