Hi guys! I am trying to solve the Poisson equation in a sphere. I am doing that to check the results from the numerical calculation with NDSolveValue against the known solution in the sphere. If that works out fine, I will solve the Poisson equation in a more general spheroid or ellipsoid.

See the attached Notebook (Mathematica 11, but I don't use special features) for the solver I am using. I have two problems:

• The numerical solution does not agree completely with the analytical solution. In a previous version it used to agree in one part in 10^3, but by refining the mesh at the boundary I took it down to 1 part in 10^6. Not bad, but can I take it down to at least 1 part in 10^8 without refining the mesh to crazy levels?

• I ask for 32 digits of precision to the final result of NDSolve, and the result does not have the required precision. What is the problem? I read the manuals and tutorials, and I still need help here.

In the .nb you'll see

>> b0=1;

>> uval = NDSolveValue[{Inactive[Laplacian][u[x, y, z], {x, y, z}] == 
    Rationalize[3/b0], DirichletCondition[u[x, y, z] == f1 /. r -> 1, True]}, u, {x, y, z} \[Element] mesh, WorkingPrecision -> 32, 
  AccuracyGoal -> 32, PrecisionGoal -> 32, MaxSteps -> Infinity, InterpolationOrder -> All]

>> N[uval[0.1, 0, 0],16]

-1.49498

while th theoretical value is -1.495

Any suggestion?

Thanks

Attachments: