Problems with 1D heat diffusion with the Crank Nicholson method

Hi everyone, I am trying to implement the code for 1D heat diffusion with the Crank Nicholson method. but it turned out to an error at the temperature matrix. could you please help me with the code? I would really appreciate it.

% Heat Diffusion in 1D within the Crank-Nicholson Method
%du/dt=cond*d2u/d2t
%dU/dt is the net change in the total energy of the system
%the flux of heat is proportional to the gradient of the temperature

clear all;
clc;

% Parameters to define space and time

L = 0.05; % Lenght/Thickness [m]

Tmax =900.; % Final time [s]

% Parameters needed to solve the equation

Nt = 1000; % Number of time steps

t = Tmax/Nt; %time differential

N = 25.; % Number of space steps

h = L/N; %space differential

cond = 0.4; % Conductivity

rho= 900; % density

ce= 800; %specific heat

b = cond*t/(rho*ce*h*h); % Parameter of the method 

%solution "u" is temperature , 
%that is why I use specific heat (ce) and density (rho)

% Initial temperature 
Tamb= 25; % constant ambient temperature

for i0 = 1:N+1

      dx(i0) =(i0-1)*h; %position vector (used in graphics)

      u(i0,1) =Tamb; % initial ambient temperature 

end

% BC: Temperature at the boundary
T1=25; %inital temperature for boundary condition
T2=35; %final temperature for boundary condition

T0= linspace(T1,T2,Nt+1); 
%linearly gradient between 2 known temperatures in x=0

TL= linspace(T1,T2,Nt+1);
%linearly gradient between 2 known temperatures in x=L

for j0=1:Nt+1

      u(1,j0) = T0(j0); %BC in the left extreme x=0

      u(N+1,j0) = TL(j0); %BC in the right extreme x=L

      dt(j0) = (j0-1)*t; %time vector (used in graphics)

end

% Defining the Matrices M_right and M_left in the CN method

aal(1:N-2)=-b;

bbl(1:N-1)=2.+2.*b;

ccl(1:N-2)=-b;

MMl=diag(bbl,0)+diag(aal,-1)+diag(ccl,1);

aar(1:N-2)=b;

bbr(1:N-1)=2.-2.*b;

ccr(1:N-2)=b;

MMr=diag(bbr,0)+diag(aar,-1)+diag(ccr,1);

% Implementation of the Crank-Nicholson method

for j0=2:Nt+1             % Time Loop

      uu=u(2:N,j0-1);

      u(2:N,j0)=inv(MMl)*MMr*uu;

end

The problem is that the solution vector "u" seems to end with an error (illogical solution), or divided by 0.

Thank you in advanced for the help.

Answer 1

[Torsten]

rho and ce are undefined. Thus b is undefined.

Best wishes

Torsten.

[Royito]

Thank you Torsen for the comment. However, that is not the error since I defined them already in my code but I did not write it down here.

Sorry for the misunderstood. I edit the question so it is not confusing for anyone else.

[Torsten]

In the for-loop at the end of your code, you don't use u(1,:) and u(N+1,:) (thus the boundary conditions for your problem). Consequently, the resulting solution must be wrong.

Best wishes

Torsten.

[Royito]

I am not wrong, the calculation with the crank nicholson method is applied only to the internal nodes u(2:N). That is why i did not used u(1,:) and u(N+1,:). Is that correct?

[Torsten]

No, it can't be correct if you don't use the boundary conditions for the calculation of u in the interior of the interval.

Best wishes

Torsten.

[Royito]

I managed to change the code in order to insert also the boundary conditions. THANK YOU!

[kolade aboluwade]

please can you drop the code for crank nicolson method