2
$\begingroup$

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:

Analitical solution check

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$.

Comparison

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
$\endgroup$

1 Answer 1

3
$\begingroup$

The answer is quite simple:

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

Add following line (fifth line):

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:

enter image description here

$\endgroup$
1
  • $\begingroup$ 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. $\endgroup$ Dec 4, 2023 at 22:58

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.