Numerical solution for the heat equation does not match the exact solution

I write a matlab code for the problem, and i want to compare the numerucal and the exact solution at t=0.5 but the graph of the soultions are very different.

The problem is in the following:

% Parameters
L = 5; % Length of the rod
T = 1;  % Total time
Nx = 100; % Number of spatial points
Nt = 500; % Number of time steps
alpha = 0.01; % Thermal diffusivity

% Discretization
dx = 0.1; % Spatial step size
dt = 0.04; % Time step size
r = alpha * dt / dx^2; % Stability parameter for the scheme


% Initial condition
x = linspace(0, L, Nx);% generates Nx points. The spacing between the points is (L-0)/(Nx-1).
u = cos(x-0.5); %  initial condition



% Preallocate solution matrix
U = zeros(Nx, Nt+1); %  returns an Nx-by-Nt matrix of zeros
U(:,1) = u; %extracts all the elements from the first column of matrix

% Time-stepping loop
for n = 1:Nt
    % Update the solution using explicit scheme for the Heat Equation
    for i = 2:Nx-1
        U(i,n+1) = U(i,n) + r * (U(i+1,n) - 2*U(i,n) + U(i-1,n));
    end
end

Numerical_sol = U


% Plot results
t = linspace(0, T, Nt+1);
[X, T] = meshgrid(x, t);

% Exact solution
Exact= @(x, t) exp(-t) .* cos(pi * x - 0.5 * pi);


figure;
mesh(X, T, U');
xlabel('x');
ylabel('t');
zlabel('u(x,t)');
title('Heat Equation Solution');


figure
plot(x, U(:,1), 'o');

[John D'Errico]

You accepted the answer to this question. Then you asked it just again, 42 minutes ago. I closed the duplicate question, since it asks nothing you should not have learned here.

[Torsten]

@John D'Errico

I didn't know that you are the Website Master taking the right to delete questions and answers that are actually in flow.

Answer 1

[Torsten]

I cannot believe that the value for alpha has no influence on the solution.

So I guess that the exact solution maybe holds for alpha = 1.

[Erm]

if i set alpha =1 there is no results for any graphs !

[Torsten]

Why don't you use the values for L, dx and dt from your assignment ?

% Parameters
L = 1; % Length of the rod
T = 1;  % Total time
% Discretization
dx = 0.1; % Spatial step size
dt = 0.04; % Time step size
Nx = L/dx + 1; % Number of spatial points
Nt = T/dt; % Number of time steps
alpha = 1; % Thermal diffusivity
r = alpha * dt / dx^2; % Stability parameter for the scheme
% Initial condition
x = linspace(0, L, Nx);% generates Nx points. The spacing between the points is (L-0)/(Nx-1).
u = cos(pi*(x-0.5)); %  initial condition
% Preallocate solution matrix
U = zeros(Nx, Nt+1); %  returns an Nx-by-Nt matrix of zeros
U(:,1) = u; %extracts all the elements from the first column of matrix
% Time-stepping loop
for n = 1:Nt
    % Update the solution using explicit scheme for the Heat Equation
    for i = 2:Nx-1
        U(i,n+1) = U(i,n) + r * (U(i+1,n) - 2*U(i,n) + U(i-1,n));
    end
end
plot(x,[U(:,1),U(:,2),U(:,3),U(:,4),U(:,5)])

[Erm]

somtimes I change the interval up to 5 to see the bahavior of the graph. however , how we compare the exact solution with the numerical at specific time since they are differents.

[Torsten]

The solution you gave is not the exact solution.

Just differentiate exp(-t)*cos(pi*(x-0.5)) with respect to t and two times with respect to x and see that the results are different.

[Erm]

[Erm]

% Parameters
L = 1; % Length of the rod
T = 1;  % Total time
Nx = 100; % Number of spatial points
Nt = 500; % Number of time steps
alpha = 0.01; % Thermal diffusivity

% Discretization
dx = 0.1; % Spatial step size
dt = 0.04; % Time step size
r = alpha * dt / dx^2; % Stability parameter for the scheme


% Initial condition
x = linspace(0, L, Nx);% generates Nx points. The spacing between the points is (L-0)/(Nx-1).
u = cos(x-0.5); %  initial condition



% Preallocate solution matrix
U = zeros(Nx, Nt+1); %  returns an Nx-by-Nt matrix of zeros
U(:,1) = u; %extracts all the elements from the first column of matrix

% Time-stepping loop
for n = 1:Nt
    % Update the solution using explicit scheme for the Heat Equation
    for i = 2:Nx-1
        U(i,n+1) = U(i,n) + r * (U(i+1,n) - 2*U(i,n) + U(i-1,n));
    end
end

Numerical_sol = U


% Plot results
t = linspace(0, T, Nt+1);
[X, T] = meshgrid(x, t);

% Exact solution
Exact= @(x, t) exp(-t) .* cos(pi * x - 0.5 * pi);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Define the range for x
X = linspace(0, 1, 100); % 100 points between 0 and 1

% Define a specific value for t
t = 1; % You can change this value to see different plots

% Compute the function values
f = exp(-t) .* cos(pi * X - 0.5 * pi);



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

figure;
mesh(X, T, U');
xlabel('x');
ylabel('t');
zlabel('u(x,t)');
title('Heat Equation Solution');


figure
plot(x, U(:,1), 'o');
hold on 

plot(X, f, 'LineWidth', 2);
xlabel('x');
ylabel('f(x)');
title(['Plot of f(x) = exp(-t) .* cos(\pi x - 0.5 \pi) for t = ' num2str(t)]);
Hold off

[Erm]

above i update the code to plote the exact solution

[Torsten]

The blue dotted curve cannot be correct because u(0,t) = u(1,t) = 0 must hold for all times t.

And as said:

The solution you call "exact" is not a solution of the differential equation. Thus comparing with the numerical solution makes no sense.

[Torsten]

The exact solution for the differential equation

du/dt = alpha * d^2u/dx^2

with initial condition

u(0,x) = cos(pi*(x-0.5))

and boundary conditions

u(t,0) = u(t,1) = 0

is given by

uexact(t,x) = exp(-alpha*pi^2*t) * cos(pi*(x-0.5))

That's the function you have to compare the numerical solution with.

syms x t alpha pi
u = exp(-alpha*pi^2*t) * cos(pi*(x-0.5));
diff(u,t)

ans = 

alpha*diff(u,x,2)

ans = 

[Erm]

Thanks a lot