Crank Nicolson for non homogeneous heat equation

3 次查看（过去 30 天）

Ethan Boyce 2020-12-7

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/684908-crank-nicolson-for-non-homogeneous-heat-equation

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?