# Unable to improve accuracy of MethodOfLines

I'm studying a PDE which roughly models a rope subjected to driving on one end. The coordinate $$v$$ is related to physical time $$t$$ via the relation $$v=t-x$$. I'd like to solve the initial value problem using the method of lines. The following Mathematica code returns a solution, but the accuracy is rather low. Passing the WorkingPrecision option to NDSolveValue throws the following error:

The initialization of the method NDSolveStateSpace failed

How can one systematically improve the solution accuracy?

Clear["Global*"]

vmax = 10;
f[v_] := Sin[v]^2;

method := {
"MethodOfLines",
"TemporalVariable" -> v,
"SpatialDiscretization" -> {
"TensorProductGrid",
"DifferenceOrder" -> "Pseudospectral", MaxPoints -> 100}};

{timing, uinterp} = Timing@NDSolveValue[{
2 Derivative[1, 1][u][x, v] == Derivative[2, 0][u][x, v],
u[x, 0] == 0,
u[0, v] == f[v],
u[1, v] == 0}
, u, {x, 0, 1}, {v, 0, vmax}, Method -> method, MaxStepSize -> 0.1]

Animate[Plot[If[x < t, uinterp[x, t - x], 0], {x, 0, 1}, PlotRange -> {-1, 1.1}, AxesLabel -> {x, u}], {t, 0, vmax}]


• what do you write [u][x, v] when v is supposed to be time? Why not use normal notation and use t there, just like you did when you defined f[t_]? Also in the plot you changed to t. ofcourse the computer does not care what letter you use, but it makes your code confusing to read. Commented Dec 2, 2023 at 16:26
• When I removed all these options and methods you had added, I get smooth solution. Did you try that? I always start with nothing added. And only add options and change settings if something is not right. Try this NDSolveValue[{2 Derivative[1,1][u][x,v]==Derivative[2,0][u][x,v],u[x,0]==0,u[0,v]==f[v],u[1,v]==0},u,{x,0,1},{v,0,T}] Commented Dec 2, 2023 at 16:30
• @Nasser I've addressed your question about v vs t in the revised question. Thanks - it's interesting that removing the method significantly improved the quality of the solution, but my question is specifically about improving the accuracy of MethodOfLines. The PDE in the OP is a minimal model for a more complicated PDE that fails without method specification. Commented Dec 2, 2023 at 17:09
• You can check these question which seem to be related mathematica.stackexchange.com/questions/118194/… and see end of this post community.wolfram.com/groups/-/m/t/2933255 might be behind the problem, and that gave NDSolve The initialization of the method NDSolve StateSpace failed. That might be the hint that someone needs to let you know how to change methods. They show same error message you got. Commented Dec 2, 2023 at 17:30

Once again, rule of thumb: if NDSolve doesn't work well on solving PDE, and you're sure you haven't made any simple mistake, then usually the problem lies in spatial discretization.

So the WorkingPrecision should be the last thing to touch. Let's adjust the spatial and temporal grid a bit:

mol[n : _Integer | {_Integer ..}, o_ : "Pseudospectral"] := {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid", "MaxPoints" -> n,
"MinPoints" -> n, "DifferenceOrder" -> o}}

vmax = 10;
f[v_] := Sin[v]^2;

{timing, uinterp} =
Timing@NDSolveValue[{2  Derivative[1, 1][u][x, v] == Derivative[2, 0][u][x, v],
u[x, 0] == 0, u[0, v] == f[v], u[1, v] == 0}, u, {x, 0, 1}, {v, 0, vmax},
Method -> mol[200, 2], MaxStepSize -> {Automatic, 0.001}]
(* {0.65625, …} *)

Plot3D[If[x < t, uinterp[x, t - x], 0], {x, 0, 1}, {t, 0, vmax},
PlotRange -> {-1, 1.1}, PlotPoints -> 50, Exclusions -> None]


Animate[Plot[If[x < t, uinterp[x, t - x], 0], {x, 0, 1}, PlotRange -> {-1, 1.1},
AxesLabel -> {x, u}], {t, 0, vmax}]


In example above I've used second order spatial discretization. "Pseudospectral" method can be used if you like, just set e.g. mol[101]. Interestingly, unlike the problem in your previous question, fourth order grid e.g. mol[100, 4] also produces reasonable result.

You can adjust the parameters related to discretization further to see what will happen.

• Thanks for your answer and useful advice. However, I'm still left wondering why WorkingPrecision causes MethodOfLines to fail. Perhaps implementing my own method of lines routine is the only workaround in Mathematica. Commented Dec 3, 2023 at 16:38
• I was wondering if you could also explain the meaning of this line: MaxStepSize -> {Automatic, 0.001}? I found it to be crucial for getting good results with the pseudospectral method Commented Dec 3, 2023 at 17:08
• @phonon It's max step sizes for the 2 directions of the domain. Automatic is the max step size for x direction (the step size in x direction is already controled by "MaxPoints" and "MinPoints" in mol, so we set it to Automatic to avoid possible conflict, you can try setting it to an arbitrary number to see what has happened), and 0.001 is for v direction. I've also explained this a bit in this post: mathematica.stackexchange.com/a/229184/1871 Commented Dec 4, 2023 at 0:12