I succesfully created the following script for 1D conduction using the explicit discretization thanks to the natural way to write in on paper using matrix operations. I'm trying to write down (using pen and paper!) the matrix operations for the implicit discretization (the problem is the same for Crank Nicolson)... Obviously a change on the time step doesn't work... e.g. fourier * (T(i-1, k+1) -...)
%% Finite differences - Explicit discretization
close all; clear; clc;
%% Problem's input
alpha = input('Thermal diffusivity [mc/s] = ');
L = input('Lenght of the bar [m] = ');
Tleft = input('Imposed temperature on the left [K] = ');
Tright = input('Imposed temperature on the right [K] = ');
T0 = input('Temperature before the perturbation [K] = ');
m = input('How many nodes for along x? ');
n = input('How many time steps? ');
tmax = input('How long do you want to see? ');
%deltat = tmax / n;
deltax = L / m;
%fourier = alpha*deltat/deltax^2 % Fron page 222 Mills: Fourier has to be minor than 0.5 to avoid divergent oscillation
% That means 0.5*deltax^2/alfa = deltat to avoid divergent oscillation
% Fix fourier 0.4
fourier = 0.4;
deltat = ( fourier * deltax^2 ) / alpha;
% At the beginning, we have a matrix of the temperature with the quite
% temperature everywhere along x!
T = ones(m, n+1)*T0;
% Live plot
figure;
h = mesh(linspace(0, L, m), linspace(0, tmax, n+1), T');
xlabel('Position [m]');
ylabel('Time [s]');
zlabel('Temperature [K]');
title('Temperature in the bar',LineWidth=19);
%% Iterative calculations
for k = 1:n
% For step k, it is calculated the temperature at node i for time step
% k+1
for i = 2:m-1
T(i, k+1) = T(i, k) + fourier * (T(i-1, k) - 2*T(i, k) + T(i+1, k));
end
% Border conditions
T(1, k+1) = Tleft;
T(end, k+1) = Tright;
% Updated mesh plot
set(h, 'ZData', T');
title(['Temperature distribution at time t = ', num2str(k * deltat)]);
drawnow;
end
%% Final plot
figure(1)
plot(linspace(0, L, m), T(:, end));
xlabel('Position [m]');
ylabel('Temperature [K]');
title(['Temperature at time t = ', num2str(tmax,2)]);
pause;
