2
$\begingroup$

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$: Exact solution for $v$ at $t = 1$

Approximate solution computed by MATLAB: Approximate solution computed by MATLAB

This function doesn't even satisfy the boundary conditions.

$v(x,t)$ graph:

$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:

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?

$\endgroup$
1
$\begingroup$

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]

enter image description here

| cite | improve this answer | |
$\endgroup$
  • 1
    $\begingroup$ I wrote almost the same code earlier. You've read my mind! :) $\endgroup$ – jokersobak Feb 17 '16 at 11:10
  • $\begingroup$ My post here was inspired by and plagiarized from the same Wolfram tutorial document we read. $\endgroup$ – LCFactorization Feb 17 '16 at 12:34
0
$\begingroup$

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.

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ 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$. $\endgroup$ – jokersobak Feb 17 '16 at 11:04

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.