2
$\begingroup$

I try to find a general approach to plot the time evolution of horizontal distances from maximum to maximum in a solution of PDE. The solution u[x,t] normally have multiple maximum and minimum in space x, which move in space x and evolve in time t.

Here is a simple example, in which the maxima and minima are periodic. But in my real problem they are not periodic and the distances between different pairs of adjacent max are different at a given t, also the distances between two adjacent max can change with t.

sol = NDSolve[{D[u[x, t], t] + u[x, t] D[u[x, t], x] + D[u[x, t], x, x] + 
0.4*D[u[x, t], {x, 3}] + D[u[x, t], {x, 4}] == 0,
u[-4 \[Pi], t] == u[4 \[Pi], t], u[x, 0] == 0.1*Sin[x]}, u, {t, 0, 20},
{x, -4 \[Pi], 4 \[Pi]}]

Plot3D[Evaluate[u[x, t] /. First[sol]], {t, 0, 10}, {x, -4 Pi, 4 Pi}, PlotRange -> All, PlotPoints -> 100]

enter image description here

I have tried to use Table[FindMaximum[Evaluate[u[x, t] /. First[sol]], {x, x0}][[2, 1, 2]], {t,0,tend,0.01}] with an initial position x0 to find a local maximum. But I don't know how to find two adjacent maxima simultaneously in order to plot the time evolution of their distance.

$\endgroup$

1 Answer 1

1
$\begingroup$

A possible approach is :

  • calculate the norm of the gradient of the interpolating function (to see where it is 0).
  • transform the resulting function in a MeshRegion
  • Apply RegionIntersection to the MesgRegion and a plane near z=0.

The result is a collection of lines, which you can further process.

Your code slightly modified for convenience :

sol1 = NDSolveValue[{D[u[x, t], t] + u[x, t] D[u[x, t], x] + 
      D[u[x, t], x, x] + 0.4*D[u[x, t], {x, 3}] + 
      D[u[x, t], {x, 4}] == 0, u[-4 \[Pi], t] == u[4 \[Pi], t], 
    u[x, 0] == 0.1*Sin[ x]}, 
   u[x, t], {t, 0, 20}, {x, -4 \[Pi], 4 \[Pi]}];
Plot3D[sol1, {t, 0, 10}, {x, -4 Pi, 4 Pi}, PlotRange -> All, 
 PlotPoints -> 100]  

enter image description here

The norm of the gradient, and a plot :

gradNorm = Norm[Grad[sol1, {x, t}]]
Plot3D[{sol1, gradNorm}, {t, 0, 10}, {x, -4 Pi, 4 Pi}, 
 PlotRange -> All, PlotPoints -> 100]  

enter image description here

The region and its plot with Mesh-> All. The option Mesh->All permits to understand what we will obtain later.

region = DiscretizeGraphics @ 
   Plot3D[gradNorm, {t, 0, 10}, {x, -4 Pi, 4 Pi}, PlotRange -> All, 
    PlotPoints -> 100];
RegionPlot3D[region, Mesh -> All, BoxRatios -> {1, 1, .2}, 
 Axes -> True]  

enter image description here

The plane at z=0.01. The value 0.01 comes from what we see in the plot just above (To be tuned), and finally the intersection.

zcut = 0.01;
myPlane = InfinitePlane[{{0, 0, zcut}, {1, 0, zcut}, {0, 1, zcut}}]

intersection00 = 
  RegionIntersection[
      myPlane, #] & /@ (MeshPrimitives[region, 2]) // 
   DeleteCases[#, 
     EmptyRegion[
      3]] &; (* takes 80 seconds on my 9 years old intel I7 2630QM *)
\

intersection00[[1]] (* just to see what it is, out of curiosity *)
Graphics3D[intersection00]  

Line[{{9.82633, -11.7335, 0.01}, {9.89899, -11.7277, 0.01}}]

enter image description here

Notes :

  • this gives the maxima and the minima

  • the idea tu use RegionIntersection[.. [ MeshPrimitives comes frome here

$\endgroup$

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.