how to solve 1D heat equation by Crank-Nicolson method

I need to solve a 1D heat equation by Crank-Nicolson method . The tempeture on both ends of the interval is given as the fixed value u(0,t)=2, u(L,t)=0.5. I solve the equation through the below code, but the result is wrong. Attached figures are the correct result. I don't know why? Could you please anyone offer me a hand? Thanks a lot.

clear all;
L=1;
dx=0.1;
imax=L/dx+1;
X=linspace(0,L,imax);
% inital value
f0 = 2-1.5.*X+sin(pi.*X);
dt=0.05;
nstep=15;
D=1.44;
alpha=D*dt/(dx)^2;


e  = ones(imax,1);
B  = [-alpha*e 2*(1+alpha)*e -alpha*e];
Lx = spdiags(B,[-1 0 1],imax,imax);
B  = [alpha*e 2*(1-alpha)*e alpha*e];
Rx = spdiags(B,[-1 0 1],imax,imax);

%%CN method:
u=zeros(imax,nstep+1);
u(:,1)=(f0)';

for n=2:nstep+1
    u(:,n)=Lx\(Rx*u(:,n-1));
    % boundary value
    u(1,n)=2;
    u(end,n)=0.5;
end
% display contour plot
t=[0:nstep];
[Xg,tg] = meshgrid(X,t);
figure;
contourf(Xg,tg,u.');

[reyhdh saleh]

I have problem :

Undefined function or variable 'Crank'

[reyhdh saleh]

what is mean

[reyhdh saleh]

%This program is meant to test out the Crank-Nicolson scheme using a simple

%nonlinear diffusion scheme.

n=5000;m=30;

t=linspace(0,20,n);% set time and distance scales

x=linspace(0,1,m);

dx=x(2)-x(1);dt=t(2)-t(1); %Increments

s=dt/(2*dx^2);%Useful for the solution.

u=zeros(n,m); %set up solution matrix

p=s*ones(1,m-1); q=-(1+2*s)*ones(1,m);

Z=diag(p,1)+diag(p,-1)+diag(q);

A=Z(2:m-1,2:m-1);

%Add in intial condition:

u(1,:)=exp(-5*(x-0.5).^2);

v=zeros(m-2,1);

%Solve the system

for i=2:n-1

%Construct the RHS for the solution

for j=2:m-1

v(j-1)=s*u(i-1,j+1)-(2*s+1)*u(i-1,j)-s*u(i-1,j-1);

end

%Solve for the new time step

w=A\v;

u(i,2:m-1)=w;

u(i,1)=u(i,2);

u(i,end)=u(i,end-1);

end

I have problem ....

Undefined function or variable 'Crank'

what does it mean

[Torsten]

We don't know since the string "Crank" (most probably for Crank-Nicholson) does not appear in the code you posted.

You can run the code from above just as it is as a script in MATLAB.

Answer 1

[Jiali]

I figure out that the boundary is dealt with incorrectly. After re-deriving the matrix to include the boundary value, I finally get the correct result.

[darova]

Please post it

Answer 2

[Jiali]

Dear Darova,

   The below code is What I figure out. Hopefully it helps someone.
   
   

u=zeros(length(X),nstep+1);
u(:,1)=(f0)';
% boundary value
u(1,:)=2;
u(end,:)=0.5;

imax=length(X)-2;   % remove boundary point
e  = ones(imax,1);
A  = [-alpha*e 2*(1+alpha)*e -alpha*e];
Lx = spdiags(A,[-1 0 1],imax,imax);
B  = [alpha*e 2*(1-alpha)*e alpha*e];
Rx = spdiags(B,[-1 0 1],imax,imax);

%% CN method:
for n=2:nstep+1
    % adjust the right matrix based on the fixed value on boundary
    temp=zeros(imax,1);
    temp(1)=alpha*2*u(1,n);
    temp(end)=alpha*2*u(end,n);    
    u(2:end-1,n)=Lx\(Rx*u(2:end-1,n-1)+temp);

end