# Perform a numerical integration?

GROUPS:
 I am getting error while evaluating following integral in Mathematica. Even though I am getting result I am sure it is wrong.  ClearAll["GLOBAL'*"] {r, a, b, Z, A, B,phi } = {4.087, 1.205, 0.3812, 0, 345.0527, 606741.04395, Pi/3} // Rationalize[#, 0] &; JJ[n_] := r NIntegrate[1/((t Sin[phi])^2 + r^2 + (t Cos[phi] + Z - z)^2 - 2 r t Sin[phi] Cos[eta])^ n, {eta, 0, 2 Pi}, {z, 0, 50}, {t, -5/3, 5/3}, MinRecursion -> 10, MaxRecursion -> 35, WorkingPrecision -> 20]; a b (-A JJ[3] + B JJ[6]) NIntegrate::eincr: The global error of the strategy GlobalAdaptive has increased more than 2000 times. The global error is expected to decrease monotonically after a number of integrand evaluations. Suspect one of the following: the working precision is insufficient for the specified precision goal; the integrand is highly oscillatory or it is not a (piecewise) smooth function; or the true value of the integral is 0. Increasing the value of the GlobalAdaptive option MaxErrorIncreases might lead to a convergent numerical integration. NIntegrate obtained 0.0144321840040547434887970543697143839268866630195143253550482722378718970. and 1.03407772966344851904649524874741779744326972276591906470617709498489670.*^-10 for the integral and error estimates. >> NIntegrate::eincr: The global error of the strategy GlobalAdaptive has increased more than 2000 times. The global error is expected to decrease monotonically after a number of integrand evaluations. Suspect one of the following: the working precision is insufficient for the specified precision goal; the integrand is highly oscillatory or it is not a (piecewise) smooth function; or the true value of the integral is 0. Increasing the value of the GlobalAdaptive option MaxErrorIncreases might lead to a convergent numerical integration. NIntegrate obtained 6.36390097912654357958699898191423875303498751269284441958249877991594170.*^-6 and 1.24379856161842596853385313946812168626482690745490628686126496952390770.*^-13 for the integral and error estimates. >> -2.100045775865940530 Symbolic integration doesn;t give answer in Mathematica. Is there anyway to speed up while doing symbolic integration in this case/ Attachments:
1 month ago
11 Replies
 Neil Singer 1 Vote It integrates quickly without errors if you delete the WorkingPrecision->20. Do you really need 20 digits?You can quickly and easily get 8 digits of precision without changing the WorkingPrecision. JJ[n_] := r NIntegrate[ 1/((t Sin[phi])^2 + r^2 + (t Cos[phi] + Z - z)^2 - 2 r t Sin[phi] Cos[eta])^n, {eta, 0, 2 Pi}, {z, 0, 50}, {t, -5/3, 5/3}, PrecisionGoal -> 8]; 
1 month ago
 Thank you. The answer is -2.00015 which is not matching with answer in the source from which I have taken this problem. The answer from their calculation is -4.4. I dont need 20 digits 5 or 6 is ok. but I am really confused with what is causing such a difference. Is it possible to evaluate it symbolically using Integrate command rather than numerical integration.
1 month ago
 Thank you. The answer is -2.00015 which is not matching with answer in the source from which I have taken this problem. The answer from their calculation is -4.4. I dont need 20 digits 5 or 6 is ok. but I am really confused with what is causing such a difference. Is it possible to evaluate it symbolically using Integrate command rather than numerical integration.
1 month ago
 Neil Singer 1 Vote I tried several algorithms (using method->) and get the same answer. I believe you likely have a small discrepancy with your problem source (or they have a typo). I doubt that Mathematica is wrong here because the function is smooth -- there are no numerical instabilities (such as discontinuities). I would be surprised if the answer is wrong.
1 month ago
 I am also thinking the same. But when \phi=0 I get the same value as theirs which can rule out the possibility of a typo. How you came to a conclusion that it is smooth. Did you plot the function by keeping one of the integration variable as a constant. Is it possible to do a symbolic integration rather than numerical integration. I tried it but it takes too long which forced me to kill the job. For \phi=0 the computation time is very less compared for non zero \phi's and without any error messages.
1 month ago
 Neil Singer 1 Vote The expression is of the form 1/(y^n). The expression only can blow up when y is close to zero. I solved for y= 0: NSolve[((t Sin[phi])^2 + r^2 + (t Cos[phi] + Z - z)^2 - 2 r t Sin[phi] Cos[eta]) == 0, {t, z, eta}, Reals] Which gave me three solutions: {{t -> ConditionalExpression[-4.71926, C[1] \[Element] Integers], z -> ConditionalExpression[-2.35963, C[1] \[Element] Integers], eta -> ConditionalExpression[1. (-3.14159 + 6.28319 C[1]), C[1] \[Element] Integers]}, {t -> ConditionalExpression[-4.71926, C[1] \[Element] Integers], z -> ConditionalExpression[-2.35963, C[1] \[Element] Integers], eta -> ConditionalExpression[1. (3.14159 + 6.28319 C[1]), C[1] \[Element] Integers]}, {t -> ConditionalExpression[4.71926, C[1] \[Element] Integers], z -> ConditionalExpression[2.35963, C[1] \[Element] Integers], eta -> ConditionalExpression[6.28319 C[1], C[1] \[Element] Integers]}} Note that all of the solution are outside of your integration range so unless I missed something, your expression is continuous in your region. I also fixed eta and n-> 3 and plotted the expression: Plot3D[1/((t Sin[phi])^2 + r^2 + (t Cos[phi] + Z - z)^2 - 2 r t Sin[phi] Cos[eta])^n /. eta -> 2 Pi /. n -> 3, {z, 0, 50}, {t, -5/3, 5/3}, PlotRange -> All] It spikes to .003 and is mostly very small but the area under the curve seems calculable and smooth -- You can try different values and integrate it in a different program to check. Are you sure the original source does not have a slightly different parenthesis or a subtle difference?
1 month ago
 Thank you.I tried plotting the integrand for \phi=0 and \phi=pi/2. There is a spike for \phi=pi/2 which is missing in case of \phi=0. Is that spike causing problems. Please see attachment. Attachments:
1 month ago
 Neil Singer 1 Vote I don't think the shape would cause problems -- the functions are smooth and continuous. I was working in Matlab today so I checked your integral in Matlab and got the same answer -- you should look for a problem in your expression or the published answer -- I think Mathematica and Matlab must be numerically correct.Regards
 Neil Singer 1 Vote I looked at this one more time but I believe that the MMA numerical result is correct. You can integrate the expression in eta for a given value of n and get a closed form solution. I added assumptions based on your range of interest until the equation simplified without Conditionals and without complex terms: exp = 1/((t Sin[phi])^2 + r^2 + (t Cos[phi] + Z - z)^2 - 2 r t Sin[phi] Cos[eta])^n int1 = Integrate[ exp /. phi -> Pi/3 /. Z -> 0 /. n -> 3, {eta, 0, 2 Pi}, Assumptions -> z > 0 && z < 50 && t > 0 && r \[Element] Reals && t \[Element] Reals && z \[Element] Reals && -r^4 + r^2 (t^2 + 2 t z - 2 z^2) - (t^2 - t z + z^2)^2 < 0 && r != 0 ] The result is (\[Pi] (2 r^4+2 (t^2-t z+z^2)^2+r^2 (7 t^2-4 t z+4 z^2)))/(r^4-r^2 (t^2+2 t z-2 z^2)+(t^2-t z+z^2)^2)^(5/2) If you numerically integrate this in z and t you get the same answer. I think you can also now analytically integrate this wrt t and z. I did the indefinite integral quickly but the definite integral was taking a long time. Either way, you can plot the curve and see the area that you are trying to integrate: Plot3D[r*int1 /. MapThread[ Rule, {{r, a, b, Z, A, B, phi}, {4.087, 1.205, 0.3812, 0, 345.0527, 606741.04395, Pi/3}}], {z, 0, 50}, {t, -5/3, 5/3}, PlotRange -> All] which gives you this:Note the curve is smooth and continuous -- there is nothing there that would cause integration problems. I know you do not want to accept this answer, however, it is unlikely for MMA, Matlab, and this analysis to all be simultaneously wrong with the same wrong answer. Isn't it more likely that the paper that was published has a typographical error? If you really want to integrate the expression analytically, you can take the result that I got (which is a ratio of polynomials) and experiment some more in MMA (maybe take the indefinite integral and evaluate it yourself). Another viable option is to look up the integration formulas (for t and then z, or vice-versa). A good place to start is Gradshteyn and Ryzhik -- they will have many integrals that MMA cannot do. You can apply the formulas in MMA and analytically get your answer. I hope this helps.Regards