Skip to main content

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



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.*)

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

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.


Comments

Post a Comment

Popular posts from this blog

front end - keyboard shortcut to invoke Insert new matrix

I frequently need to type in some matrices, and the menu command Insert > Table/Matrix > New... allows matrices with lines drawn between columns and rows, which is very helpful. I would like to make a keyboard shortcut for it, but cannot find the relevant frontend token command (4209405) for it. Since the FullForm[] and InputForm[] of matrices with lines drawn between rows and columns is the same as those without lines, it's hard to do this via 3rd party system-wide text expanders (e.g. autohotkey or atext on mac). How does one assign a keyboard shortcut for the menu item Insert > Table/Matrix > New... , preferably using only mathematica? Thanks! Answer In the MenuSetup.tr (for linux located in the $InstallationDirectory/SystemFiles/FrontEnd/TextResources/X/ directory), I changed the line MenuItem["&New...", "CreateGridBoxDialog"] to read MenuItem["&New...", "CreateGridBoxDialog", MenuKey["m", Modifiers-...
Post a Comment
Read more

How to thread a list

I have data in format data = {{a1, a2}, {b1, b2}, {c1, c2}, {d1, d2}} Tableform: I want to thread it to : tdata = {{{a1, b1}, {a2, b2}}, {{a1, c1}, {a2, c2}}, {{a1, d1}, {a2, d2}}} Tableform: And I would like to do better then pseudofunction[n_] := Transpose[{data2[[1]], data2[[n]]}]; SetAttributes[pseudofunction, Listable]; Range[2, 4] // pseudofunction Here is my benchmark data, where data3 is normal sample of real data. data3 = Drop[ExcelWorkBook[[Column1 ;; Column4]], None, 1]; data2 = {a #, b #, c #, d #} & /@ Range[1, 10^5]; data = RandomReal[{0, 1}, {10^6, 4}]; Here is my benchmark code kptnw[list_] := Transpose[{Table[First@#, {Length@# - 1}], Rest@#}, {3, 1, 2}] &@list kptnw2[list_] := Transpose[{ConstantArray[First@#, Length@# - 1], Rest@#}, {3, 1, 2}] &@list OleksandrR[list_] := Flatten[Outer[List, List@First[list], Rest[list], 1], {{2}, {1, 4}}] paradox2[list_] := Partition[Riffle[list[[1]], #], 2] & /@ Drop[list, 1] RM[list_] := FoldList[Transpose[{First@li...
Post a Comment
Read more

mathematical optimization - Minimizing using indices, error: Part::pkspec1: The expression cannot be used as a part specification

I want to use Minimize where the variables to minimize are indices pointing into an array. Here a MWE that hopefully shows what my problem is. vars = u@# & /@ Range[3]; cons = Flatten@ { Table[(u[j] != #) & /@ vars[[j + 1 ;; -1]], {j, 1, 3 - 1}], 1 vec1 = {1, 2, 3}; vec2 = {1, 2, 3}; Minimize[{Total@((vec1[[#]] - vec2[[u[#]]])^2 & /@ Range[1, 3]), cons}, vars, Integers] The error I get: Part::pkspec1: The expression u[1] cannot be used as a part specification. >> Answer Ok, it seems that one can get around Mathematica trying to evaluate vec2[[u[1]]] too early by using the function Indexed[vec2,u[1]] . The working MWE would then look like the following: vars = u@# & /@ Range[3]; cons = Flatten@{ Table[(u[j] != #) & /@ vars[[j + 1 ;; -1]], {j, 1, 3 - 1}], 1 vec1 = {1, 2, 3}; vec2 = {1, 2, 3}; NMinimize[ {Total@((vec1[[#]] - Indexed[vec2, u[#]])^2 & /@ R...
Post a Comment
Read more