# Where am I making a mistake in solving the heat equation using the spectral method (Chebyshev's differentiation matrix)?

I would like to numerically solve the following heat equation problem:

$$u_t = \Bigg(2{a \over l}\Bigg)^2 u_{xx} \tag 1$$

$$x \in [ -1, 1 ] \tag 2$$ $$u(x, 0) = 0 \tag 3$$ $$u(1, t) = A \sin \Bigg({a^2 \over l^2}t\Bigg) \tag 4$$ $$u_x(-1, t)=0 \tag 5$$ $$u(x, 0) = 0 \tag 6$$

where $$u(x,t)$$ is the heat function, $$l$$ is the rod's length, $$a$$ is the heat constant and $$A$$ is the amplitude of the boundary condition signal. The analitical solution of this problem is:

$$u(x,t) = A \sin \Bigg({a^2 \over l^2}t\Bigg) + \sum_{n=0}^\infty {4 \sin ^2 \big({w_n \over 2}\big) \over w_n + w_n^5} \sin \Bigg( w_n {1 - x \over 2} \Bigg) \Bigg ( w_n^2 e^{-w_n^2 {a^2 \over l^2}t} - w_n^2 \cos \bigg({a^2 \over l^2}t\bigg) - \sin \bigg({a^2 \over l^2}t\bigg) \Bigg ) \tag 7$$

where $$w_n = n \pi + \pi / 2$$. I graphically checked the analitical solution:

However, my numerical solution seems to fail. Using Matlab, here is what I did.

First, I simply defined my simulation time parameters using the code listed below.

%% Simulation time parameters

nt    = 2000001;                % Number of time points
t_max = 1;                      % Simulation time
t     = linspace(0, t_max, nt);
dt    = t(2) - t(1);            % Time step


Second, I defined the constants I am using for the simulation:

%% Constants

a = 1;            % Heat constant
A = 100;          % Amplitude of the boundary condition signal
l = a / sqrt(pi); % Rod length


Next, I formulated my numerical solution using the code listed below. Here, I define the number of spatial points and I get the Chebyshev's differentiation matrix, as well as the spatial vector defined at Chebyshev's points. After that, I square the differentiation matrix because I need the second spatial derivative and I also change the last row of it because I need to satisfy the Neumann boundary condition. Then, I define my linear differential operator $$L$$ and just scale it down for the finite difference time loop. In the end, I solve the temporal ODE using the finite difference method while enforcing the boundary condition inside the time marching loop.

%% Numerical solution

N = 30;

[D, x] = cheb(N);     % D is the Chebyshev's differentiation matrix
nx     = length(x);   % Number of spatial points

D2 = D^2;
D2(N + 1, :) = D(N + 1, :);     % Neumann BC

L = (2 * a / l)^2 * D2;         % Linear differential operator: u_t = L(u)
dL = dt * L;                    % Scaling for the finite difference loop

u_n = zeros(nx, nt);            % Numerical heat function

u_0 = A * sin((a^2 / l^2) * t); % Boundary condition

for i = 1 : nt - 1

u_n(1, i)     = u_0(i);      % Enforce the boundary condition
u_n(:, i + 1) = u_n(:, i) + dL * u_n(:, i);

end


The $$cheb()$$ function is from the book Spectral Methods in Matlab by Lloyd N. Trefethen. I didn't want to format that code so it would be the same as from the book:

% CHEB compute D = differentiation matrix, x = Chebyshev grid
function [D,x] = cheb(N)
if N==0, D=0; x=1; return, end
x = cos(pi*(0:N)/N)';
c = [2; ones(N-1,1); 2].*(-1).^(0:N)';
X = repmat(x,1,N+1);
dX = X-X';
D = (c*(1./c)')./(dX+(eye(N+1))); % off-diagonal entries
D = D - diag(sum(D')); % diagonal entries
end


Unfortunetly, my numerical solution differes from the analitical solution a lot. This can be seen on the picture below. The analitical solution is $$u_a$$, while the numerical is $$u_n$$.

My question is, where are the mistakes in my numerical solution? Can this problem be solved using Chebyshev's spectral method at all? Could the problem be that I don't have enough spatial points, so I am not catching the gradient correctly?

I listed my complete code below in case it is needed.

%% Heat equation

close all
clear
clc

%% Simulation time parameters

nt    = 2000001;                % Number of time points
t_max = 1;                      % Simulation time
t     = linspace(0, t_max, nt);
dt    = t(2) - t(1);            % Time step

%% Constants

a = 1;            % Heat constant
A = 100;          % Amplitude of the boundary condition signal
l = a / sqrt(pi); % Rod length

%% Numerical solution

N = 30;

[D, x] = cheb(N);     % D is the Chebyshev's differentiation matrix
nx     = length(x);   % Number of spatial points

D2 = D^2;
D2(N + 1, :) = D(N + 1, :);     % Neumann BC

L = (2 * a / l)^2 * D2;         % Linear operator: u_t = L(u)
dL = dt * L;                    % Scaling for the finite difference loop

u_n = zeros(nx, nt);            % Numerical heat function

u_0 = A * sin((a^2 / l^2) * t); % Boundary condition

for i = 1 : nt - 1

u_n(1, i)     = u_0(i);      % Enforce the boundary condition
u_n(:, i + 1) = u_n(:, i) + dL * u_n(:, i);

end

%% Analitical solution

u_a = zeros(nx, nt); % Analitical heat function

for n = 0 : 100

w_n = (0.5 + n) * pi;

u_a = u_a + sin(w_n * 0.5 * (1 - x)) * ( 4 * (sin(0.5 * w_n))^2 / (w_n + w_n^5) ) * ...
( w_n^2 * ( exp(-(w_n^2) * (a^2 / l^2) * t) - cos((a^2 / l^2) * t) ) - sin((a^2 / l^2) * t) );

end

u_a = A * ones(nx, 1) * sin((a^2 / l^2) * t) + A * u_a;

u_a(1, :) = A * sin((a^2 / l^2) * t);

%% Limits

u_a_max = max( max(u_a) );
u_a_min = min( min(u_a) );

u_n_max = max( max(u_n) );
u_n_min = min( min(u_n) );

u_max = max(u_a_max, u_n_max);
u_min = min(u_a_min, u_n_min);

%% Graphs

figure('Name', 'Heat function');

k = [(nt - 1) / 4, (nt - 1) / 2, (nt - 1) / 2 + (nt - 1) / 4, (nt - 1)];
k = k + 1;

for i = 1 : 4

subplot(2, 2, i);
plot(x, u_a(:, k(i)), 'r', x, u_n(:, k(i)), 'b', 'linewidth', 3);
set(gca,'XDir','reverse');
xlabel('x');
ylabel('u');
ylim([u_min, u_max]);
grid on
title(num2str(t(k(i)), 't : %.2f'));
legend('u_a', 'u_n');

end


• You have to set the Neumann boundary condition $$u_x(-1,x)=0$$ explicitly

for i = 1 : nt - 1
% Enforce Dirichlet BC
u_n(1, i)     = u_0(i);
% Enforce Neumann BC
u_n(end, i)   = u_n(end, i)-(D(end,:)*u_n(:,i))/D(end,end);
% Solve
u_n(:, i + 1) = u_n(:, i) + dL * u_n(:, i);
end


Solution:

• Thank you very much @ConvexHull, the solution really was simple! As for the direction of the coordinate system, It is a consequence of my coordinate transformation. I like to map $z = 0$ of my original system to $\theta = 0$ which means $\cos \theta = 1$. I also map $z = l$ to $\theta = \pi$. So if I revert to my original coordinate system, it goes from zero on the left to $l$ on the right. Dec 4, 2023 at 22:58