Solving a PDE with variable coefficients

Question

With the help from the expert Nasser I am able to write the code for a PDE with variable coefficients and mixed boundary conditions, defined on a square domain:

ClearAll[y, x1, x2];
a = x1 + x2;
pde = D[a D[y[x1, x2], x1], x1] + D[a D[y[x1, x2], x2], x2] - 4;
bc = {y[x1, 2] == 2 + x1, y[x1, 3] == 3 + x1};
sol = 
  NDSolve[
    {pde == NeumannValue[-1*a, x2 == 2] + NeumannValue[1*a, x2 == 3], bc}, 
    y, {x1, 2, 3}, {x2, 2, 3}]
Plot3D[Evaluate[y[x1, x2] /. sol], {x1, 2, 3}, {x2, 2, 3}, 
  PlotRange -> All, AxesLabel -> {"x1", "X2", "y[x1,x2]"}, BaseStyle -> 12]

where the Neumann boundary condition is defined for this problem as

$\qquad a(x)\,dy/dn,$

which is equal to

$\qquad a(x)\,(grad(y, x1)\,n1 + grad(y, x2)\,n2,$

which is equal, in this problem (using the known solution just to check the accuracy, to

$\qquad y = x1 + x2 = a(x)\,(n1 + n2),$

where $n1, n2$ are components of the normal vector.

Can we

1. Evaluate the numerical values of the solution on the domain?

2. Evaluate the relative error between the numerical and the exact solution which is equal to $y = x1 + x2$?

Answer 1

I am having hard time understanding your Neumman BC

why have a(x) on both sides, it then cancels? Also directional derivative is always normal to the surface in the outwards direction. So why write n1(x)+n2(x) there?

I took your Neumman as -1 on left edge and +1 on right edgs. But you can change that to zero as I do not understand your notation.

ClearAll[y, x1, x2];
a = x1 + x2;
pde = D[a D[y[x1, x2], x1], x1] + D[a D[y[x1, x2], x2], x2] - 4;
bc = {y[x1, 2] ==2+x1, y[x1, 3] == 3+x1};
sol = NDSolve[{pde == NeumannValue[-1, x2 == 2] + NeumannValue[1, x2 == 3]  , bc}, 
            y, {x1, 2, 3}, {x2, 2, 3}]

Plot3D[Evaluate[y[x1, x2] /. sol], {x1, 2, 3}, {x2, 2, 3}, 
 PlotRange -> All, AxesLabel -> {"x1", "X2", "y[x1,x2]"}, 
 BaseStyle -> 12]

If Neumann meant to be zero, then the result is

sol = NDSolve[{pde == NeumannValue[0, x2 == 2] + NeumannValue[0, x2 == 3]  , bc}, 
            y, {x1, 2, 3}, {x2, 2, 3}]

But feel free to adjust the above code if you think the Neumann should be something else.