6
$\begingroup$

I am facing an approximation issue during a plot from numerical integration through NIntegrate. The code is:

zz[r_, r0_, \[Omega]_] = ((r0/r)^(2 \[Omega]) - 1)^(-(1/2)); 
r0 = 1;
\[Omega] = -0.5;

f111[x_, r0_, \[Omega]_] := NIntegrate[zz[r, r0, \[Omega]], {r, 0, x}]

p111 = Plot[{Re[f111[x, r0, \[Omega]]]}, {x, 0, 5}] // Quiet

The plot looks like this: Numerical Plot

Although the integration can be solved analytically in Mathematica, and the code and plot looks like this:

In[1]:= \[Integral]((r0/r)^(2 \[Omega]) - 1)^(-(1/2)) \[DifferentialD]r
Out[2]= (r Sqrt[1 - (r0/r)^(2 \[Omega])]
  Hypergeometric2F1[1/2, -(1/(2 \[Omega])), 
  1 - 1/(2 \[Omega]), (r0/r)^(2 \[Omega])])/Sqrt[-1 + (r0/r)^(
 2 \[Omega])]

In[3]:= A[r_, r0_, \[Omega]_] := (
  r Sqrt[1 - (r0/r)^(2 \[Omega])]
    Hypergeometric2F1[1/2, -(1/(2 \[Omega])), 
    1 - 1/(2 \[Omega]), (r0/r)^(2 \[Omega])])/
  Sqrt[-1 + (r0/r)^(2 \[Omega])];
r0 = 1;
\[Omega] = -0.5;
Plot[Re[A[r, r0, \[Omega]]], {r, 0, 5}]

Analytical solution plot

So, there must present some errors in the numerical solution, and it must be avoided by approximation. Can anyone help with the code?

$\endgroup$
4
  • $\begingroup$ I have updated the codes. $\endgroup$
    – A_Dutta
    Jan 10, 2023 at 7:36
  • $\begingroup$ Adding the Method->"LocalAdaptive" option for NIntegrate alleviates the numerical inaccuracies on my machine. $\endgroup$ Jan 10, 2023 at 7:46
  • $\begingroup$ Can you please show me the code? @JulienKluge $\endgroup$
    – A_Dutta
    Jan 10, 2023 at 7:48
  • $\begingroup$ That's just adding an option. Include in the list of arguments to NIntegrate, e.g. NIntegrate[func, range, "Method" -> "LocalAdaptive"] $\endgroup$
    – user170231
    Jan 10, 2023 at 16:15

2 Answers 2

8
$\begingroup$

Exclude the poles. And do not need Quite anymore also. V 13.2 on windows. Basically you have an improper integral due to discontinuity at $x=1$

ClearAll["Global`*"]
zz[r_, r0_, ω_] := 1/Sqrt[(r0/r)^(2*ω) - 1]
r0 = 1;
ω = -1/2;
f[x_, r0_, ω_] := NIntegrate[zz[r, r0, ω], {r, 0, x}, Exclusions ->{1}]
Plot[Re[f[x, r0, ω]], {x, 0, 5}] 

Mathematica graphics

You can see the problem like this

Mathematica graphics

This seems to have caused NIntegrate a problem as it hit on $x=1$ value in the domain of integration.

$\endgroup$
0
5
$\begingroup$

Alternative: Use NDSolve

F = NDSolveValue[{int'[r] == 
Simplify[ComplexExpand[Re[zz[r, r0, \[Omega]]]], r > 0]
, int[0] == 0}, int, {r, 0, x} , Method -> "StiffnessSwitching" ]
Plot[ F [x]  , {x, 0, 5} ]

enter image description here

or even simpler

Derivative[-1][Function[r,Simplify[ComplexExpand[Re[zz[r, r0, \[Omega]]]], r > 0]]][r]

$\begin{cases} 0 & r\leq 1 \\ 2 \sqrt{r-1} & \text{True} \end{cases}$

$\endgroup$
0

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.