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?

  • $\begingroup$ You should add more details. For example: Order of the derivative ? Are you trying to produce phase plots ? $\endgroup$ – Sektor Feb 7 '14 at 17:14
  • $\begingroup$ The derivative is first order in time, and I understand that the function is not analytical, and neither will the derivative be. Nevertheless I want to plot both the function given above as a function of time and the numerical value of its derivative (albeit separately)... $\endgroup$ – NSERC Protester Feb 7 '14 at 17:44

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]

Mathematica graphics

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]

Mathematica graphics

And you can, of course, overlay them

Show[{dens, contour}]

Mathematica graphics


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:

 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}] 
  • $\begingroup$ There's a typo in your code; just below ".. check the "quality" of that ..." - energy[t] instead of energy6m[t] $\endgroup$ – Sektor Feb 9 '14 at 22:24
  • $\begingroup$ @Sektor: thanks for pointing out, I actually changed that line so it shows the data that I was explicitly generating before already, which shows more clearly what I try to point out. Are you aware of the fact that you can in such cases edit the answer yourself? I think such corrections are common and welcomed on this site... $\endgroup$ – Albert Retey Feb 10 '14 at 9:08
  • $\begingroup$ Yes, of course, I am aware of such a functionality, but I thought I should just point you, so you can preserve the consistency of the post :) $\endgroup$ – Sektor Feb 10 '14 at 16:42
  • $\begingroup$ @Sektor: I of course appreciate your politeness :-) $\endgroup$ – Albert Retey Feb 11 '14 at 21:05
  • $\begingroup$ You are welcome :) $\endgroup$ – Sektor Feb 11 '14 at 21:08

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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