MATLAB: incorrect computation with pdepe for parabolic-elliptic system

Question

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');
shading interp
colorbar
set(gca,'YDir','reverse');

figure;
surf(x,t,u2);
title('u2(x,t)');
xlabel('Distance x');
ylabel('Time t');
shading interp
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];

instead of

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?

Answer 1

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[NDSolve`FiniteDifferenceDerivative[Derivative[#],nx, "DifferenceOrder" -> order]["DifferentiationMatrix"]&,{1,2}] ;
u = Array[u$[#][t]&,npts]; v = Array[v$[#][t]&,npts];

eqn1=Thread[D[u,t]==d2dx2.u+2];
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;

ic1=Thread[u==0]/.t->0;
ic2=Thread[v==0]/.t->0;

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]

Answer 2

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.