Hi Xavier,
To change from the explicit method to the implicit method for calculating temperature distribution across a slab, here’s a high-level explanation of the steps needed to implement the implicit finite difference method:
- Define Parameters and Set Initial Conditions: Set up the slab length, time steps, and boundary conditions, same as in the explicit method.
- Construct the Coefficient Matrix (A): Build a tridiagonal matrix based on lambda values, which will serve as the system's coefficient matrix for the implicit approach.
- Set Up the Right-Hand Side Vector (b) at Each Time Step: Use the temperature values from the previous time step to set up the vector b, incorporating boundary conditions.
- Solve the System and Update Temperatures: Solve the matrix equation '
' to get the temperature distribution, updating the results for each time step.
%% Define Parameters and Set Initial Conditions
Length = 1;
delta_t = 0.001;
total_time = 0.1;
N = total_time/delta_t;
M = 20;
delta_x = Length/M;
Lamda = (delta_t)/(delta_x^2);
% Initialize temperature arrays
T = zeros(M+1, N+1);
% Initial conditions
T(1,:) = 350; % Left boundary
T(M+1,:) = 410; % Right boundary
% Initial temperature distribution
for i = 1:M+1
x(i) = (i-1)*delta_x;
T(i,1) = 200*(i-1)*delta_x*sin(pi*(i-1)*delta_x);
end
%% Create Coefficient Matrix for Implicit Solution
r = Lamda;
A = zeros(M-1, M-1);
for i = 1:M-1
A(i,i) = 1 + 2*r; % main diagonal
if i > 1
A(i,i-1) = -r; % lower diagonal
end
if i < M-1
A(i,i+1) = -r; % upper diagonal
end
end
%% Iterate with Matrix Solutions at Each Time Step
for n = 1:N
% Create the right-hand side vector b
b = zeros(M-1, 1);
for i = 1:M-1
b(i) = T(i+1, n);
if i == 1
b(i) = b(i) + r * T(1, n+1); % Left boundary condition
end
if i == M-1
b(i) = b(i) + r * T(M+1, n+1); % Right boundary condition
end
end
solution = A\b;
% Update temperature array
for i = 2:M
T(i, n+1) = solution(i-1);
end
end
figure(1)
plot(x, T(:,end), 'b-', 'LineWidth', 2)
title('Temperature Distribution (Implicit Method)')
xlabel('Length of Slab')
ylabel('Temperature in Kelvin @ Time = 0.1 @ \Delta t = 0.001')
grid on
- For more details, refer to the following MathWorks documentation on 'mldivide' https://www.mathworks.com/help/matlab/ref/double.mldivide.html
