The only thing that changes is the way we build a new row from the previous row. To do this using Crank Nicolson we call the tridiagonal system solver tridiag.m. The program below is the same as for Euler, but with the addition of the definition of 3 vectors a,b,c used in tridiag, and a change in the way we go from one row to the next. The program takes h,k as input, returning (1) the approximate value of u(0,1) and (2) the time of unimodality:
function [v1 v2]=t5q5(h,k)
% set up the grid
T = 2 ;
Xmin = -2 ;
Xmax = 2 ;
N = (Xmax-Xmin)/h ; % assume this is an integer
M = T/k ; % assume this is an integer
x = Xmin + (0:N)*h ;
% t = (0:M)*k ;
% vectors making the tridiagonal matrix:
a=(1+k/h^2)*ones(1,(N-1));
b=-k/2/h^2*ones(1,(N-2));
c=b;
% prepare matrix to hold the solution
u = zeros(M+1,N+1);
% no need to impose boundary conditions as they are zero
u(1,:) = x.^4 - 2*x.^2 - 8;
% Crank-Nicolson method to evolve
for j=1:M
u(j+1,2:N) = tridiag(a,b,c,u(j,2:N) + k/2/h^2*( u(j,3:(N+1)) - 2*u(j,2:N) + u(j,1:(N-1)) ) ) ;
if ( u(j+1,N/2+1)-u(j+1,N/2) ) * ( u(j,N/2+1)-u(j,N/2) ) < 0
v2=(j+1)*k;
end
end
% plot the solution at times 0, 0.5, 1, 2
% color blue,red,green and black
figure(1)
hold on
plot(x,u(1,:),'b');
plot(x,u(M/4+1,:),'r');
plot(x,u(M/2+1,:),'g');
plot(x,u(M+1,:),'k');
v1=u(M/2+1,N/2+1);
Some results: For h=1/10, 1/20, 1/40 it was possible to use k=1/10, k=1/8 and k=1/6 and get sensible looking plots. Here are some approximations of u(0,1):
h=1/10 k=1/10 -5.66284 error 203 h=1/20 k=1/20 -5.66132 51 h=1/40 k=1/40 -5.66094 13 h=1/80 k=1/80 -5.66084 3 h=1/1000 k=1/1000 -5.66081
Taking the last result as accurate to 6 significant figures (which it is), we can find the errors in the other results (listed) and we see that when we halve both h and k we quarter the error, consistent with the error being O(h^2) + O(k^2).
Crank Nicolson gives much better results for much less effort. We see that we can get 4 figure accuracy with h=1/20 and k=1/20, and we could not get this with the Euler method with k=1/4000. However - if we wish to discuss the first time at which the solution is unimodal and we work with k=1/20, then using the method we used before we will never get the time more accurately than 1/20! So sometimes large time steps is not an advantage.