# Extracting data from InterpolatingFunction 2D surface plot

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).

Needs["NDSolveFEM"]
Needs["DifferentialEquationsInterpolatingFunctionAnatomy"];
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}]

• You are not solving a 1D diffusion equation but a 2D convection diffusion equation. See my answer. Jan 8 '20 at 8:15

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}}
];

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

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}]]

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

• 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? Jan 7 '20 at 18:46
• @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. Jan 7 '20 at 19:23
• @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? Jan 7 '20 at 19:25
• 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. Jan 7 '20 at 19:47
• @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. Jan 7 '20 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:

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.

• Thank you. That clears a few issues up. Jan 8 '20 at 16:44