3
$\begingroup$

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.*)
$\endgroup$
  • $\begingroup$ I guess don't try to use non-WorkingPrecision when it's explicitly stated that WorkingPrecision is the only level of precision allowed? $\endgroup$ – Feyre Jun 11 '16 at 12:32
  • $\begingroup$ NDSolve insists on using the IDA package for this particular PDE, and IDA evidently does not accept WorkingPrecision as an option, even WorkingPrecision -> $MachinePrecision. $\endgroup$ – bbgodfrey Jun 11 '16 at 13:09
  • $\begingroup$ 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 $\endgroup$ – bbgodfrey Jun 11 '16 at 13:51
  • $\begingroup$ I'm afraid "PrecisionGoal works fine" isn't related to the issue here: mathematica.stackexchange.com/q/118249/1871 $\endgroup$ – xzczd Jun 7 '18 at 7:47
4
$\begingroup$

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}]

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.

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.