How to make a single-variable out of a NDSolve solution?

Question

So I've got this line that contains the solution of a partial, non-extraordinary differential equation (because Mathematica doesn't handle extraordinary partial differential equations):

phi6m = NDSolveValue[{D[u[t, x], t, t] - 
 D[u[t, x], x, x] == -6 u[t, x]^5 + 10.5 u[t, x]^3 - 4.5 u[t, x],    u[0, x] == Tanh[x], Derivative[1, 0][u][0, x] == 0, u[t, -7] == -1,
u[t, 7] == 1}, u, {t, 0, 6}, {x, -7, 7}]

And now I wish to compute the following as a function of time, as well as its derivative with respect to time:

$ E(t)=\int_{-7}^{7} dx (\frac{\partial u}{\partial x})^2+(u(t,x)^2-1)^2(u(t,x)^2-0.625)$

where u is the solution given by the aforementionned line of code. But now I wish to take the function, as well as its derivative, and plot it into a graph... what I tried allowed me to take the function at a fixed time value, here t=0.1:

phi62 = phi6m[0.1, x];
energiephi6 = D[phi62, x]^2 + (phi62^2 - 1)^2 (phi62^2 - 0.625);
energietotale = NIntegrate[energiephi6, {x, -7, 7}]

But not to make something plottable as a function of time, much less a derivative of the aforementionned function with respect to time. Is there anything else I can do to resolve the issue?

Answer 1

I don't think this is what you want, but nevertheless I think it is more of an answer rather than a comment, so will be bold and post it.

First, you can't obtain an analytic function as a solution to a definite integral. Second, there are some problems with the PDE itself. To overcome them we can do something like that.

phi6m = NDSolveValue[{D[u[t, x], t, t] - D[u[t, x], x, x] == -6 u[t, x]^5 +
               10.5 u[t, x]^3 - 4.5 u[t, x], u[0, x] == Tanh[x],
               Derivative[1, 0][u][0, x] == 0, u[t, -7] == -1,u[t, 7] == 1},
               u, {t, 0, 6}, {x, -7, 7}, 
               Method -> {"PDEDiscretization" -> {"MethodOfLines",
                  "SpatialDiscretization" -> {"TensorProductGrid", "MinPoints" -> 10000}}}]

Which produces a nice looking plot

Plot3D[phi6m[x, y], {x, 0., 6.}, {y, -7., 7.}, ImageSize -> 500]

contour = ContourPlot[{phi6m[x, y], phi6m'[x, y]}, {x, 0.`, 6.`}, {y, -7.`, 7.`},
                       PlotPoints -> 500];

dens = DensityPlot[{phi6m[x, y], phi6m'[x, y]}, {x, 0.`, 6.`}, {y, -7.`, 7.`},  
                   ColorFunction -> "CandyColors", PlotPoints -> 500]

And you can, of course, overlay them

Show[{dens, contour}]

Answer 2

You are almost there. Just make a function from the code you wrote for a single time value:

energyTotale[t_?NumericQ] := Module[{phi62, energiephi6},
  phi62 = phi6m[t, x];
  energiephi6 = D[phi62, x]^2 + (phi62^2 - 1)^2 (phi62^2 - 0.625);
  NIntegrate[energiephi6, {x, -7, 7}]
  ]

You can now in principle use that as a function to be plotted. You can not as easily get a derivative of that. I would thus suggest to create a new interpolating function by evaluating the above function for several values in the range of interest. Of course that means you have to make a decision about how many points you need for a result good enough for your purposes. This visually looked pretty good for me, but you might want to invest some extra computing time in a finer grid:

energy6mData = Table[{t, energyTotale[t]}, {t, 0, 1, 0.025}]

ListLinePlot[energy6mData, Frame -> True, Axes -> False]

From such data, it is easy enough to create an InterpolatingFunction. Here I used the standard polynomial interpolation with InterpolationOrder 3 as that will creat a "nice" continuous and differentiable derivative, but of course choosing a decent interpolation order and method might depend on the details of your problem, you might want to read the Interpolation documentation and "tutorial/ApproximateFunctionsAndInterpolation" for details:

energy6m = Interpolation[energy6mData, InterpolationOrder -> 3]

You can check the "quality" of that interpolation:

Show[
 ListLinePlot[energy6mData, Frame -> True, Axes -> False],
 Plot[energy6m[t], {t, 0.1, 1}, PlotStyle -> Red]
 ]

You can also try FunctionInterpolation as an automation of the above which I think tries to do something smart (probably similar to what Plot does) and not necessarily uses equidistant grids. It might evaluate quite often and thus be very expensive in some cases, though. You probably need to try out what works best for your case:

energy6m = FunctionInterpolation[energyTotale[t],{t, 0, 1}]

From such an InterpolatingFunction, it is straightforward to get a derivative (I prefer to treat those as "function objects" which is why I only use the head without the argument):

energy6mD = Head@D[energy6m[t], t]

EDIT: I just remembered that there is of course a more elegant way to get that derivative, namely:

energy6mD = Derivative[1][energy6m]

This is now again a function which you can evaluate, plot, integrate...:

Plot[energy6mD[t], {t, 0, 1}]