differential equations - 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[[1]] /. x -> 3 /. y -> 3

gives error: NDSolve::initf: The initialization of the method NDSolve`StateSpace 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[[1]]/. x -> 3 /. y -> 3
(* 324.*)

Answer

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,
   please find it in the link above. *)
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.