I would like to directly extract the surface plot data from the solution to a differential equation (using NDSolve). While there are many examples of this I haven't been able to find an example in 2D.

Here I am solving a 1D diffusive heat equation and plotting the results (the solution is very poorly behaved with significant artifacts). I would like to be able to extract the heat distributions along the x axis (i.e. the temperature distribution in the rod) for each moment in time and save to something like a text file(s) so I can have direct access to the simulation data and analyze it using a different program. This is easiest done with InterpolatingFunctionAnatomy but I am not sure how to extract and package it.

If anyone has suggestions for improving the accuracy of the solution I would also appreciate it. (I know I can decrease MaxCellMeasure but it greatly increases computation time).

parameters = {v -> -1, DT -> 0.000000143, rho -> 1000, c -> 4200, 
eta -> 0.001002};
HEATIMP = (-D[u[r, t], t] - 
 v D[u[r, t], r] + (DT/r^2) D[r^2. D[u[r, t], r], 
   r] + (12. eta)/(rho c) (v/r) (v/r)) /. parameters;
BCE1 = DirichletCondition[u[r, t] == 1, r == 2];
BCE2 = DirichletCondition[u[r, t] == 0, r == 1];
BCTot = {BCE1, BCE2};
ic = {u[r, 0] == 0};
heatdist = 
NDSolve[{HEATIMP == 0, BCTot, ic}, u, {r, 1, 2}, {t, 0.00000, 1}, 
Method -> {"FiniteElement", 
"MeshOptions" -> {MaxCellMeasure -> 0.005}}]
Plot3D[Evaluate[First[u[r, t] /. %]], {r, 1, 2}, {t, 0.00000, 1}, 
 MaxRecursion -> 6, WorkingPrecision -> 20,  
PlotRange -> Full, AxesLabel -> {RodPosition, Time, Temperature}]

Very Poor Solution

  • $\begingroup$ You are not solving a 1D diffusion equation but a 2D convection diffusion equation. See my answer. $\endgroup$
    – user21
    Jan 8, 2020 at 8:15

2 Answers 2


Decreasing the MaxCellMeasure to $5\times10^{-5}$ improves the quality of the solution considerably while still taking only a few seconds to calculate, so I would suggest that as a viable compromise:

heatdist = 
  NDSolve[{HEATIMP == 0, BCTot, ic}, u, {r, 1, 2}, {t, 0.00000, 1}, 
    Method -> {"FiniteElement", "MeshOptions" -> {MaxCellMeasure -> 5*^-5}}

  Evaluate[u[r, t] /. First@heatdist], {r, 1, 2}, {t, 0.00000, 1},
  PlotRange -> All, AxesLabel -> {RodPosition, Time, Temperature}

plot of improved solution

Note that I've removed MaxRecursion -> 6 and WorkingPrecision -> 20 from your plotting code because the problems were with the solution from NDSolve; more recursion or precision in plotting could not improve the solution, but would only slow the plot down.

Rather than digging within the structure of the interpolating function, would it be acceptable to simply evaluate the solution from NDSolve at a range of times, and for each time at a range of positions across the rod (i.e. along the rod's profile)?

Something like this (using your definitions):

u[r, t] /. First@heatdist;
data = Table[{r, t, u[r, t] /. First@heatdist}, {t, 0.00000, 1, 0.05}, {r, 1, 2, 0.01}];

As defined above, data contains 21 lists, each of 101 triplets; each triplet is structured as {position along the rod, time, temperature}. Each of the 21 lists corresponds to a point in time; each of the 101 triplets corresponds to a position along the rod.

As a sample of what data would contain, here are a series of temperature profiles across the rod over time:

ListLinePlot[#, PlotRange -> All] & /@ data[[All, All, {1, 3}]]

several temperature profiles at different times

You could then export data to a file, or modify the Table to generate appropriate subsets of interest and then export each one, etc.

  • $\begingroup$ Thanks, that is exactly what I was looking for. I still have several questions comments. I put MaxRecursion -> 6 in there because for other solutions I was getting jagged edges and other artifacts and doing this improved the solution (albeit at a significant performance penalty). I had changed MaxCellMeasure but even with smaller cells there are still rather noticeable issues with numerical accuracy (specifically the overshoot above a value of 1). Last question: How would I go about exporting only the height values (third element) of the table? $\endgroup$
    – I am tired
    Jan 7, 2020 at 18:46
  • $\begingroup$ @Iamtired 1) MaxRecursion allows Plot to more faithfully reproduce any fine structure that is present in the function being plotted, but does not "smooth" the function in any way; that has to be done in NDSolve, practically most easily by reducing the size of the FEM cells. 2) Regarding the overshoot above a value of 1, that is a rather small artifact, but again you would want to look into the suggestions provided by NDSolve in relation with the NDSolve::femcscd warning that you get. $\endgroup$
    – MarcoB
    Jan 7, 2020 at 19:23
  • $\begingroup$ @Iamtired 3) To export only the temperature, modify the table to include only that? E.g. data = Table[u[r, t] /. First@heatdist}, {t, 0.00000, 1, 0.05}, {r, 1, 2, 0.01}]? and then use Export appropriately. This depends on how you want your data exported, e.g. do you want one file per time step? $\endgroup$
    – MarcoB
    Jan 7, 2020 at 19:25
  • $\begingroup$ I would like for each row in the file to correspond to each plotted temperature distribution. For the plots, it would look something like this (plots resemble a step function) 0 0 ... 0 0 1 // 0 0 ... 0 1 1 // 0 0 ... 1 1 1 // as the system warms up. I would guess this involves using "Flatten" somehow to remove the innermost bracer and then plotting each list of numbers as its own row. $\endgroup$
    – I am tired
    Jan 7, 2020 at 19:47
  • $\begingroup$ @Iamtired With the definitions above, Export["yourfilename.txt", Chop@data[[All, All, 3]], "Table"] should give what you want, if I understood your requirements. The "Table" format will give you a tab-separated file; if you want a comma-separated one, you could use the "CSV" format instead. See also Importing and Exporting Data. $\endgroup$
    – MarcoB
    Jan 7, 2020 at 19:53

This is not an answer to your question but the solution you get is not what you think it is: When you specify

Method -> {"FiniteElement", 
  "MeshOptions" -> {MaxCellMeasure -> 0.005}}

You are solving this a 2D spatial problem but you want to solve this as a 1D spatial and time dependent problem. The easiest way to do so is by using:

heatdist = 
 NDSolve[{HEATIMP == 0, BCTot, ic}, u, {r, 1, 2}, {t, 0, 1}]

And a plot:

enter image description here

If you want to refine the spatial resolution you can do so:

heatdist = 
 NDSolve[{HEATIMP == 0, BCTot, ic}, u, {r, 1, 2}, {t, 0, 1}, 
  Method -> {"PDEDiscretization" -> {"MethodOfLines", 
      "SpatialDiscretization" -> {"FiniteElement", 
        "MeshOptions" -> {MaxCellMeasure -> 0.001}}}}]

Have a look at the introductory FEM tutorial in the documentation.

  • $\begingroup$ Thank you. That clears a few issues up. $\endgroup$
    – I am tired
    Jan 8, 2020 at 16:44

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.