Crank Nicolson for non homogeneous heat equation

8 views (last 30 days)
Hey folks, I have ran into a problem. I am coding for the c-n method and using this base code with the given starting values for I get proper results
%u_t-u_xx=0
h=0.1;
k=0.01;
alpha=1; lambda=(alpha^2)*k/(h^2)
a=0; b=1; %0<x<1
c=0; d=1; %0<t<1
x=a:h:b; t=c:k:d;
n=length(x); m=length(t);
w=zeros(n,m);
w(1,:)=zeros(1,m); w(end,:)=zeros(1,m); %boundary conditions
w(:,1)=sin(pi*x); %initial conditions
f=zeros(n,m);
% u=exact solution
for i=1:n
for j=1:m
u(i,j) = exp(-(pi^2)*t(j))*sin(pi*x(i));
end
end
%C-N
A=diag((1+lambda)*ones(n-2,1))+diag((-lambda/2)*ones(n-3,1),-1)+diag((-lambda/2)*ones(n-3,1),1);
B=diag((1-lambda)*ones(n-2,1))+diag((lambda/2)*ones(n-3,1),-1)+diag((lambda/2)*ones(n-3,1),1);
for j=2:m
b=B*w(2:end-1,j-1);
w(2:end-1,j)=A\b;
end
figure(1)
mesh(u)
figure(2)
mesh(w)
My problems start when I move to a nonhomogenous heat equation. such as
%u_t-u_xx=f
h=0.2;
k=0.01;
alpha=1; lambda=(alpha^2)*k/(h^2);
a=0; b=1; c=0; d=1;
x=a:h:b; t=c:k:d;
n=length(x); m=length(t);
w=zeros(n,m);
w(1,:)=zeros(1,m); w(end,:)=sin(1)*cos(t); %boundary conditions
w(:,1)=sin(x); %initial conditions
f=zeros(n,m);
for i=1:n %Right hand side of function
for j=1:m
f(i,j)=-sin(x(i))*sin(t(j))+sin(x(i))*cos(t(j));
end
end
%solution
u=zeros(n,m);
for i=1:n
for j=1:m
u(i,j) = cos(t(j))*sin(x(i));
end
end
%C-N
A=diag((1+lambda)*ones(n-2,1))+diag((-lambda/2)*ones(n-3,1),-1)+diag((-lambda/2)*ones(n-3,1),1);
B=diag((1-lambda)*ones(n-2,1))+diag((lambda/2)*ones(n-3,1),-1)+diag((lambda/2)*ones(n-3,1),1);
b=zeros(n-2,1);g1=w(1,:);g2=w(end,:);
for j=2:m
b(1)=w(2,j-1)+k*f(2,j)+lambda*g1(j);
for i=2:length(b)-1
b(i)=w(i+1,j-1)+k*f(i+1,j);
end
b(end)=w(length(b)+1,j-1)+k*f(length(b)+1,j)+1.0455*lambda*g2(j); %1.0455 scalor was a manually and arbitrarily chosen value as I try to mend this abomination of a code
z=B*b;
w(2:end-1,j)=A\(z);
end
figure(1)
mesh(u)
figure(2)
mesh(w)
I have tried multiple attempts that have not gotten quite as close. This will also break down at varying time steps which should absolutely not be the case. Does anyone have any tips?

Answers (0)

Products


Release

R2019b

Community Treasure Hunt

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

Start Hunting!