Numerical problems with calculation of function

Question

I have a problem with a calculation which I try to do with Mathematica (8.0.1).

I have the following function which describes a cone (with half opening angle theta) which can be tilted (with the angles thetaTilt and phiTilt) around its base point lG. For some unimportant reasons the tip of this cone (if not tilted) is at the position –rout/Tan[theta]. So, rout describes in the following function the diameter of the cone tip at the position z=0 (again for the not tilted case):

coneFunction[x_, y_, theta_, thetaTilt_, phiTilt_, lG_, rout_] := 
1/(Cos[2 theta] - 
Cos[2 thetaTilt]) (rout Cos[thetaTilt] Sin[2 theta] + 
 2 lG (1 + Cos[2 theta] + 2 Cos[thetaTilt]) Sin[thetaTilt/
   2]^2 - \[Sqrt](Cos[theta]^2 (lG^2 + x^2 + y^2 + 
      rout Cot[theta] (2 lG + rout Cot[theta]) +  2 x y Sin[2 phiTilt] + 
      Cos[2 thetaTilt] (-lG^2 + x^2 + y^2 -  rout Cot[theta]
 (2 lG + rout Cot[theta]) -  2 x y Sin[2 phiTilt]) + 
      8 (x Cos[phiTilt] + y Sin[phiTilt]) Sin[theta] (rout Cos[theta] +
lG Sin[theta]) Sin[thetaTilt] + 
      2 (x - y) (x + y) Cos[2 phiTilt] Sin[thetaTilt]^2 + 
      Cos[2 theta] (-lG^2 - 2 (x^2 + y^2) + 
         lG^2 Cos[2 thetaTilt] - 2 rout Cot[theta] (2 lG + rout Cot[theta])
 Sin[thetaTilt]^2))) + (x Cos[phiTilt] + y Sin[phiTilt]) Sin[2 thetaTilt]);

Since it is difficult to understand what this function does (and what I mean in my text), I suggest that one plots an cross-section of this function, e.g. with

theta=3.5*10^-3;(*half opening angle of cone in rad*)
rout=0.65*10^-6;(*diameter at tip of cone*)
lG=10*10^-3;(*basepoint of the cone*)
Plot[coneFunction[x, 0, theta, 0, 0, lG, rout], 
{x, -0.000005,0.000005}, 
PlotRange -> {1.2*{-0.000005, 0.000005}, {-0.0002, 0.001}}, 
PerformanceGoal -> "Quality", PlotPoints -> 120]

Or a 3D plot with e.g.:

thetaTilt = 0;
phiTilt = 0;
Show[Plot3D[{coneFunction[x, y, theta, thetaTilt, phiTilt, lG, rout]}, 
{x, -0.001, 0.001}, {y, -0.001, 0.001}, 
Exclusions -> None, PlotPoints -> 60, PerformanceGoal -> "Quality", 
Mesh -> False, PlotStyle -> Directive[Opacity[0.4], Gray]],
Axes -> True, AxesLabel -> {"x/m", "y/m", "z/m"}, 
BoxRatios -> {1, 1, 1}, 
PlotRange -> {{-0.0001, 0.0001}, {-0.0001, 0.0001}, {-0.001, 0.012}}]

If one sets thetaTilt to small values such as 10^-4, it can be seen easily, that the cone is being tilted by this value.

If the cone is not tilted, the function is independent of lG. This can be seen easily by:

Simplify[coneFunction[x, y, theta, 0, 0, lG, rout]]

So, regardless what x or y values I chose, I should get back always the same result for different values of the base point lG.

But that is not the case! If I use e.g.

testTable = 
Table[SetPrecision[
coneFunction[-0.000005, 0, theta, 0, 0, lG, rout], 10], 
{lG, 0.0053, 0.0055, 0.0000001}];

I get as output most time 0.001242850670, but also some other values (e.g. 0.001242846719). So, if I have bad luck and chose a “bad” lG value, I would get back wrong results.

I assume that this has something to do with the numeric of Mathematica. I tried to get rid of this effect with some attempts (e.g. SetPrecision, Rationalize, Evaluate), but none of them worked.

So, who can give me some advice? How can I manage that my function gives back for all lG the right result (also for the tilted case)?

Since I am clueless I would appreciate any support!

Answer 1

With the current setup you get three different answers:

testTable // Union

(* ==> {0.001242846719, 0.001242850670, 0.001242854621} *)

The problem is that you haven't sufficiently increased the precision of all numbers, including those in the loop bounds.

theta = SetPrecision[3.5*10^-3, 50];
rout = SetPrecision[0.65*10^-6, 50];

testTable = 
  Table[
    SetPrecision[
      coneFunction[-0.000005`50, 0, theta, 0, 0, lG, rout], 
      10
    ], 
    {lG, 0.0053`50, 0.0055`50, 0.0000001`50}
  ];

testTable // Union

(* ==>  {0.001242852068} *)

The following might also be worth to note:

SetPrecision[Sin[1000], 10] // Precision

(* ==> 7.167413092 *)

The expression is evaluated with 10 digits of precision, but the Sin operation loses quite some precision.