Crank Nicolson for non homogeneous heat equation
3 visualizzazioni (ultimi 30 giorni)
Mostra commenti meno recenti
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?
0 Commenti
Risposte (0)
Vedere anche
Categorie
Scopri di più su Geometry and Mesh in Help Center e File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!