Crank-Nicolson method for solving a simple diffusion/heat problem with time-dependence.

Crank-Nicolson.m

%% Startup

clc
clear all
close all

%% Definitions

% Axis
dx = 0.1;
dt = 0.01;
num_x_points = 21;
num_t_points = 11;

x_axis = linspace(0, 2, num_x_points);
t_axis = linspace(0, 0.1, num_t_points);
z_axis = linspace(-2, 2, 20);

% Conditions
% A sample initial condition is given, apply as needed
initial = sin(x_axis*pi) + sin(3*pi*x_axis);

% Left and right temperature boundaries
T_l = 0;
T_r = 0;

%% Crank-Nicolson Method

% Set up answer matrix
T = zeros(num_t_points, num_x_points);
% N.B. The columns are x coordinates, the rows are t coordinates

% Apply boundary conditions to matrix
T(1,:) = initial;
T(:,1) = T_l;
T(:,end) = T_r;

% Set up solver
a = 1;
r = (a*dt)/(2*dx^2);

% Diagonals for M
central = ones(num_x_points - 2, 1) * (1 + 2*r);
upper = ones(num_x_points - 3, 1) * (-r);
lower = ones(num_x_points - 3, 1) * (-r);

M = diag(central) + diag(upper, 1) + diag(lower, -1);

% b
% Using the sneaky, offset matrix way
% Iterate, keep solving for the temperature in the row
for i = 2:num_t_points
    t_t = [0; T(i-1, 2:end-1)'; 0];

    t_reg = t_t(2:end-1);
    t_minus = t_t(1:end-2, :);
    t_plus = t_t(3:end, :);

    b = r*t_plus + (1-2*r)*t_reg + r*t_minus;
    
    b(1) = b(1) + r*T_l;
    b(end) = b(end) + r*T_r;

    % Solve the equations by x = M \ b as usual
    t = [T_l; M \ b; T_r]';
    T(i, :) = t;
end

% Plot the solution
surf(x_axis, t_axis, T, 'EdgeColor', 'None')

xlabel('x')
ylabel('Time (s)')
zlabel('Temperature')