Hi Markus,
Could you provide an inefficient, but accurate, solution pointing out the key required features to capture? When I look at your solutions, it is clear that the low resolution case missed many spikey features. Do you care about those features? Does that imply that the high resolution solution is accurate? There is no point trying to speed up an inaccurate solution. Is there an area that you would like to resolve? Most of the action appears to be centered around low absolute
$\tau$. What are you trying to resolve and has your approach resolved it or do we need to refine further?
I used the following plot functions to view your cases of MinPoints->600 and 2000, respectively.
plt = Plot3D[
Evaluate[Abs[Us[\[Tau], \[Xi]]] /.
Last[sol]], {\[Tau], -\[Tau]max, \[Tau]max}, {\[Xi], 0, \[Xi]max},
PlotPoints -> {100, 100}, MaxRecursion -> 4,
ColorFunction ->
Function[{x, y, z}, Directive[ColorData["DarkBands"][z]]],
ColorFunctionScaling -> False, MeshFunctions -> {#3 &}, Mesh -> 50,
PlotRange -> All, Boxed -> False, Axes -> False, ImageSize -> Large]
sz = 300;
Grid[{{Show[{plt}, ViewProjection -> "Orthographic",
ViewPoint -> Front, ImageSize -> sz,
Background -> RGBColor[0.84`, 0.92`, 1.`], Boxed -> False],
Show[{plt}, ViewProjection -> "Orthographic", ViewPoint -> Right,
ImageSize -> sz, Background -> RGBColor[0.84`, 0.92`, 1.`],
Boxed -> False]}, {Show[{plt}, ViewProjection -> "Orthographic",
ViewPoint -> Top, ImageSize -> sz,
Background -> RGBColor[0.84`, 0.92`, 1.`], Boxed -> False],
Show[{plt}, ViewProjection -> "Perspective",
ViewPoint -> {Above, Right, Front}, ImageSize -> sz,
Background -> RGBColor[0.84`, 0.92`, 1.`], Boxed -> False]}},
Dividers -> Center]
Hi Res Grid
I indicated a "Hair Line" on the figure. I would have assumed it was a simulation artifact, but you might be trying to resolve it given your statement about "fission", but how would one know? I zoomed in on the right and left side.
Hi Res Right Zoom
Hi Res Left Zoom
Low Res Grid
If we zoom in on the
$\tau$ scale, the instability for the low resolution is clearer.
Hi Res Zoom
Low Res Zoom
This may look like fission, but it is probably numerical. If we look at the underside of the Hi Res Zoom, we see some very sharp spikes that the solver needs to detect while solving before it can resolve.
You can view the adaptive mesh actions the solver took to resolve the spikes by selecting the Mesh->All option in Plot3D. If we zoom in on a spike, we can see the refinements.
In the low res case, the spike probably existed within a mesh element.
Applying a Window Function
For
$\tau>5$, the Gaussian function is numerical insignificant relative to your range as alluded to by Alexander. Could we apply a window function using (UnitStep[4-$\tau>5$]-UnitStep[4+$\tau>5$])? When I do that with DifferenceOrder->4, all sorts of splitting happens (4 distinct peaks). Is that real? Also, the solution on my machine is about 3.5 x faster (20 seconds) than without the window function.
Clear["`*"]; ClearAll; ClearAll["Global\[OpenCurlyQuote]*"]; \
Remove["Global`*"];
SetDirectory[NotebookDirectory[]];
\[Delta]3 = 0.00404040404040404;
\[Delta]4 = -4.2087542087542085*^-6;
NS = 4.483404754211932;
\[Xi]max = 1;
\[Tau]max = 30;
Us1[\[Tau]_] =
Exp[-1.3862943611198906*\[Tau]^2] (UnitStep[5 - \[Tau]] +
UnitStep[5 + \[Tau]]);
Plot[Us1[\[Tau]], {\[Tau], -\[Tau]max, \[Tau]max}, Frame -> True,
PlotStyle -> Darker[Blue], PlotRange -> All, Filling -> Bottom];
eqs = {D[Us[\[Tau], \[Xi]] , \[Xi]] -
I/2*D[Us[\[Tau], \[Xi]], {\[Tau], 2}] - \[Delta]3*
D[Us[\[Tau], \[Xi]], {\[Tau], 3}] -
I*\[Delta]4*D[Us[\[Tau], \[Xi]], {\[Tau], 4}] -
I*NS^2*Abs[Us[\[Tau], \[Xi]]]^2*Us[\[Tau], \[Xi]] == 0};
bc = {Us[\[Tau]max, \[Xi]] == 0,
Us[-\[Tau]max, \[Xi]] ==
0, (D[Us[\[Tau], \[Xi]], {\[Tau], 1}] /. \[Tau] -> \[Tau]max) ==
0, (D[Us[\[Tau], \[Xi]], {\[Tau], 1}] /. \[Tau] -> -\[Tau]max) ==
0};
ic = {Us[\[Tau], 0] == Us1[\[Tau]]};
Print["seconds to obtain solution: ",
time = AbsoluteTiming[
Monitor[
sol =
NDSolve[Join[eqs, ic, bc],
Us, {\[Tau], -\[Tau]max, \[Tau]max}, {\[Xi], 0, \[Xi]max},
EvaluationMonitor :> (\[Xi]current = \[Xi]),
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
"MinPoints" -> 2000, "DifferenceOrder" -> 4}},
PrecisionGoal -> 5, AccuracyGoal -> 5,
MaxSteps -> \[Infinity]],
"current \[Xi] (\!\(\*SubscriptBox[\(\[Xi]\), \(max\)]\)=" <>
ToString[\[Xi]max] <> "): " <> ToString[ \[Xi]current]];][[
1]]] ;
Print["grid points along {\[Tau],\[Xi]}-directions: ",
Length /@ InterpolatingFunctionCoordinates[Us /. sol[[1]]]];
(UnitStep[4-$\tau$]-UnitStep[4+$\tau$])
(UnitStep[5-$\tau$]-UnitStep[5+$\tau$])
There seems to be minor change going from 4 to 5 on the window function, indicating little sensitivity in that range.
Zoomed
(UnitStep[6-$\tau$]-UnitStep[6+$\tau$])
Here is an animation along the
$\xi$ coordinate for an alternative visualization.
Summary
We may or may not be making progress. Adding a window function when the Gaussian becomes insignificant, seems to induce splitting. We will need some more definition from you about your requirements and what the expected results are to make any additional progress. If we are on the right path, then a 20 seconds solve time should enable 4000 simulations per day.
