1D heat equation finite difference code improvement

I'm trying to solve a problem of a 1D heat equation using finite difference. I just would like to know if this code is well-written or can be improved?

Any help will be appreciated.

clear; clc; close all;


L = 1;
N = 5;                     % number of interior nodes
T0 = 100;                  % boundary at x=0
Tend = 200;                % boundary at x=L

h = L/(N+1);               % grid spacing
x_interior = (1:N)*h;      % interior x positions


e = ones(N,1);
A = (1/h^2) * spdiags([e -2*e e], -1:1, N, N); 


b = zeros(N,1);

b(1) = -T0 / h^2;

b(end) = -Tend / h^2;


T_interior = A \ b;


x_full = [0; x_interior'; L]';
T_full = [T0; T_interior; Tend];


fprintf('Grid spacing h = %.6f\n', h);

Grid spacing h = 0.166667

fprintf('Interior positions and temperatures:\n');

Interior positions and temperatures:

for i = 1:N
    fprintf(' x = %.6f,  T_%d = %.6f\n', x_interior(i), i, T_interior(i));
end

 x = 0.166667,  T_1 = 116.666667
 x = 0.333333,  T_2 = 133.333333
 x = 0.500000,  T_3 = 150.000000
 x = 0.666667,  T_4 = 166.666667
 x = 0.833333,  T_5 = 183.333333

fprintf('\nFull solution (including boundaries):\n');

Full solution (including boundaries):

for i = 1:length(x_full)
    fprintf(' x = %.6f,  T = %.6f\n', x_full(i), T_full(i));
end

 x = 0.000000,  T = 100.000000
 x = 0.166667,  T = 116.666667
 x = 0.333333,  T = 133.333333
 x = 0.500000,  T = 150.000000
 x = 0.666667,  T = 166.666667
 x = 0.833333,  T = 183.333333
 x = 1.000000,  T = 200.000000

figure;
plot(x_full, T_full, '-o', 'LineWidth', 1.4);
xlabel('x (m)');
ylabel('Temperature (^\circC)');
title('Steady-state 1D heat: finite-difference solution (N=5 interior)');
grid on;

T_analytic = 100 + 100*x_full; % since T(x) = 100 + 100 x
fprintf('\nAnalytic (linear) solution at grid points:\n');

Analytic (linear) solution at grid points:

for i=1:length(x_full)
    fprintf(' x = %.6f,  T_analytic = %.6f\n', x_full(i), T_analytic(i));
end

 x = 0.000000,  T_analytic = 100.000000
 x = 0.166667,  T_analytic = 116.666667
 x = 0.333333,  T_analytic = 133.333333
 x = 0.500000,  T_analytic = 150.000000
 x = 0.666667,  T_analytic = 166.666667
 x = 0.833333,  T_analytic = 183.333333
 x = 1.000000,  T_analytic = 200.000000

[Torsten]

I usually add the boundary conditions

T0 = 100
Tend = 200

as algebraic equations to the matrix A so that A becomes 7x7 instead of 5x5.

This way, I solve for all temperatures simultaneously and don't need to define

b(1) = -T0 / h^2;
b(end) = -Tend / h^2;

which - at first view - doesn't look obvious.

Instead, b becomes

b = [T0;0;0;0;0;0;Tend] 

which is easy to interpret.

But you can do as you like - your code looks fine as is.

Maybe a few more comments in the computational part can help to understand what you are doing if you look at the code in 3 months.

[Cesar]

Thank you, I'm trying is construct A and b and solve for the interior temperatures, and plot the temperature distribution along the rod including the boundary points. Also the linear system in matrix form AT = b for the interior unknownsT = [T1,T2,T3,T4,T5]T including the boundary conditions in the right-hand side vector b.

[Torsten]

I know this, and it is correct.

I only suggested to include the temperatures in the boundary points in the matrix A to make the equations more transparent:

T0 = 100;
Tend = 200;
A = [1 0 0 0 0 0 0;
     1 -2 1 0 0 0 0;
     0 1 -2 1 0 0 0;
     0 0 1 -2 1 0 0;
     0 0 0 1 -2 1 0;
     0 0 0 0 1 -2 1;
     0 0 0 0 0 0 1];
b = [T0;0;0;0;0;0;Tend];
T = A\b

T = 7×1

100.0000 116.6667 133.3333 150.0000 166.6667 183.3333 200.0000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

Answer 1

[Ayush]

Hi Cesar,

Your code seems correct and efficiently solves for the interior temperatures by incorporating the boundary conditions into the right-hand side vector.

Alternatively, as discussed in comments, you could include the boundary temperatures as unknowns and set up a larger system where the boundary conditions appear directly as equations in the matrix. This makes the system more transparent and the boundary values explicit in the solution vector, but is less common for numerical efficiency.

Both methods are valid and should yield the same temperature distribution along the rod.

Hope this helps!