Consider the heat equation problem in the theory lecture:

We can plot the initial condition

x = [0:0.1:4];

y = 2 - abs(x-2);

plot(x,y,'linewidth', 4)

Now lets use the Crank-Nicolson scheme (the gold standard for finite differences for the heat equation) to solve this. We saw in the theory lecture that our finite difference scheme simplifies to

this indicates we have a tridiagonal system on both sides of the equation, but since we know we simply do the matrix multiplication to give us a vector on the right hand side. For more detail take a look at the theory lecture.

dt = 0.01;

dx = 0.01;

x = [0:dx:4];

X = length(x);

k = 2;

u = 2 - abs(x-2);

u = u';

mu = k*dt/(2*dx^2);

u_a = 0; u_b = 0; %Dirichlet boundary conditions

main_diag = (1+2*mu)*ones(1,X-2);

second_diag = -mu*ones(1, X-3);

A = diag(main_diag) + diag(second_diag,1) + diag(second_diag, -1);

main_diag = (1-2*mu)*ones(1,X-2);

B = diag(main_diag) - diag(second_diag,1) - diag(second_diag, -1);

plot(x,u)

axis([0 4 0 2])

pause(0.1)

hold on

for t = dt:dt:1

b = B*u(2:end-1);

b(1) = b(1) + mu*u_a;

b(end) = b(end) + mu*u_b;

u(2:end-1) = A\b;

plot(x,u)

axis([0 4 0 2])

pause(0.1)

end

hold off

Lets try it with an initial condition that doesn't have a kink.

dt = 0.01;

dx = 0.01;

x = [0:dx:4];

X = length(x);

k = 2;

u = 2 - (x-2).^2;

u = u';

mu = k*dt/(2*dx^2);

u_a = 0; u_b = 0; %Dirichlet boundary conditions

main_diag = (1+2*mu)*ones(1,X-2);

second_diag = -mu*ones(1, X-3);

A = diag(main_diag) + diag(second_diag,1) + diag(second_diag, -1);

main_diag = (1-2*mu)*ones(1,X-2);

B = diag(main_diag) - diag(second_diag,1) - diag(second_diag, -1);

plot(x,u)

axis([0 4 0 2])

pause(0.1)

hold on

for t = dt:dt:1

b = B*u(2:end-1);

b(1) = b(1) + mu*u_a;

b(end) = b(end) + mu*u_b;

u(2:end-1) = A\b;

plot(x,u)

axis([0 4 0 2])

pause(0.1)

end

hold off