# Question on how MATLAB's pdepe solver works

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?

• Can you put axes on your plots and maybe explain a bit more? Your numerical problem is definined on $[0,10]$ but your PDE has a left inflow condition at $x=1$? Also, you are taking $\Delta x$ quite large. Taking it smaller and smaller and still seeing nondecreasing error is a good indication that your spatial scheme is messed up somewhere, and the $x$ location where the errors are high might also tell you if one of the boundary conditions is messed up Feb 9, 2020 at 0:45
• @whpowell96 Thanks a lot for the response. When you say ,"but your PDE has a left inflow condition at x=1?" do you mean there is a mistake in this line dC(1,1) = 0? I set dC(1,1)=0 since MATLAB indexing starts from 1. Please let me know if you are point to any other line in the code Feb 10, 2020 at 2:00

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

• Thanks a ton for the detailed solution. From what I understand, the truncation error of the backward difference is O( $\Delta x$) and that of centered difference is O($\Delta x^2$). The errors obtained from the numerical scheme is of the order 0.04 for backward + central difference and 1e-13 for central alone. I am not able to clearly understand what leads to this drastic difference. Could you please elaborate a bit on this? Feb 10, 2020 at 16:34
• That is really a separate question about the accuracy of these different numerical schemes and is covered in detail in many texts. Your mesh in the spatial direction is quite coarse. If you use one or two orders of magnitude more points, you can also get solutions with the two schemes that are close to each other. But your question basically was "How does pdepe work?" and that is what I answered. Feb 10, 2020 at 16:46
• Thanks a lot for the response. I've created a new post here Feb 11, 2020 at 16:44