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?

$\endgroup$
  • 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 '15 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 '15 at 15:41
  • 1
    $\begingroup$ Piecewise should be capitalized, like all built-in Mathematica functions. $\endgroup$ – Michael Seifert Jul 1 '15 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 '15 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 '15 at 16:23
2
$\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, privacy policy and cookie policy

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