Solution diverges for 1D heat equation using Crank-Nicholson

6 ビュー (過去 30 日間)
Matthew Hunt
Matthew Hunt 2018 年 7 月 3 日
コメント済み: Torsten 2018 年 7 月 4 日
I am trying to solve the 1D heat equation using the Crank-Nicholson method. I have managed to code up the method but my solution blows up. I'm using Neumann conditions at the ends and it was advised that I take a reduced matrix and use that to find the interior points and then afterwards. Would this work? My code is blowing up and it should work but for some reason it isn't and I can't see the reason why. Any suggestions? Code below:
%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

回答 (1 件)

Torsten
Torsten 2018 年 7 月 3 日
%This program is meant to test out the Crank-Nicolson scheme using a simple
%nonlinear diffusion scheme.
n=5000;m=150;
t=linspace(0,10,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
A=zeros(m,m);
A(1,1) = 1;
A(1,2) = -1;
for j=2:m-1
A(j,j-1) = -s;
A(j,j) = (1+2*s);
A(j,j+1) = -s;
end
A(m,m-1)= -1;
A(m,m) = 1;
%Add in intial condition:
u(1,:)=exp(-5*(x-0.5).^2);
v=zeros(m,1);
%Solve the system
for i=2:n-1
%Construct the RHS for the solution
for j=2:m-1
v(j)=s*u(i-1,j+1)+(1-2*s)*u(i-1,j)+s*u(i-1,j-1);
end
%Solve for the new time step
w=A\v;
u(i,:)=w;
end
end
  4 件のコメント
Matthew Hunt
Matthew Hunt 2018 年 7 月 4 日
Actually, I saw the error in my code which still exists in the code you gave me. The problem wasn't in the allocation of the matrix A which you seemed to think but was in the definition of v which you didn't change.
Torsten
Torsten 2018 年 7 月 4 日
The definition of v has changed - take a closer look.

サインインしてコメントする。

カテゴリ

Help Center および File ExchangeProgramming についてさらに検索

製品


リリース

R2018a

Community Treasure Hunt

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

Start Hunting!

Translated by