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

조회 수: 284 (최근 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
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 CenterFile Exchange에서 Ordinary Differential Equations에 대해 자세히 알아보기

제품


릴리스

R2015a

Community Treasure Hunt

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

Start Hunting!

Translated by