0
$\begingroup$

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 ϕ.

$\endgroup$
6
  • 1
    $\begingroup$ You haver not shown the NIntegrate expression that you are used. We cannot reproduce your problem without all of the relevant code. $\endgroup$
    – Bob Hanlon
    Jul 21, 2018 at 17:33
  • $\begingroup$ I have now included it. $\endgroup$
    – Morgan
    Jul 21, 2018 at 17:46
  • 1
    $\begingroup$ What happens when you replace all your machine 0.0, 1.0 etc. with exact 0, 1, etc.? $\endgroup$
    – John Doty
    Jul 21, 2018 at 18:00
  • 2
    $\begingroup$ As @JohnDoty said, machine precision numbers (e.g., 0.0, 1.0, 2.0) "contaminate" the precision of everything they touch. Either use exact numbers, high precision arbitrary precision numbers, or Rationalize the numbers. You have not defined \[Alpha] so the numeric integration cannot be performed. $\endgroup$
    – Bob Hanlon
    Jul 21, 2018 at 18:02
  • $\begingroup$ I actually changed them from 1 to 1.0 etc. in an attempt to fix the problem. $\endgroup$
    – Morgan
    Jul 21, 2018 at 18:06

3 Answers 3

3
$\begingroup$

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.

$\endgroup$
1
$\begingroup$

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]

enter image description here

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.

$\endgroup$
2
  • $\begingroup$ It gives the correct answer up to 10^-12, where it diverges from the true value. $\endgroup$
    – Morgan
    Jul 21, 2018 at 19:28
  • $\begingroup$ You should probably set WorkingPrecision to twice the precision of the desired result. See also PrecisionGoal. $\endgroup$
    – John Doty
    Jul 22, 2018 at 1:37
0
$\begingroup$

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
$\endgroup$
10
  • $\begingroup$ Okay, so the problem is that my function is more complicated than that, and involves constants, etc. When I do NIntegrate with WorkingPrecision->100, all I get is "The precision of the argument function [...] is less than WorkingPrecision (100.`)." $\endgroup$
    – Morgan
    Jul 21, 2018 at 16:42
  • $\begingroup$ Can you show how you numerically define constants? $\endgroup$ Jul 21, 2018 at 17:09
  • $\begingroup$ I have edited my question to include the function I am trying to integrate numerically. $\endgroup$
    – Morgan
    Jul 21, 2018 at 17:25
  • $\begingroup$ How the matrix H is defined? H=R? $\endgroup$ Jul 21, 2018 at 17:38
  • $\begingroup$ I have defined it better now, thanks. $\endgroup$
    – Morgan
    Jul 21, 2018 at 17:48

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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