How can I numerically evaluate a non-trivial integral in MuPAD?

3 views (last 30 days)
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
Tommaso Isolabella
Tommaso Isolabella on 3 Sep 2015
as I wrote in my post, c1,h and epsilon are constants I defined before integrating

Sign in to comment.

Accepted Answer

Walter Roberson
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
Tommaso Isolabella
Tommaso Isolabella on 4 Sep 2015
You are right; I didn't intend to be rude and I apologize if I answered inappropriately-I just didn't check my question to see if I mentioned epsilon (which I tought I did). Anyways thank you for the answer; epsilon is indeed in the allowed interval, so the problem is not there.
Torsten
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.

Sign in to comment.

More Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!