# NDsolve with discontinuous non-numeric PDE coefficients

I want to solve the heat diffusion equation with variable thermal conductivity. So I have the following code

Clear["Global'*"];
Needs["NDSolveFEM"];
mesh = ToElementMesh[Rectangle[{0, 0}, {10, 10}]];
k[x_, y_] := piecewise[{{10, x^2 + y^2 <= 2}, {5, x^2 + y^2 > 2}}];

gT = Grad[w[x, y], {x, y}];
tC = {{k, 0}, {0, k}};
pde = Div[tC.gT, {x, y}];

sol = First@NDSolve[{pde == 0,
DirichletCondition[w[x, y] == 100, y == 0],
DirichletCondition[w[x, y] == 400, y == 10],
DirichletCondition[w[x, y] == 0, x == 0],
DirichletCondition[w[x, y] == 0, x == 10]},
{w}, {x, y} \[Element] mesh];

ContourPlot[sol[x, y], {x, 0, 10}, {y, 0, 10}]


However I get this error

The PDE coefficient {{k,0},{0,k}} does not evaluate to a numeric \
matrix of dimensions {2,2}.


And I have not been able to figure out why this occurs and I've already looked at the following which didn't help.

Any help would be appreciated. Thanks.

Also, as another problem I'm experiencing, is that if I take away the First@ in First@NDsolve. I get the following errors.

0.000714286 cannot be used as a variable
0.000714286 cannot be used as a variable


Why does this occur?

• k and k[x,y] are different things as far as Mathematica is concerned. Have you tried replacing the former with the latter in your definition for tC? Jul 1, 2015 at 15:38
• If I do that I get The PDE coefficient {{piecewise[{{10,x^2+y^2<=2},{5,x^2+y^2>2}}],0},{0,\ piecewise[{{10,x^2+y^2<=2},{5,x^2+y^2>2}}]}} does not evaluate to a \ numeric matrix of dimensions {2,2}. and 0.0007142857142857143 is not a valid variable Jul 1, 2015 at 15:41
• Piecewise should be capitalized, like all built-in Mathematica functions. Jul 1, 2015 at 15:45
• There's one more error in your code (it's not an uncommon error for Mathematica newbies, and completely understandable.) See my full answer below. Jul 1, 2015 at 16:01
• @MonkSphere The plot command should be ContourPlot[w[x, y] /. sol // Evaluate, {x, 0, 10}, {y, 0, 10}]. The solution sol is a replacement rule for w. (Type sol to see.) BTW, you can set up the mesh with an interior boundary that matches the discontinuity in k. See the "Element Mesh Creation" tutorial. (In this example, the discontinuity is practically insignificant, though.) Jul 1, 2015 at 16:23

As noted in the comments, there are three main errors:

1. k and k[x, y] are not the same thing in Mathematica.
2. Piecewise should be capitalized.
3. NDSolve returns a list of "rules" for the various solutions of the equations, which need to be "applied" (using /.) to be plotted. Alternately, if you know that there's only going to be one solution of the equations, you can use NDSolveValue instead of NDSolve, which returns the function itself rather than the rules.

Corrected code and results:

Clear["Global'*"];
Needs["NDSolveFEM"];
mesh = ToElementMesh[Rectangle[{0, 0}, {10, 10}]];
k[x_, y_] := Piecewise[{{10, x^2 + y^2 <= 2}, {5, x^2 + y^2 > 2}}];

gT = Grad[w[x, y], {x, y}];
tC = {{k[x, y], 0}, {0, k[x, y]}};
pde = Div[tC.gT, {x, y}];

sol = First@NDSolve[{pde == 0,
DirichletCondition[w[x, y] == 100, y == 0],
DirichletCondition[w[x, y] == 400, y == 10],
DirichletCondition[w[x, y] == 0, x == 0],
DirichletCondition[w[x, y] == 0, x == 10]},
{w}, {x, y} \[Element] mesh];

ContourPlot[w[x, y] /. sol, {x, 0, 10}, {y, 0, 10}]


Alternately,

sol = First@NDSolveValue[{pde == 0, ...

ContourPlot[sol[x, y], {x, 0, 10}, {y, 0, 10}]


gives the same results.