# MATLAB: incorrect computation with pdepe for parabolic-elliptic system

I have found that MATLAB solves 1D parabolic-elliptic system incorrectly by using pdepe function. Here is a system: $$u_t = u_{xx} + 2,$$ $$0 = v_{xx} + u.$$ Boundary conditions: $$(-u_x+u)|_{x=0} = 2t, \;\; (u_x+u)|_{x=1} = 2t,$$ $$(-v_x+v)|_{x=0} = 0, \;\; (v_x+v)|_{x=1} = 0.$$ Initial conditions: $$u|_{t=0} = 0,\;\;v|_{t=0} = 0.$$ MATLAB documentation informs that pdepe can solve such systems.

Source code:

function a

tt = 1;

m = 0;
x = linspace(0,1,100);
t = linspace(0,tt,100);

sol = pdepe(m,@apde,@aic,@abc,x,t);
u1 = sol(:,:,1);
u2 = sol(:,:,2);

figure;
surf(x,t,u1);
title('u1(x,t)');
xlabel('Distance x');
ylabel('Time t');
colorbar
set(gca,'YDir','reverse');

figure;
surf(x,t,u2);
title('u2(x,t)');
xlabel('Distance x');
ylabel('Time t');
colorbar
set(gca,'YDir','reverse');

figure
plot(x,u1(end,:))
title('u1 at t = tt')
xlabel('Distance x')
ylabel('u1(x,tt)')

figure
plot(x,u2(end,:))
title('u2 at t = tt')
xlabel('Distance x')
ylabel('u2(x,tt)')

figure
plot(x,tt*(1+x-x.^2))
title('TRUE u2 at t = tt')
xlabel('Distance x')
ylabel('u2(x,tt)')

% --------------------------------------------------------------

function [c,f,s] = apde(x,t,u,DuDx)
c = [1; 0];
f = [1; 1] .* DuDx;
s = [2; u(1)];

% --------------------------------------------------------------

function u0 = aic(x)
u0 = [0; 0];

% --------------------------------------------------------------

function [pl,ql,pr,qr] = abc(xl,ul,xr,ur,t)
pl = [ul(1)-2*t; ul(2)];
ql = [-1; -1];
pr = [ur(1)-2*t; ur(2)];
qr = [1; 1];


Exact solution: $$u(x,t) = 2t,\;\;v(x,t) = t(1+x-x^2).$$

MATLAB computes $u$ correctly. Exact solution for $v$ at $t = 1$:

Approximate solution computed by MATLAB:

This function doesn't even satisfy the boundary conditions.

$v(x,t)$ graph:

Let us replace our elliptic equation with a parabolic one, that is use the vector

c = [1; 1e-100];


c = [1; 0];


It means that we use the equation $$10^{-100}v_t = v_{xx} + u$$

Now MATLAB computed a correct solution:

Thus, this example demonstrates that MATLAB solves parabolic-elliptic systems with Robin boundary conditions incorrectly.

If we use Dirichlet boundary conditions, the solution is correct.

Is my reasoning correct?

The latest version Mathematica solves this problem very well:

Clear[u,v,nx,eqn1,eqn2,ic1,ic2,vars,eqn1,eqn2,ic1,ic2,vars,sol,PDEsol,order,npts];
npts = 100; order=12;
nx = Range[0,1,1/(npts-1)];
{ddx,d2dx2} = Map[NDSolveFiniteDifferenceDerivative[Derivative[#],nx, "DifferenceOrder" -> order]["DifferentiationMatrix"]&,{1,2}] ;
u = Array[u$[#][t]&,npts]; v = Array[v$[#][t]&,npts];

eqn2 = Thread[d2dx2.v+ u==0];

eqn1[[1]]=u[[1]]==ddx[[1]].u+2t;
eqn1[[-1]]=u[[-1]]==2t-ddx[[-1]].u;
eqn2[[1]]=v[[1]]==ddx[[1]].v;
eqn2[[-1]]=v[[-1]]==-ddx[[-1]].v;

vars = Flatten[{u, v}] /. x__[t] :> x ;
PDEsol = First[NDSolve[{eqn1, eqn2, ic1, ic2}, vars, {t, 0, 30}]];
time = Flatten@First[vars[[1]] /. PDEsol];
{usol, vsol} =
Flatten[Table[
Transpose[{nx, ConstantArray[ti, npts], # /. PDEsol}] /.
t -> ti, {ti, time[[1]], time[[2]], 0.1}], 1] & /@ {u, v};

GraphicsRow[ ListPlot3D[#, PlotRange -> All, PlotTheme -> "Classic",
Mesh -> {50, 50}, Lighting -> "Neutral"] & /@ {usol, vsol},
ImageSize -> 800]


• I wrote almost the same code earlier. You've read my mind! :) – jokersobak Feb 17 '16 at 11:10
• My post here was inspired by and plagiarized` from the same Wolfram tutorial document we read. – LCFactorization Feb 17 '16 at 12:34

No. You have not changed the boundary conditions, the boundary conditions are still Robin type. what you changed is the DEs to both parabolic and used the same boundary conditions. It is in fact surprising to see this. Reasoning is not correct though.

• Why can't boundary conditions be of Robin type in DAE system? Moreover, the boundary and initial conditions are consistent in this problem and BCs for $v$ don't depend on $t$. – jokersobak Feb 17 '16 at 11:04