Why is my boundary condition for the heat transfer problem not satisfied using the PDEPE function in MATLAB?

The left boundary of the heat transfer problem that I am trying to solve is not correct. I imposed a Neuman boundary condition on the left bound with "du/dx=0". The solution should look flat with a slope of zero on the left boundary but the solution plot shows a positive slope.

Answer 1

[MathWorks Support Team]

We have verified that there is a limitation in MATLAB 7.0 (R14)+ in the way that the PDEPE function handles boundary conditions with certain values.

PDEPE imposes Boundary Conditions in terms of flux "f()":

    p() + q()*f() = 0

At the left bound, "x=0". If the flux term is in a form such that it is zero for any specified boundary condition, the boundary condition will be ignored.

For example, "f() = (A+B*x)*du/dx" where "A" and "B" are constants. If "A = 0", flux is "(0+B*x)*du/dx".

You specify the boundary condition on the left bound to be "du/dx =0".

On the left bound, "B*(0)*du/dx = 0". This boundary condition is always satisfied regardless the value of "du/dx" and MATLAB will disregard the previously specified boundary condition.

Hence this particular combination:

     f() = x*du/dx

     0 <= x 

     Boundary Condition at 0,  du/dx = 0

gives an incorrect boundary result.

To work around this issue, shift the value of "x" from 0 to a very small number called delta that is close to 0. Define another point next to delta as follow:

x = [delta, 1.01*delta, linspace(2*delta,1000,100)];

The boundary condition on the left side should also be changed from:

pl = 0; 

ql = 1;

pr = 0; 

qr = 1;

to

pl = 0; 

ql = 1/(A+Bx*xl);

pr = 0; 

qr = 1;