How can I modify my code from the Explicit Method to the Implicit Method? Thank you

How would I go about changing this from Explicit Method to the Implicit Method?

Thank you

Answer 1

[Shantanu Dixit]

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