4

I am solving PDE by means of NDSolve and I am interested in a highly accurate value for the solution at a certain point. For some reason, the solution obtained by using NDSolveValue evaluates only to 6 digits (even if increasing WorkingPrecision and PrecisionGoal!).

You can find four very simple examples describing the issue below...

Particularly, I noticed the following behaviour: when solving an ODE by using DirichletCondition, it causes the following:

sol1 = NDSolveValue[{D[u[x], {x, 2}] == 2 , u[0] == 0, u[1] == 1},     
   u, {x, 0, 1}, WorkingPrecision -> 60, PrecisionGoal -> 30, 
   AccuracyGoal -> 30];
N[sol1[1/3], 20]

sol2 = NDSolveValue[{D[u[x], {x, 2}] == 2 , 
    DirichletCondition[u[x] == 0, x == 0], 
    DirichletCondition[u[x] == 1, x == 1]}, u, {x, 0, 1}, 
   WorkingPrecision -> 60, PrecisionGoal -> 30, AccuracyGoal -> 30];
N[sol2[1/3], 20]

which returns as a result:

0.11111111111111111111
0.111111

While for a PDE, this does not make any difference, obtaining more than 6 digits does not seem possible ??:

sol4 = NDSolveValue[{D[u[x, y], {x, 2}] + D[u[x, y], {y, 2}] - 
      2 u[x, y] == 0, u[0, y] == 0, u[x, 1] == 1}, 
   u, {x, 0, 1}, {y, 0, 1}, WorkingPrecision -> 60, 
   PrecisionGoal -> 30, AccuracyGoal -> 30];
N[sol4[1, 0], 20]

sol3 = NDSolveValue[{D[u[x, y], {x, 2}] + D[u[x, y], {y, 2}] - 
      2 u[x, y] == 0, DirichletCondition[u[x, y] == 0, x == 0], 
    DirichletCondition[u[x, y] == 1, y == 1]}, 
   u, {x, 0, 1}, {y, 0, 1}, WorkingPrecision -> 60, 
   PrecisionGoal -> 30, AccuracyGoal -> 30];
N[sol3[1, 0], 20]

which returns:

0.297163
0.297163

Any help or insight would be appreciated!

AnSa
  • 41
  • 2
  • Yep, it's not possible. Mathematica's numerical methods are not designed to give 30 digits of precision. You can investigate small effects my including distorsions and solving linearized equation for them though. – Vsevolod A. Mar 15 '18 at 17:34
  • Interesting observation. FYI to view the full machine precision result do sol4[1, 0] // FullForm Also be aware the result of NDSolve is an interpolation function, so whatever result you get suffers a further loss of precision if you seek a value that doesnt happen to be a numerical grid point. – george2079 Mar 15 '18 at 18:03
  • FWIW I get the same 60 digit precision result for cases 1 and 2.. It is just the PDE that runs only machine precision. – george2079 Mar 15 '18 at 18:07
  • @george2079 : Do I understand you correctly, you obtain a 60 digit precision result of the function sol2[1/3] defined above ? If yes, could you please post the line of code? – AnSa Mar 18 '18 at 12:42
  • I get the same 0.11111111111111111111 for both 1 and 2 using your exact code. sol2[1/3]//Precision shows 59.xxx.... (version 10.1) You should specify your version. – george2079 Mar 18 '18 at 13:41
  • I use version 11.1.1.0 and it returns 59.011 for sol1[1/3]//Precision but MachinePrecision for sol2[1/2]//Precision. So it could be that in the newer version, the use of DirichletCondition for an ODE will indeed invoke the FEM solver. – AnSa Mar 19 '18 at 10:18

1 Answers1

3

FEM is restricted to machine precision (see, e.g., this). It is automatically invoked when you use DirichletCondition.

Mathematica graphics

Michael E2
  • 235,386
  • 17
  • 334
  • 747