1-D Transient Heat Conduction With No Heat Generation [FDM] [Crank-Nicholson]

I am currently coding a 1 D Transient Heat Conduction using Crank Nicholson method and I would like an expert opinion as to the accuracy of the result

%1-D Transient Heat Conduction With No Heat Generation [FDM][CN]
clear all
clc
clf
%% Variable Declaration
n = 21;                         %number of nodes
L = 1;                          %length of domain

A = sparse(n,n);                %initilizing Space
B = sparse(n,n);
Bx = zeros(n,1);
Ta = zeros(n,1);
Tb = zeros(n,1);

dx = L/(n-1);                    %domain element
dt = 20;                         %time step
tmax = 4000;                     %total Time steps (s)
t = 0:dt:tmax;
Tl = 100;                        %temperature at left face °C
Tr = 100;                        %temperature at right face °C

x =  linspace(0,L,n);            %linearly spaced vectors x direction
alpha = 1e-4;                    %thermal diffusivity (m^2/s)
r = alpha * dt /(2 *dx^2);       %for stability, must be 0.5 or less

%% Set Up Matrix
for i = 2 : n-1
    A(i,i-1) = -r;
        A(i,i ) = (1+2*r);       %Implicit Matrix
            A(i,i+1) = -r;
end
            A(1 ,1 ) = 1;
            A(n ,n ) = 1;

for j = 2 : n-1
    B(j,j-1) = r;
        B(j,j ) = (1-2*r);       %Explicit Matrix
            B(j,j+1) = r;
end
            B(1 ,1 ) = 1;
            B(n ,n ) = 1;
            
%% Boundry Condition
Bx(1,1) = Tl;            %Left Wall (Dirichlet conditions)
Bx(n,1) = Tr;            %Right Wall(Dirichlet conditions)

%% Solution
for k = 2 : length(t)   %time steps
    
        Bx(2 : n-1) = Tb(2 : n-1);
            Tb = B * Bx;
                        
        fprintf('Time t=%d\n',k-1);
    
end
            Ta = A\Tb;   %Solve CN Matrix
   
 %% Plot
 pos1 = [0.1 0.17 0.4 0.7];
    subplot('Position',pos1)
        plot(x,Ta);
        title('Nodes  Vs Temperature');
        xlabel('Nodes (i)');
        ylabel('Temperature °C');    
        grid on
        grid minor
    
    a=0:L;
    d=0:0;
    
 pos2 = [0.51 0.45 0.45 0.1];
    subplot('Position',pos2) 
    ax = gca;
    imagesc(a,d,Ta', 'Parent', ax);
    ax.YAxis.Visible = 'off';
        
        title('Temperature Gradient');
        xlabel('Nodes (i)');
        h = colorbar;
        ylabel(h, 'Temperature °C')
        colormap jet

[darova]

please show the main equation

[Leroy Coelho]

[darova]

Can you show original equation (in LaTeX form)?

[Leroy Coelho]

-r\cdot{T}_{i-1,j+1} + (1+2\cdot{r})\cdot{T}_{i,j+1} -r\cdot{T}_{i+1,j+1}=r\cdot{T}_{i-1,j} + (1-2\cdot{r})\cdot{T}_{i,j} +r\cdot{T}_{i+1,j}

r= \frac {alpha\cdot{dt}}{2\cdot{dx^2}}

Answer 1

[Leroy Coelho]

Code Update and cross checked with other sources for accuracy

%% Solution
for k = 2 : length(t)   %time steps
    
     
            Tb = B * Bx;
            
                    Ta = A\Tb;   %Solve CN Matrix    
            
                            Bx(2 : n-1) = Ta(2 : n-1); 
                                    
                                    fprintf('Time t=%d\n',k-1);
    
end
   

Answer 2

[Bjorn Gustavsson]

My suggestion is that you add a line:

Tall(:,k) = Tb;

at the end of your "time-step"-loop. Then after you're done with it all you can simply calculate the time-derivative and second-order spatial derivative of your solution. Then you can check how well those satisfy the heat-equation by simply calculating the required gradients:

subplot(2,2,1)                    
imagesc(t,x,Tall)                          
colorbar
subplot(2,2,2)
imagesc(t,x,dT1)                           
caxis([-1 1]/100)
caxis([0 .2])    
colorbar
[d2T21,d2T22] = gradient(dT2,t,x);
subplot(2,2,3)
imagesc(t,x,d2T22)                         
colorbar           
subplot(2,2,4)
imagesc(t,x,(dT1-alpha*d2T22)./Tall)
caxis([-1 1]/100)  
colorbar

Here I've not bothered with calculating the gradients properly centred for your comparison - but that surely is your job...

HTH

[Leroy Coelho]

thank you for your response ,much appreciated