Heat equation by Crank Nicolson

80 Views Asked by At

I was doing the comparison between the exact solution $u (x, t)$ and the approximate solution $w (x, t) $ of the heat equation, at the point $ x = 0.5$ and making $t$ variate in the interval $ [0,1]$ only that the values of $w$ only the first value matches and the others come out differently. If anyone could help me, I appreciate in advance.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function w=Crank_Nicolson(xl,xr,yb,yt,M,N) %% Equation by Crank-Nicolson

% input : [xl,xr], [yb,yt],  M,  N
f=@(x) exp(-0.5*x);
l=@(t)exp(t);
r=@(t)exp(t-0.5);
D=4;
h=(xr-xl)/M; k=(yt-yb)/N; m=M-1; n=N;

sigma=(D*k)/(h*h);
a=diag((2+2*sigma)*ones(m,1))+diag(-sigma*ones(m-1,1),1);
a=a+diag(-sigma*ones(m-1,1),-1);% matriz a tridiagonal

b=diag((2-2*sigma)*ones(m,1))+diag(sigma*ones(m-1,1),1); % matrix b trid
b=b+diag(sigma*ones(m-1,1),-1);

lside=l(yb+(0:n)*k);rside=l(yb+(0:n)*k);
w(:,1)=f(xl+(1:m)*h)';
u=@(x,t) exp(t-(0.5).*x); % solution exact
t=yb;

fprintf('\n t_j \t \t  u(0.5,j) \t w(0.5,j)  \n');

for j=1:n 
    sides=[lside(j)+lside(j+1);zeros(m-2,1);rside(j)+rside(j+1)];
    w(:,j+1)=a\(b*w(:,j)+sigma*sides); 
end

w=[lside;w;rside]; 

for j=1:n+1 
  fprintf('\n %5.4f \t %5.4f \t %5.4f \n',t,u(0.5,t),w(6,j));
  t=yb+j*k; 
end

>>w=Crank_Nicolson(0,1,0,1,10,10);