2
$\begingroup$

The following code works but it give several warnings that makes me feel nervous about the reliability and accuracy of the output:

NDSolve: Warning: estimated initial error on the specified spatial grid in the direction of independent variable x exceeds prescribed error tolerance.

How can I fix the warnings? Take note that I need an accurate and high resolution output.

The code below solves a partial differential equation for the propagation of a scalar field in 2D:

Clear["Global`*"]

size = 30; (* max distance from origin *)
simulationtime = 20; (* time of the simulation *)
alpha = 2; (* non-linear interaction parameter. *)

fieldequation = (D[field[t, x, y], t, t] - D[field[t, x, y], x, x] - D[field[t, x, y], y, y] + 2 alpha^2 (field[t, x, y]^2 - 1) field[t, x, y] == 0);

skeleton = Table[{x, y, RandomReal[{-1, 1}]}, {x, -size, size, 2}, {y, -size, size, 2}];

(* smooth initial conditions for the scalar field : *)

fluctuations = Interpolation[Flatten[skeleton, 1], Method -> "Spline"];
phi0[t_, x_, y_] = fluctuations[x, y]; (* static initial scalar field *)

bordersconditions = {
    field[t, -size, y] == phi0[t, -size, y],
    field[t, size, y] == phi0[t, size, y],
    field[t, x, -size] == phi0[t, x, -size],
    field[t, x, size] == phi0[t, x, size]
};

initialconditions = {
    field[0, x, y] == phi0[0, x, y],
    (D[field[t, x, y], t] /. t -> 0) == D[phi0[t, x, y], t] /. t -> 0
};

fieldsolution = NDSolve[
    Flatten@{fieldequation, initialconditions, bordersconditions},
    field,
    {t, 0, simulationtime},
    {x, -size, size},
    {y, -size, size},
    Method -> {
    "MethodOfLines",
    "SpatialDiscretization" -> {
        "TensorProductGrid",
        "MaxPoints" -> 200, (* slow compilation if > 200! *)
        "MinPoints" -> 200, (* less reliable if < 100. *)
        "DifferenceOrder" -> 4  (* At least 2.  Very slow if > 6. *)
    }}];

Manipulate[Plot3D[Evaluate[field[t, x, y] /. fieldsolution /. t -> time],
    {x, -10, 10}, {y, -10, 10},
    PlotPoints -> ControlActive[20, 60],
    PlotRange -> {{-10, 10}, {-10, 10}, {-3, 3}}],
    {{time, 0, "t"}, 0, simulationtime, 0.1}
]

Preview of what this code is doing:

enter image description here

$\endgroup$
9
  • $\begingroup$ What happens if you have a smoother initial field (not random)? Does NDSolve still have trouble? $\endgroup$
    – Michael E2
    Jul 30, 2023 at 17:24
  • $\begingroup$ @MichaelE2, I tried three very different random initial conditions, and I get the same warnings. $\endgroup$
    – Cham
    Jul 30, 2023 at 17:33
  • $\begingroup$ I was thinking of something like {x, y, Sin[x/30] Sin[y/30]}. I suspect it's the interpolation of random initial conditions that leads to the NDSolve::??? message. I was suggesting you test take. But if the minimal-fluctuation fields gives the same message, then it's probably something else. $\endgroup$
    – Michael E2
    Jul 30, 2023 at 17:37
  • $\begingroup$ @MichaelE2, I just tried a smooth (non random) initial condition, and still get similar warnings. $\endgroup$
    – Cham
    Jul 30, 2023 at 17:38
  • $\begingroup$ @MichaelE2, I also tried by turning off the non-linear part of the field equation (using alpha = 0). I get less warnings, but some are still there. I suspect I need to fine-tune the parameters "MaxPoints" -> 200, "MinPoints" -> 200, "DifferenceOrder" -> 4 $\endgroup$
    – Cham
    Jul 30, 2023 at 17:42

2 Answers 2

1
$\begingroup$

No error message if you use "FiniteElement"

reg = Rectangle[{-size, -size}, {size, size}]; 
fieldsolution = 
 NDSolveValue[
  Flatten@{fieldequation, initialconditions, bordersconditions}, 
  field, {t, 0, simulationtime}, Element[{x, y }, reg], 
  Method -> {"MethodOfLines", 
    "SpatialDiscretization" -> {"FiniteElement" , 
      "MeshOptions" -> {"MaxCellMeasure" -> 1, 
        "MeshElementType" -> "TriangleElement", "MeshOrder" -> 1}  }}]

Manipulate[
 Plot3D[Evaluate[ fieldsolution[time, x, y]], {x, -10, 10}, {y, -10, 
   10}, PlotPoints ->100,Mesh->None, 
  PlotRange -> {{-10, 10}, {-10, 10}, {-3, 3}}], {{time, 0, "t"}, 0, 
  simulationtime, 0.1}] 

enter image description here

addendum

As recommended by @MichaelE2 here the solution for smooth initial conditions Cos[Pi/2 x/30] Cos[Pi/2 y/30] FEM shows traveling waves as expected

Clear[phi0]
phi0[t_, x_, y_] := Cos[Pi/2 x/30] Cos[Pi/2 y/30];

(*static initial scalar field*)(*bordersconditions={field[t,-size,y]\
\[Equal]phi0[t,-size,y],field[t,size,y]\[Equal]phi0[t,size,y],field[t,\
x,-size]\[Equal]phi0[t,x,-size],field[t,x,size]\[Equal]phi0[t,x,size]}\
;*)
bordersconditions = 
 DirichletCondition[field[t, x, y] == phi0[t, x, y], True]
initialconditions = {field[0, x, y] == phi0[0, x, y], 
   Derivative[1, 0, 0][field][0, x, y] == 
    Derivative[1, 0, 0][phi0][0, x, y] };
reg = Rectangle[{-size, -size}, {size, size}]; 
fieldsolution = 
 NDSolveValue[
  Flatten@{fieldequation, initialconditions, bordersconditions}, 
  field, {t, 0, simulationtime}, Element[{x, y }, reg], 
  Method -> {"MethodOfLines", "TemporalVariable" -> t, 
    "SpatialDiscretization" -> {"FiniteElement" }}
  ]

enter image description here

The mesh in this simulation isn't fine, FEM only needs 400 quad-elements

$\endgroup$
14
  • $\begingroup$ It works, but the output is very ugly with all these triangular shapes. How to rise up the resolution considerably? $\endgroup$
    – Cham
    Jul 30, 2023 at 16:38
  • $\begingroup$ I see that using "MaxCellMeasure" -> 0.1 raises the resolution, but it's still crude and takes more time to compile than the original version. $\endgroup$
    – Cham
    Jul 30, 2023 at 16:44
  • $\begingroup$ You asked for a code without error message! $\endgroup$ Jul 30, 2023 at 16:46
  • 1
    $\begingroup$ Well, yes, but not by changing the output that much! I need an accurate and high resolution output. $\endgroup$
    – Cham
    Jul 30, 2023 at 16:55
  • 1
    $\begingroup$ @Cham "MeshOrder" -> 2 should give a more accurate interpolation of the solution (as well as somewhat more accurate mesh elements). To get a better solution, use "MaxCellMeasure" (as you've already tried and rejected on the grounds that it took longer than your unsatisfactory solution). Your initial field oscillates between values every $\Delta x = \Delta y = 2$. So you need a mesh/grid that's a quite a bit smaller than this, I would think. Eight points per cycle (picked out of a hat) gives "MinPoints" -> 240, around what your setting is in the OP. Maybe it needs to be more than 8. $\endgroup$
    – Michael E2
    Jul 30, 2023 at 17:48
0
$\begingroup$

You can just use Quiet@ or Quiet[] and the code will work without the warning

$\endgroup$
3
  • 4
    $\begingroup$ I guess that using this would prevent to see errors or important warnings about accuracy. $\endgroup$
    – Cham
    Jul 30, 2023 at 17:27
  • $\begingroup$ Normally, yes - I'd agree. But the code seems to be working in this case $\endgroup$
    – Navvye
    Jul 30, 2023 at 17:29
  • 3
    $\begingroup$ @Cham I was going to suggest this in a COMMENT, but with a ":)", since it's not a serious answer. $\endgroup$
    – Michael E2
    Jul 30, 2023 at 17:29

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.