Question on how MATLAB's pdepe solver works

Question

I'm solving the following 1D transport equation in MATLAB's pdepe solver.

$$\frac{\partial C}{\partial t} = D\frac{\partial ^2 C}{\partial x^2}-v\frac{\partial C}{\partial x}$$

At the inlet (left boundary), Dirichlet boundary condition is applied $C(1) = C_L$ . (1 is the inlet node number)

At the outlet (right boundary), diffusive flux is ignored. $-D \frac{dC}{dx} = 0$

Implementation of the above boundary conditions in MATLAB's pdepe solver.

function [pl,ql,pr,qr] = bcfun(xl,cl,xr,cr,t)
    pl = cl - 3;
    ql = 0;
    pr = 0;
    qr = 1;
end

From what I understand, the spatial direction is discretized and the resulting ode's are solved using ode15s in pdepe solver.

I'd tried to do the same in my own version of code that implements whats is done in pdepe solver. However, my results don't agree with the pdepe solver. I've used the backward difference scheme for the first derivative and centered difference scheme for the second derivative.I'm not sure about the scheme implemented in MATLAB's pdepe solver.

I've implemented the boundary conditions in the following way.

dC(1) = 0
dC(nnode,1) = -v*(C(nnode) - C(nnode-1))/delx + (D/delx^2)*2*(C(nnode-1) - C(nnode))

Right boundary condition: $-D \frac{dC}{dx} = 0$ $$ \frac{C_{N+1} - C_{N-1}}{2 \Delta x} = 0$$

At last node,

dC(nnode,1) = -v*(C(nnode) - C(nnode-1))/delx + (D/delx^2)*(C(nnode-1) - C(nnode) +C(nnode+1))

is equal to

dC(nnode,1) = -v*(C(nnode) - C(nnode-1))/delx + (D/delx^2)*2*(C(nnode-1) - C(nnode))

The complete code is

function sol=so()
format short
global D nnode init_co find_index v
m = 0;
delx = 0.25;
xend = 10; 
D = 500;
v = 200;
x = 0:delx:xend;
find_index  = x;
tspan = 0:0.00001:1;
init_co = [3 ; zeros(length(x)-1,1)];
nnode = length(x);

%% pdepe solver
sol = pdepe(m,@pdefun,@icfun,@bcfun,x,tspan)
figure(1)
subplot(1,2,2)
plot(tspan,sol)
xlabel('time')
ylabel('c')
xlim([-0.01 0.5])
ylim([2.995 3.005])
title('MATLAB - pdepe')
grid on


function [g,f,s] = pdefun(x,t,c,DcDx)
g = 1;
f = D*DcDx;
s = -v*DcDx;
end

function c0 = icfun(x)
c0 = init_co(find(find_index==x));
end

function [pl,ql,pr,qr] = bcfun(xl,cl,xr,cr,t)
    pl = cl - 3;
    ql = 0;
    pr = 0;
    qr = 1;
end

function dC= mysolver(t,C)
    dC(1,1) = 0;
    for i = 2:nnode-1
        dC(i,1) = -v*(C(i) - C(i-1))/delx + D/delx^2*(C(i-1) -2*C(i) + C(i+1)); 
    end
    dC(nnode,1) = -v*(C(nnode) - C(nnode-1))/delx + (D/delx^2)*2*(C(nnode-1) - C(nnode)); % DdC/dx = 0
end

%% my solver
[tspan C]  = ode15s(@(t,s) mysolver(t,s), tspan , init_co);
figure(1)
subplot(1,2,2)
plot(tspan,C)
xlabel('time')
ylabel('c')
xlim([-0.01 0.5])
ylim([2.995 3.005])
title('My solver')
grid on


% figure(2)
% plot(tspan, abs(sol - C))
% title('Absolute error')
end

The resulting absolute error( pdepe solution - my implementation) is

Also, the absolute error increases with an increase in mesh size (increasing delx from 0.25 to 1).

I am not sure why the absolute error increases. Is it because of the backward and centered difference schemes that I 've used or is it due to the way in which my boundary condition has been implemented?

Any suggestions?

Answer 1

The main difference between pdepe and your finite difference code is that pdepe basically uses a central difference approximation consistently where your code uses a combination of backward and central difference approximations.

When I run my modified code, shown below, with your original spatial discretization, the maximum difference between the pdepe solution and finite difference is on the order of 1e-13.

You'll notice I made a few other changes to your code. One of these was to reduce the time span of the solution; very little happens to the solution after about .1 seconds. Also, I changed the way the Dirichlet constraint is applied at the left end. You used a "rate form" to specify this constraint. Theoretically that is fine but numerically it introduces a slight error in the solution; this is known as "constraint drift." I prescribed this constraint using an algebraic equation; this is also how pdepe applies this type of constraint. If you run my code with the rate form of the constraint by setting useRateFormDirichlet=true, the maximum difference is 1e-10; so the effect of this change is small.

function cse_02_09_20
m = 0;
delx = 0.25;
xend = 10; 
D = 500;
v = 200;
x = 0:delx:xend;
find_index  = x;
tf=.1;
tspan=linspace(0,tf,100);
init_co = [3 ; zeros(length(x)-1,1)];
nnode = length(x);
fdRHS = @(t,x) mysolver(t,x,v,D,delx);
useRateFormDirichlet=false;
fdRHS = @(t,x) cdRHS(t,x,v,D,delx,useRateFormDirichlet);
% make ode solver tolerances very small so we can
% better see effects of spatial discretoization differences
opts=odeset('abstol', 1e-10, 'reltol', 1e-9);
%% pdepe solver
sol = pdepe(m,@pdefun,@icfun,@bcfun,x,tspan,opts);

  function [g,f,s] = pdefun(x,t,c,DcDx)
    nx=length(x);
    g = ones(1,nx); 
    f = D*DcDx;
    s = -v*DcDx;
  end

  function c0 = icfun(x)
    c0 = init_co(find(find_index==x));
  end

  function [pl,ql,pr,qr] = bcfun(xl,cl,xr,cr,t)
    pl = cl - 3;
    ql = 0;
    pr = 0;
    qr = 1;
  end

%% finite difference solver
e = ones(nnode,1);
jpat = spdiags([e e e],-1:1,nnode,nnode);
% improve performance by specifying a pattern for the jacobian
opts=odeset(opts, 'jpattern', jpat);
if ~useRateFormDirichlet
  opts=odeset(opts, 'mass', spdiags([0; ones(nnode-1,1)], 0, nnode, nnode));
end
tic
[tspan, C]  = ode15s(fdRHS, tspan , init_co, opts);
toc

solutionDifference=abs(sol-C);

figure; plot(tspan, sol(:,end), tspan, C(:,end)); grid;
title 'end C as a function of time'
legend('pdepe', 'finite difference');

figure; plot(x, sol(end,:), x, C(end,:)); grid;
title 'C at final time';
legend('pdepe', 'finite difference');

figure; plot(tspan, solutionDifference(:,end)); grid;
title 'tip difference as a function of time'

maxSolDiff=max(solutionDifference(:));
fprintf('Maximum difference between pdepe and finite difference=%g\n', ...
  maxSolDiff);

end

function dC= mysolver(t,C,v,D,delx)
N=size(C,1);
dC=zeros(N,1);
i = 2:N-1;
dC(i) = -v*(C(i) - C(i-1))/delx + D/delx^2*(C(i-1) -2*C(i) + C(i+1));
dC(N) = -v*(C(N) - C(N-1))/delx + (D/delx^2)*2*(C(N-1) - C(N)); % DdC/dx = 0
end

function dC=cdRHS(t,C,v,D,delx,useRateFormDirichlet)
N=size(C,1);
dC=zeros(N,1);
if ~useRateFormDirichlet
  dC(1)=C(1)-3;
end
i = 2:N-1;
dC(i) = -v*(C(i+1) - C(i-1))/(2*delx) + D/delx^2*(C(i-1) -2*C(i) + C(i+1));
dC(N) = -v*(C(N) - C(N-1))/delx + 2*D/delx^2*(C(N-1) - C(N)); % DdC/dx = 0
end