0
$\begingroup$

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

Clear["Global'*"];
Needs["NDSolve`FEM`"];
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.

http://reference.wolfram.com/language/FEMDocumentation/tutorial/SolvingPDEwithFEM.html NDSolve for axisymmetric problem

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?

Improve this question
$\endgroup$
5
  • 1
    $\begingroup$ 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? $\endgroup$
    – Michael Seifert
    Jul 1, 2015 at 15:38
  • $\begingroup$ 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 $\endgroup$
    – MonkSphere
    Jul 1, 2015 at 15:41
  • 1
    $\begingroup$ Piecewise should be capitalized, like all built-in Mathematica functions. $\endgroup$
    – Michael Seifert
    Jul 1, 2015 at 15:45
  • $\begingroup$ 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. $\endgroup$
    – Michael Seifert
    Jul 1, 2015 at 16:01
  • $\begingroup$ @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.) $\endgroup$
    – Michael E2
    Jul 1, 2015 at 16:23

1 Answer 1

Reset to default
3
$\begingroup$

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["NDSolve`FEM`"];
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}]

enter image description here

Alternately,

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

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

gives the same results.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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