I want to get the raw data of solutions of partial differential equation in mathematica

Question

I have one partial differential equation.

sol = NDSolve[{ 
    D[T[x, y, t], {t, 1}] == 
     (D[T[x, y, t], {x, 2}] + D[T[x, y, t], {y, 2}]), 
    T[x, y, 0] == 
     400 - 400 Exp[-100000 x^2 y^2 (x - 1)^2 (y - 1)^2] + 
      350 Exp[-100000 x^2  (x - 1)^2],
    T[x, 0, t] == 350 Exp[-100000 x^2  (x - 1)^2], 
    T[x, 1, t] == 350 Exp[-100000 x^2  (x - 1)^2], T[0, y, t] == 350, 
    T[1, y, t] == 350}, {T[x, y, t]}, {x, 0, 1}, {y, 0, 1}, {t, 0, 
    100}, PrecisionGoal -> 2];

This sol is an InterpolationgFunctiion. I want to get the raw numerical data from this InterpolatingFunction.

Do I use InterpolatingFunctionValuesOnGrid or InterpolatingFunctionCoordinates, I don't understand these tools.

How do I extract data from above function?

Welcome! I suggest the following:

As you receive help, try to give it too, by answering questions in your area of expertise.

Take the tour and check the faqs!

When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! — , Mar 26 '16 at 08:37

score 3 · Answer 1 · answered Mar 26 '16 at 06:26

3

Try this

Needs["DifferentialEquations`InterpolatingFunctionAnatomy`"];
sol = T/.NDSolve[{D[T[x,y,t], {t,1}]==(D[T[x,y,t], {x,2}] + D[T[x,y,t], {y,2}]), 
  T[x,y,0]==400-400 Exp[-100000 x^2 y^2(x-1)^2 (y-1)^2]+ 350 Exp[-100000 x^2 (x-1)^2], 
  T[x,0,t]==350 Exp[-100000 x^2 (x-1)^2], 
  T[x,1,t]==350 Exp[-100000 x^2 (x-1)^2],
  T[0,y,t]==350,
  T[1,y,t]==350},
  T, {x,0,1}, {y,0,1}, {t,0,100}, PrecisionGoal->2][[1]];
InterpolatingFunctionValuesOnGrid[sol]

The output is large enough that it shows a shortened version. You can then click on show more (repeatedly) or click on show all and wait..... for it.

That should show you a large grid of result values.

You can assign that grid of values to a variable and use subscripting to extract individual items.

answered Mar 26 '16 at 06:26

Bill

12,001
12
13

Alternatively, sol["ValuesOnGrid"] – J. M.'s missing motivation Mar 26 '16 at 06:28
What does [[1]] mean? – Shakariky Mar 26 '16 at 07:10
@Shakariky Click here -> https://reference.wolfram.com/language/ref/Part.html – Bill Mar 26 '16 at 07:23
What does T/. mean? Sorry. – Shakariky Mar 26 '16 at 07:42
@Shak, whenever you encounter an unfamiliar symbol in Mathematica, just highlight it and press F1. – J. M.'s missing motivation Mar 26 '16 at 07:55

score 0 · Answer 2 · edited Apr 13 '17 at 12:56

Bill and J.M. has already answered OP's question, but after some communication with OP in the comment under this post, I think it's worth explaining the usage of InterpolatingFunctionValuesOnGrid[sol] / sol["ValuesOnGrid"] and InterpolatingFunctionCoordinates[sol] / sol["Coordinates"] a bit further.

As the function names suggest, sol["ValuesOnGrid"] is a list of the function value on the grid. (For OP's case, it's a list of value of $T$.) sol["Coordinates"] is a list of the coordinate of grid. (For OP's case, it's a list of value of $x$,$y$,$t$.) If you want to build a list of $(x,y,t,T)$ that can be used in e.g. ListDensityPlot3D, you can:

sol = T /. First@
    NDSolve[{D[T[x, y, t], {t, 1}] == (D[T[x, y, t], {x, 2}] + D[T[x, y, t], {y, 2}]), 
      T[x, y, 0] == 
       400 - 400 Exp[-100000 x^2 y^2 (x - 1)^2 (y - 1)^2] + 
        350 Exp[-100000 x^2 (x - 1)^2], T[x, 0, t] == 350 Exp[-100000 x^2 (x - 1)^2], 
      T[x, 1, t] == 350 Exp[-100000 x^2 (x - 1)^2], T[0, y, t] == 350, 
      T[1, y, t] == 350}, T, {x, 0, 1}, {y, 0, 1}, {t, 0, 100}, PrecisionGoal -> 2];

lstT = sol["ValuesOnGrid"];
lstxyt = {lstx, lsty, lstt} = sol["Coordinates"];
{lenx, leny, lent} = Dimensions@lstT;

(* solution 1, not that fast but easy to understand *)
lstxytT = Table[{lstx[[i]], lsty[[j]], lstt[[k]], lstT[[i, j, k]]}, 
                {i, lenx}, {j, leny}, {k, lent}];

(* solution 2, a bit harder to understand but fast *)
lstxytT2 = Transpose[{Transpose[ConstantArray[lstx, {leny, lent}], {2, 3, 1}], 
                      Transpose[ConstantArray[lsty, {lenx, lent}], {1, 3, 2}], 
                                ConstantArray[lstt, {lenx, leny}], 
                                              lstT}, 
                     {4, 1, 2, 3}];

lstxytT == lstxytT2
(* True *)

ListDensityPlot3D@lstxytT
(* I'm still in v9 so let me omit the picture here *)

Finally I'd like to mention that lstT and lstxyt can be directly used for the building of a interpolating function of $T(x,y,t)$:

solsimpler = ListInterpolation[lstT, lstxyt]

Notice solsimpler isn't the same as sol because it's only built with "ValuesOnGrid" and "Coordinates". To learn more about the information contained in InterpolatingFunction, check this and this post.

I want to get the raw data of solutions of partial differential equation in mathematica

2 Answers2