Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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.

share|improve this question
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 5 '13 at 23:27

1 Answer 1

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

share|improve this answer
    
That worked only for tables data1 and data2. Table data3 reports errors regarding NItegrate this time. –  Vaggelis_Z Mar 5 '13 at 22:21
    
@Vaggelis_Z did you try increasing WorkingPrecision for NIntegrate as well? –  Mr.Wizard Mar 5 '13 at 22:22
    
If I add WorkingPrecision -> 50 to the three NItegrate then errors appear in all three tables. –  Vaggelis_Z Mar 5 '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 5 '13 at 22:35

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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