5

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

enter image description here

With WorkingPrecision -> 5

enter image description here

With WorkingPrecision -> 10

enter image description here

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.*)
Mariusz Iwaniuk
  • 13,841
  • 1
  • 25
  • 41
  • 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: https://mathematica.stackexchange.com/q/118249/1871 – xzczd Jun 07 '18 at 07:47

1 Answers1

4

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;
ic = u[-3, y] == 0;
odeic = ic // ptoo;
bc = u[x, -3] == 0;
odebc2 = With[{sf = 10}, ptoo@diffbc[x, sf]@bc];


sollst = NDSolveValue[{Rest@ode, odeic, odebc2}, u /@ grid, {x, -3, 3}, 
   WorkingPrecision -> 16];


sol = rebuild[sollst, grid];
sol[3, 3]
(323.9999999697654)
Plot3D[sol[x, y], {x, -3, 3}, {y, -3, 3}]

Mathematica graphics

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.

xzczd
  • 65,995
  • 9
  • 163
  • 468