1

As the title says, I encounter an irritating problem which has to do with the accuracy or the precision used in the calculations. To be more specific, I have the following code:

Clear["Global`"];
Off[General::spell];

Vnuc[x_, y_, z_] := -Mn/Sqrt[x^2 + y^2 + z^2 + cn^2] /. {Mn -> 0.08, 
cn -> 0.25};
Vdisk[x_, y_, z_] := -Md/Sqrt[x^2 + y^2 + (α + Sqrt[h^2 + z^2])^2] /. {Md
-> 0.82, α -> 3, h -> 0.1};
Vbar[G_, Mb_, a_, b_, c_, x_, y_, z_] := Module[{},
 m[u_] := Sqrt[x^2/(a^2 + u) + y^2/(b^2 + u) + z^2/(c^2 + u)];
 λ = If[m[0] > 1, u /. FindRoot[m[u]^2 == 1, {u, 1}], 0];
 Δ[u_] := Sqrt[(a^2 + u) (b^2 + u) (c^2 + u)];
 w000 = NIntegrate[1/Δ[u], {u, λ, \[Infinity]}];
 w100 = NIntegrate[1/Δ[u] 1/(a^2 + u), {u, λ, \[Infinity]}];
 w001 = NIntegrate[1/Δ[u] 1/(c^2 + u), {u, λ, \[Infinity]}];
 w010 = 2/Δ[λ] - w100 - w001;
 w110 = (w010 - w100)/(a^2 - b^2); w011 = (w001 - w010)/(b^2 - c^2);
 w101 = (w100 - w001)/(c^2 - a^2);
 w200 = 1/3 (2/(Δ[λ] (a^2 + λ)) - w110 - w101);
 w020 = 1/3 (2/(Δ[λ] (b^2 + λ)) - w011 - w110);
 w002 = 1/3 (2/(Δ[λ] (c^2 + λ)) - w101 - w011);
 w111 = (w110 - w011)/(c^2 - a^2);
 w120 = (w020 - w110)/(a^2 - b^2); w012 = (w002 - w011)/(b^2 - c^2);
 w201 = (w200 - w101)/(c^2 - a^2);
 w210 = (w110 - w200)/(a^2 - b^2); w021 = (w011 - w020)/(b^2 - c^2);
 w102 = (w101 - w002)/(c^2 - a^2);
 w300 = 1/5 (2/(Δ[λ] (a^2 + λ)^2) - w210 - w201);
 w030 = 1/5 (2/(Δ[λ] (b^2 + λ)^2) - w021 - w120);
 w003 = 1/5 (2/(Δ[λ] (c^2 + λ)^2) - w102 - w012);
 cc = 15/16 G Mb;
 pot = -(cc/6) (w000 - 6 x^2 y^2 z^2 w111 + 
   x^2 (x^2 (3 w200 - x^2 w300) + 
      3 (y^2 (2 w110 - y^2 w120 - x^2 w210) - w100)) + 
   y^2 (y^2 (3 w020 - y^2 w030) + 
      3 (z^2 (2 w011 - z^2 w012 - y^2 w021) - w010)) + 
   z^2 (z^2 (3 w002 - z^2 w003) + 
      3 (x^2 (2 w101 - x^2 w201 - z^2 w102) - w001)));
 pot
];

data1 = Table[Vnuc[x, y, 0] + Vdisk[x, y, 0] + 
Vbar[1, 0.1, 7, 1.5, 0.6, x, y, 0], {x, -10, 10, 0.5}, {y, -10, 
10, 0.5}];
data2 = Table[Vnuc[x, 0, z] + Vdisk[x, 0, z] + 
Vbar[1, 0.1, 7, 1.5, 0.6, x, 0, z], {x, -10, 10, 0.5}, {z, -10, 
10, 0.5}];
data3 = Table[Vnuc[0, y, z] + Vdisk[0, y, z] + 
Vbar[1, 0.1, 7, 1.5, 0.6, 0, y, z], {y, -10, 10, 0.5}, {z, -10, 
10, 0.5}];

Sxy = ListContourPlot[data1, ContourStyle -> Black, 
ContourShading -> False, DataRange -> {{-10, 10}, {-10, 10}}]
Sxz = ListContourPlot[data2, ContourStyle -> Black, 
ContourShading -> False, DataRange -> {{-10, 10}, {-10, 10}}]
Syz = ListContourPlot[data3, ContourStyle -> Black, 
ContourShading -> False, DataRange -> {{-10, 10}, {-10, 10}}]

However, when I try to create the three data tables data1, data2 and data3 the program reports several error messages regarding the precision on the FindRoot inside the module. Even if I tested a lot of variations for MachinePrecision, WorkingPrecision, PrecisionGoal and AccuracyGoal I could not solve this problem. Is there any other mistake that eludes me?

I would be really very grateful, if someone could help me out with this issue. Many thanks in advance.

Vaggelis_Z
  • 8,740
  • 6
  • 34
  • 79
  • 1
    Make all inputs exact numbers e.g. 82/100 instead of 0.82. Then use option settings WorkingPrecision -> 50, AccuracyGoal -> 6, PrecisionGoal -> 6. This will help for data1 and data2. For data3 it looks like you might need to consider nondefault method settings for the NIntegrate uses. – Daniel Lichtblau Mar 05 '13 at 23:27

1 Answers1

4

If you Rationalize all your finite precision numbers and add WorkingPrecision -> 50 to FindRoot the messages go away. For example:

/. {Mn -> Rationalize[0.08], cn -> Rationalize[0.25]}

and:

Table[Vnuc[x, y, 0] + Vdisk[x, y, 0] + 
  Vbar[1, Rationalize[0.1], 7, Rationalize[1.5], Rationalize[0.6], x, y, 0], {x, -10, 10, 
  1/2}, {y, -10, 10, 1/2}]

See these posts for more information:

Confused by (apparent) inconsistent precision

Adding precision for the calculation of a function

Funny behaviour when plotting a polynomial of high degree and large coefficients

Increasing number of decimal places in Manipulate input

Community
  • 1
Mr.Wizard
  • 271,378
  • 34
  • 587
  • 1,371
  • That worked only for tables data1 and data2. Table data3 reports errors regarding NItegrate this time. – Vaggelis_Z Mar 05 '13 at 22:21
  • @Vaggelis_Z did you try increasing WorkingPrecision for NIntegrate as well? – Mr.Wizard Mar 05 '13 at 22:22
  • If I add WorkingPrecision -> 50 to the three NItegrate then errors appear in all three tables. – Vaggelis_Z Mar 05 '13 at 22:25
  • @Vaggelis_Z Okay, I guess there is more going on here than simply imprecise input. I don't have time to exhaustively work through this code and I'm not that good with these particular tools anyway. I hope someone else has an answer for you. In the mean time you might search for the names of the messages you are getting to see if a similar situation has a solution you can use. – Mr.Wizard Mar 05 '13 at 22:35