# WorkingPrecision in NDSolve causes failure when solving a simple PDE

I have fivefold multiple integral and I wanted a speed calculations. I came across on this Question and had already begun the problems.

Here is a toy example with simple double integral:

f[x_, y_] := x^2*y^2;
NIntegrate[f[x, y], {y, -3, 3}, {x, -3, 3}, WorkingPrecision -> 20]
(* 324.00000000000000000 *)

sol2 = NDSolve[{D[u[x, y], x, y] == f[x, y], u[-3, y] == 0,
u[x, -3] == 0}, u, {x, -3, 3}, {y, -3, 3}];
u[x, y] /. sol2 /. x -> 3 /. y -> 3

(* 324. *)


Yes works,but I change WorkingPrecision in NDSolve they begin to happen strange things:

With WorkingPrecision -> 2 With WorkingPrecision -> 5 With WorkingPrecision -> 10 I have MMA on Windows 8.1 "10.2.0 for Microsoft Windows (64-bit) (July 7, 2015)"

EDITED: 07.06.2018

On Mathematica 11.3:

f[x_, y_] := x^2*y^2;
sol2 = NDSolve[{D[u[x, y], x, y] == f[x, y], u[-3, y] == 0,
u[x, -3] == 0}, u, {x, -3, 3}, {y, -3, 3}, WorkingPrecision -> 5];
u[x, y] /. sol2 /. x -> 3 /. y -> 3


or:

f[x_, y_] := x^2*y^2;
sol2 = NDSolve[{D[u[x, y], x, y] == f[x, y], u[-3, y] == 0,
u[x, -3] == 0}, u, {x, -3, 3}, {y, -3, 3},WorkingPrecision -> 20,
PrecisionGoal -> 8, AccuracyGoal -> 8];
u[x, y] /. sol2[] /. x -> 3 /. y -> 3


gives error: NDSolve::initf: The initialization of the method NDSolveStateSpace failed..

Issue still persist on MMA 11.3.

With no options like:WorkingPrecision,PrecisionGoal, AccuracyGoal works fine.

f[x_, y_] := x^2*y^2;
sol2 = NDSolve[{D[u[x, y], x, y] == f[x, y], u[-3, y] == 0,
u[x, -3] == 0}, u, {x, -3, 3}, {y, -3, 3}];
u[x, y] /. sol2[]/. x -> 3 /. y -> 3
(* 324.*)

• I guess don't try to use non-WorkingPrecision when it's explicitly stated that WorkingPrecision is the only level of precision allowed? – Feyre Jun 11 '16 at 12:32
• NDSolve insists on using the IDA package for this particular PDE, and IDA evidently does not accept WorkingPrecision as an option, even WorkingPrecision -> \$MachinePrecision. – bbgodfrey Jun 11 '16 at 13:09
• Incidentally, DSolve cannot solve the particular PDE either, even though the solution is obvious: 81 + 3 x^3 + 3 y^3 + (x^3 y^3)/9 – bbgodfrey Jun 11 '16 at 13:51
• I'm afraid "PrecisionGoal works fine" isn't related to the issue here: mathematica.stackexchange.com/q/118249/1871 – xzczd Jun 7 '18 at 7:47

Very strange, NDSolve seems to insist on discretizing the PDE to a DAE system, but it should be able to discretize it to a set of ODEs rather than DAEs. The following is a self implementation of method of lines. I'll use pdetoode for discretization.

f[x_, y_] := x^2*y^2;
grid = Array[# &, 25, {-3, 3}];
(* Definition of pdetoode isn't included in this code piece,
ptoo = pdetoode[u[x, y], x, grid, 4];
ode = D[u[x, y], x, y] == f[x, y] // ptoo;
odeic = u[-3, y] == 0 // ptoo;
odebc2 = With[{sf = 10}, Map[# + sf D[#, x] &, ptoo[u[x, -3] == 0]]];

sol = rebuild[sollst, grid];
sol[3, 3]
(*323.9999999697654*)
Plot3D[sol[x, y], {x, -3, 3}, {y, -3, 3}] As far as I can tell, the above code is just a self implementation of WorkingPrecision -> 16, Method -> {"MethodOfLines", "DifferentiateBoundaryConditions" -> {True, "ScaleFactor" -> 10}, "SpatialDiscretization" -> {"TensorProductGrid", "MaxPoints" -> 25, "MinPoints" -> 25, "DifferenceOrder" -> 4}}, but adding these options to NDSolve won't resolve the problem. Maybe it's a bug, maybe it's due to another undocument feature of "MethodOflines"`. I suggest you to report it to WRI.

• I'm raported to WRI. – Mariusz Iwaniuk Jan 3 '17 at 12:35