Precision for function

Question

I have a notebook where I have to work with high precision (say, 100). I have tried several methods to define it for the entire notebook, including

$MinPrecision = 2 $MachinePrecision

but when I define a function, say:

f[x_] := x^3;

that function does not have the required precision. Moreover, when I try to NIntegrate that function, I get a "The precision of the argument function ... is less than WorkinPrecision".

What should I do to solve this?

Edit: Here is my actual function which I am trying to numerically integrate over the sphere:

d[n_, H_, q_] := 
 Table[Sum[
   Sum[(n[[i]] - q[[i]]) (H[[j, k]] n[[j]] n[[k]])/(2 (1 - q.n)), {k, 
      1, 3}] - 1/2 H[[i, j]] n[[j]], {j, 1, 3}], {i, 1, 3}]

The numerical integration I use is

NIntegrate[
 Sin[θ] (d[n, H, q].{1.0, 0.0, 
     0.0}) (d[m, H, q].{Cos[α], 
     0.0, -Sin[α]}), {θ, 0.0, π}, {ϕ, 
  0.0, 2.0 π}, WorkingPrecision -> 20]

Here, n and q are vectors:

n = {0.0, 0.0, 1.0};
m = {Sin[α], 0.0, Cos[α]};
q = {Sin[θ] Cos[ϕ], Sin[θ] Sin[ϕ], Cos[θ]};

and R is a matrix

R = Simplify[
   Inverse[RotationMatrix[θ, {0, 1, 0}]].Inverse[
     RotationMatrix[ϕ, {0, 0, 1}]]];

H = Transpose[R].( {
     {1.0, 1.0, 0.0},
     {1.0, -1.0, 0.0},
     {0.0, 0.0, 0.0}
    } ).R;

The integration is over θ and ϕ.

Answer 1

The precision of an exact function like f[x_] := x^3 depends on its input:

Precision[f[2]]
(* \[Infinity] *)

Control the input precision (3 digits here) to get a controlled precision result:

Precision[f[2`3]]
(* 2.52288 *)

A number with a decimal point and no precision specification is a "machine" number. Machine numbers are special: unless you explicitly prevent them, they infect calculations with their precision and disable precision tracking. Here, the machine constant 1.5 forces calculations to be in machine precision:

g[x_] := x^3 + 1.5
Precision[g[2`3]]
(* MachinePrecision *)

If you want a function to yield controlled precision output, avoid machine constants in its definition. It's best to use exact constants: 1/3, Pi, Sqrt[2], `Root[1 - #1^2 + #1^3 &, 1]`, etc.

Answer 2

Let's see if this might help

α = Pi/4;
n = {0, 0, 1};
m = {Sin[α], 0, Cos[α]};
q = {Sin[θ] Cos[φ], Sin[θ] Sin[φ], Cos[θ]};
R = Inverse[RotationMatrix[θ, {0, 1, 0}]].Inverse[RotationMatrix[φ, {0, 0, 1}]];
H = Transpose[R].({{1, 1, 0}, {1, -1, 0}, {0, 0, 0}}).R;
d[n_, H_, q_] := Table[Sum[Sum[
  (n[[i]] - q[[i]]) (H[[j, k]] n[[j]] n[[k]])/(2 (1 - q.n)), {k, 1, 3}] -
  1/2 H[[i, j]] n[[j]], {j, 1, 3}], {i, 1, 3}];

When I have a reluctant definite integral the first thing I try is Simplify or FullSimplify on the integrand and look for denominators or oscillations or infinities. That FullSimplify gives a result which is a tiny fraction of your original unsimplified expression and is half the size of just Simplify. For particularly difficult cases you can let FullSimplify know the domain of φ and θ and see if that helps it even further reduce the size of the expression.

FullSimplify[Sin[θ] (d[n, H, q].{1, 0, 0}) (d[m, H, q].{Cos[α], 0, -Sin[α]})]

which shows me

-((Sin[θ]^2 (Cos[φ] + Sin[φ]) ((-8 Cos[θ] + Sqrt[2] (3 + Cos[2 θ])) Cos[2 φ] +
8 Cos[φ] Sin[θ] - 6 Sqrt[2] Sin[θ]^2 - 4 Sin[2 θ] Sin[φ] + 
2 (3 - 2 Sqrt[2] Cos[θ] + Cos[2 θ]) Sin[2 φ]))/(32 (-Sqrt[2] + Cos[θ] + Cos[φ] Sin[θ])))

that yes you do have denominators which might go to zero.

The next thing I do is look at a Plot

Plot3D[FullSimplify[
  Sin[θ] (d[n, H, q].{1, 0, 0}) (d[m, H, q].{Cos[α], 0, -Sin[α]})],
  {θ, 0, Pi}, {φ, 0, 2 Pi}, WorkingPrecision -> 20]

This doesn't look bad, it doesn't appear to wildly oscillate near any point, doesn't appear to go to infinity, doesn't have any gaping holes were it isn't defined. But Plot can sometimes lie to me so I look for zero denominators.

And Reduce tells me that the denominator does go to zero

N[Reduce[{32 (-Sqrt[2] + Cos[θ] + Cos[φ] Sin[θ]) ==
  0, 0 <= θ <= Pi, 0 <= φ <= 2 Pi}, {θ, θ}] // Simplify]

which happens when

θ == 0.785398 && (φ == 0. || φ == 6.28319)

And yes, θ == Pi/4 with φ == 0 or φ == 2Pi give an integrand of 1/0 and those two points appear to have been missed by Plot. You might want to explore the behavior near those two points and see what you get.

So those may be the sticking points that NIntegrate is complaining about.

But I try anyway

NIntegrate[FullSimplify[
  Sin[θ] (d[n, H, q].{1, 0, 0}) (d[m, H, q].{Cos[α], 0, -Sin[α]})],
  {θ, 0, Pi}, {φ, 0, 2 Pi}, WorkingPrecision -> 20] 

And the default method used by NIntegrate complains about it is having problems getting the result to converge, but it does finally give about 20 digits of precision.

You can try giving NIntegrate some Method options to see if you can get it to stop complaining or give a better result, but this doesn't seem to be necessary.

Now you just need to decide whether that result really is correct or not.

Answer 3

Constants must be defined exactly, should be so

    n = {0, 0, 1};
m = {Sin[α], 0, Cos[α]};
q = {Sin[θ] Cos[ϕ], Sin[θ] Sin[ϕ], 
   Cos[θ]};
R = Simplify[
   Inverse[RotationMatrix[θ, {0, 1, 0}]].Inverse[
     RotationMatrix[ϕ, {0, 0, 1}]]];
H = Transpose[R].({{1, 1, 0}, {1, -1, 0}, {0, 0, 0}}).R;
t = Table[
   Sum[Sum[(n[[i]] - 
         q[[i]]) (H[[j, k]] n[[j]] n[[k]])/(2 (1 - q.n)), {k, 1, 
       3}] - 1/2 H[[i, j]] n[[j]], {j, 1, 3}], {i, 1, 3}];
t1 = Table[
   Sum[Sum[(m[[i]] - 
         q[[i]]) (H[[j, k]] m[[j]] m[[k]])/(2 (1 - q.m)), {k, 1, 
       3}] - 1/2 H[[i, j]] m[[j]], {j, 1, 3}], {i, 1, 3}];
int = Sin[θ] (t.{1, 0, 0}) (t1.{Cos[α], 
      0, -Sin[α]});
NIntegrate[
 int /. α -> Pi/4, {θ, 0, π}, {ϕ, 0, 2 π}, 
 WorkingPrecision -> 20]
0.55395347102611116853