MATLAB: How to solve 1D heat equation by Crank-Nicolson method

MATLABpartial differential equation

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.');

Best Answer

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.

Related Solutions

MATLAB: Hello, how to solve the problem for matrix that a is diagonal matrix and u is random column matrix. I would like to solve the equation to get u=a*u inside the loop. But , I got index in the position 2 is invalid answer. Could you please help me

Matlab's indexing starts at 1 not zero, so it is complaining at u(:,k) when k is zero.

MATLAB: Hi, I’m trying to solve the heat eq using the explicit and implicit methods and I’m having trouble setting up the initial and boundary conditions.The equation is : du/dt=d^2u/dx^2, initial condition u(x,0)=x, boundary conditions u(0,t)=1 du/dx(1,t)

The matlab code posted here is based on the formulation of the explicit method of the finite difference method. I tried to compare the solution to that obtained from using matlab's pdepe solver to ensure that the coding was done correctly. The both results seem to agree, but the solution based on matlab's pdepe solver came out superior in that the mesh points are much finer than the explicit method; which is true since the pdepe function is highly accurate. You can experiment with the parameters to see if you can get a finer mesh points.

The output of the code.

Output of matlab's pdepe solver.

    %%Numerical solution of the Heat equation using the explicit Finite
    % difference method
    clear variables
    close all
    %1 Function parameters
    N  = 500;
    Lx = 50;
    dx = Lx/(N-1);
    x  = 0:dx:Lx;
    alpha = .25;   % To insure stability alpha = dt/dx^2 < = 1/2
    % 2. Time vector
    M = 300;
    tf = 100;     
    dt = tf/(M-1);
    t = 0:dt:tf;
    % 4. Initial and boundary conditions
    f  = @(x) x;   % initial cond. f(x)
    g1 = @(t) 1;     % boundary conditions g1(t) and g2(t)
    g2 = @(t) 1;
    % 5. Inilization of the heat equation
    u      =  zeros(N,M);
    u(:,1) =  f(x);
    u(1,:) =  g1(t);
    % 6. Implementation of the explicit method
    for j= 1:M-1       % Time Loop 
         for i= 2:N-1  % Space Loop
           u(i,j+1) = alpha*(u(i-1,j))+(1-2*alpha)*u(i,j) + alpha*u(i+1,j);
         end
         % Vectorize the inner for loop with this line to gain some speed.
         % u(2:N-1,j+1) = alpha*(u(1:N-2,j))+(1-2*alpha)*u(2:N-1,j) + alpha*u(3:N,j);
     % Insert boundary conditions for i = 1 and i = N here.
     u(2,j+1) =  alpha*u(1,j) + (1-2*alpha)*u(2,j) + alpha*u(3,j);
     u(N,j+1) =  2*alpha*u(N-1,j)+(1-2*alpha)*u(N,j)+ 2*alpha*dx*g2(t);
    end
    % Plot results
    figure
    plot(x,u(1:end,1:30),'linewidth',2);
    a = ylabel('Temperature');
    set(a,'Fontsize',14);
    a = xlabel('x');
    set(a,'Fontsize',14);
    a=title(['Using The Explicit Method - alpha =' num2str(alpha)]);
    legend('Explicit soln')
    set(a,'Fontsize',16);
    xlim([0 1]);
    grid;
    figure
    [X, T] = meshgrid(x,t);
    s2 = mesh(X',T',u);
    title(['3-D plot of the 1D Heat Eq. using the Explicit Method - alpha = ' num2str(alpha)])
    set(s2,'FaceColor',[1 0 1],'edgecolor',[0 0 0],'LineStyle','--');
    a = title('Explicit finite difference solution of the 1D Diffusivity Equation');
    set(a,'fontsize',14);
    a = xlabel('x');
    set(a,'fontsize',20);
    a = ylabel('y');
    set(a,'fontsize',20);
    a = zlabel('z');
    set(a,'fontsize',20);
    xlim([0 2])
    zlim([0 2.5])

Best Answer

Related Solutions

MATLAB: Hello, how to solve the problem for matrix that a is diagonal matrix and u is random column matrix. I would like to solve the equation to get u=a*u inside the loop. But , I got index in the position 2 is invalid answer. Could you please help me

MATLAB: Hi, I’m trying to solve the heat eq using the explicit and implicit methods and I’m having trouble setting up the initial and boundary conditions.The equation is : du/dt=d^2u/dx^2, initial condition u(x,0)=x, boundary conditions u(0,t)=1 du/dx(1,t)

Related Question