I've been working on trying to analyze the Heat Equation in water both experimentally and theoretically.
The model goes as: there's a cuboidal bath (of say, 15x7x5 inches) filled with water, and an immersion rod is placed at one corner of it. I haven't quite reached the part of placing the boundary condition by approximating the rod as a cylindrical and have in this code, assumed a segment of the wall (3x3) is being heated at 200 degrees Celsius. The water is around 25 degrees. While the experimental part is done, I'm currently having trouble executing it theoretically. I was recommended to use Mathematica to solve the PDE and plot it (and hopefully animate it) over time, but I am stuck. Since it is my first time using it, I've borrowed and assimilated my code from various other problems, and as such, I don't know which part is wrong.

b = 1.43*0.0000001
 f = NDSolve[{D[T[x, y, z, t], t] == 
     b*(D[T[x, y, z, t], x, x] + D[T[x, y, z, t], y, y] + 
        D[T[x, y, z, t], z, z]), 
    T[x, y, z, 0] == If [0 < y < 3 && 0 < z < 3, 200, 25], 
    T[x, y, z, 0] == T[0, y, z, 0] == T[x, 0, z, 0] == 
     T[x, y, 0, 0] == T[15, y, z, 0] == 25, {x, 0, 15}, {y, 0, 7}, {z,
      0, 5}, {t, 0, 20}]; T
c = Table[
  Plot3D[T[x, y, z, t] /. f, {x, 0, 15}, {y, 0, 7}, {z, 0, 5}, 
   Mesh -> 15, PlotRange -> {{0, 5}, {0, 5}, {0.5}}, 
   ColorFunction -> Function[{x, y, z}, Hue[.3 (1 - x)]]], {3, 4, 2, 
   t}]; c

It doesn't display any error but has been running for almost an hour. What I am hoping to achieve here is the evolution of heat in 3D space, and hopefully, animate it later on. Any help on how to set an arbitrary boundary of a cylinder floating in space would also be welcome. Thanks!


1 Answer 1


There are several issue with your code, here is a version that roughly does what you want, though you still need to think about the scale of things.

b = 1.43;

Using a smaller than default mesh makes things go faster, but they are less accurate. You can play with this once the code does what you want.

region = ToElementMesh[Cuboid[{0, 0, 0}, {15, 7, 20}], 
   "MaxCellMeasure" -> 1];

I added a Monitor, then you know where in the solution process you are. I also fixed a syntax error in your posted code and re-wrote the Dirichlet boundary conditions as DirichletCondition - you can not have something like T[0,y,z,y]==something also you have several definitions for T[x,y,z,0].

Monitor[solution = 
   NDSolveValue[{D[T[x, y, z, t], t] == 
      b*Laplacian[T[x, y, z, t], {x, y, z}], 
     T[x, y, z, 0] == If[0 <= y <= 3 && 0 <= z <= 3, 200, 25],
      T[x, y, z, t] == 25, {x == 0 || x == 15 || z == 0 || y == 0}]}, 
    T, Element[{x, y, z}, region], {t, 0, 20}, 
    EvaluationMonitor :> (monitor = Row[{"t = ", CForm[t]}])], 

Visualize the solution at t=20. First get this right before you think about animations.

SliceContourPlot3D[solution[x, y, z, 20], Element[{x, y, z}, region]]

enter image description here


Your Answer

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

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