How can I numerically evaluate a non-trivial integral in MuPAD?
3 views (last 30 days)
Show older comments
Tommaso Isolabella
on 2 Sep 2015
Commented: Torsten
on 4 Sep 2015
I am trying to evaluate an integral in MuPAD using numeric::int(), numeric::quadrature and float(hold(int)) but none of these methods seems to work. The function I want to integrate is a parametric function defined as (h and c1 are numbers I defined before)
Psi: (M_p,R_p,x)-> -(h^3/c1)*(M_p/R_p)*(cos(x)/(R_p^2*(1-`ϵ`*cos(x))^2-h^2))
I then put the parameters I want into the function, leaving x as the only variable (thus defining a new function):
psi_(venus): x->Psi(4.867*10^24,108*10^9,x)
and then when I try to integrate it, with x from 0 to 2*PI, MuPAD doesn't evaluate it with any of the routines listed above.
Thanks in advance to whoever can help me solve this problem.
Tommaso
3 Comments
Torsten
on 3 Sep 2015
No poles in between 0 and 2*pi ?
Then it seems that your integral is equal to zero:
Best wishes
Torsten.
Accepted Answer
Walter Roberson
on 4 Sep 2015
You wrote, and I quote, "h and c1 are numbers I defined before". Nothing about epsilon there.
You also wrote "I then put the parameters I want into the function, leaving x as the only variable (thus defining a new function):"
That tells us that you do not intend anything other than x to be a free variable, but it does not tell us that epsilon has been given a value.
The values of h and epsilon are important, as they determine whether there are any poles in the integration. With the values you gave for the parameters, there is a pole when the denominator is 0, and that occurs when
x = pi - arccos((1/108000000000)*(h-108000000000)/epsilon)
x = arccos((1/108000000000)*(h+108000000000)/epsilon)
both of which are potentially in the range 0 to 2*pi.
We can insure that there are no real-valued roots by making the arccos() complex, which would be the case if abs() of the interior expression is greater than 1 since arccos() of a value outside -1 to +1 is complex. For any given h, the epsilon that are on the border are +/- ((1/108000000000)*h-1) and +/- (1/108000000000)*h+1 : epsilon larger in absolute value give you problems.
If there is a pole or two, and we have not been given reason to know otherwise, then the integral involves a division by 0 and becomes either infinite or undefined depending on circumstances.
2 Comments
Torsten
on 4 Sep 2015
I'd suggest you plot your integrand in the range 0<=x<=2*pi.
This will usually show where the problem for the integration stems from.
As I said: if the graph looks "regular", the integral will be zero.
Best wishes
Torsten.
More Answers (0)
See Also
Products
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!