Numerically solving a diffusion equation with a piecewise initial condition

Consider the diffusion equation given by with initial conditions and boundary conditions .

I wish to numerically solve this equation by the finite difference formula where with ∆x = 0.1 and ∆t = 0.01. The exact solution is given by

Here is the code that I have worked out so far.

clear all
D = 1/4;
dx = 0.1;
dt = 0.01;
x = 0:dx:1;
t = 0:dt:1;
Nx = length(x);
Nt = length(t);
u = zeros(Nx, Nt);
u(x<=1/2,1) = x(x<=1/2);
u(x>1/2,1) = 1 - x(x>1/2);
u(1,:) = 0;
u(Nx,:) = 0;

for j = 1:Nt-1
    for i = 2:Nx-1
        u(i,j+1) = D*(dt/dx^2)*(u(i+1,j) + u(i-1,j)) + (1-2*D*(dt/dx^2))*(u(i,j));
    end
end

% Compute the analytical solution V(x,t)
k = 1:100;
V = zeros(length(x), length(t));
for j = 1:Nt
    for i = 1:Nx
        V(i,j) = sum((4./(k.*pi).^2).*sin(k.*pi/2).* sin(k.*pi.*x(i)).*exp(-D.*t(j).*(k.*pi).^2));
    end
end

% Plot the numerical and exact solutions
figure;
[X,T] = meshgrid(x,t);
surf(X,T,u');
xlabel('x');
ylabel('t');
zlabel('u(x,t)');
title('Numerical solution');

figure;
surf(X,T,V');
xlabel('x');
ylabel('t');
zlabel('V(x,t)');
title('Exact solution');

However, there seems to be something wrong with the code. Any support would be greatly appreciated.

[Bill Greene]

Look at the pdepe function.

[Bas123]

Hi Bill, thanks for your message. I am supposed to use an explicit finite difference method for the problem.

[Alexander]

Hi! This was really helpful code for a HW assignment I'm working on. My graphs are a little wacky as well, I am using different differencing schemes, but I think the problem may lie in your parameters. I think by finding the correct dx and dt you may be able to generate smoother graphs.

However I am also pretty new to this and not sure exactly how to find the best spatial and time step sizes. Good luck!

Answer 1

[Torsten]

I changed

V(i,j) = sum((4./(k.*pi).^2).*sin(k.*pi/2).* sin(k.*pi.*x(j)).*exp(-0.5.*t(i).*(k.*pi).^2));

to

V(i,j) = sum((4./(k.*pi).^2).*sin(k.*pi/2).* sin(k.*pi.*x(i)).*exp(-D.*t(j).*(k.*pi).^2));

in your code and the results look at least similar.

You should check the results in detail, i.e. plots of u over x for some times t or plots of u over t for some positions x. A surface plot looks fine, but is not suited to find errors or unprecise results.

[Bas123]

Great, thanks a lot, as always!