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

154 次查看（过去 30 天）
Jiali
Jiali 2020-2-18
评论： Torsten 2022-1-1
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 个评论
显示 2更早的评论
reyhdh saleh
reyhdh saleh 2022-1-1
%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
Torsten 2022-1-1
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.

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Jiali
Jiali 2020-2-20
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.

更多回答（1 个）

Jiali
Jiali 2020-3-16
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

类别

在 Help Center 和 File Exchange 中查找有关 Programming 的更多信息

标签

产品


版本

R2015a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by