how to solve 1D heat equation by Crank-Nicolson method
362 views (last 30 days)
Show older comments
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.');
4 Comments
reyhdh saleh
on 1 Jan 2022
%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
on 1 Jan 2022
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.
Accepted Answer
More Answers (1)
See Also
Categories
Find more on Ordinary Differential Equations in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!