4
$\begingroup$

I integrate the following expression:

\begin{eqnarray} \int^{\pi-\alpha}_{\alpha} \sin{\theta}\sqrt{\cos^2{\theta} \left (a^2 x^2 - \frac{y^2}{\sin^2{\theta}}\right)} \, d\theta = \sqrt{a^2 x^2 \cos^2{\alpha}-y^2 \cot^2{\alpha}} \left(-\frac{\sqrt{2}\,y^2 \arctan{\left[\frac{\sqrt{2}\, a x \sin{\alpha}}{\sqrt{2y^2-a^2 x^2 + a^2 x^2 \cos{2 \alpha}}} \right]}}{ a x \sqrt{2y^2-a^2 x^2 + a^2 x^2 \cos{2 \alpha}} }- \sin{\alpha}\right) \tan{\alpha} \end{eqnarray}

I checked the result by simply differentiating the right hand side. Here, $a$, $x$ and $y$ are all constant and positive. The above result badly disagrees with the numerical integration. For instance:

$f(x=10, y=1, a= 150, \alpha=\pi/5)_{exact} = -518.232 - 0.0010472i\\ f(x=10, y=1, a= 150,\alpha=\pi/5)_{numeric} = 981.7623915187304$

I could not figure out why. Any help is appreciated. The code that I used is :

Integrate[Sqrt[Cos[θ]^2 (a^2 x^2 - y^2/Sin[θ]^2)] Sin[θ], θ]

Using the resultant expression, I took the limits $\theta \rightarrow \alpha$ and $\theta \rightarrow \pi-\alpha$ and used the simplify command:

Simplify[Limit[1/2 Sqrt[
a^2 x^2 Cos[θ]^2 - 
 y^2 Cot[θ]^2] ((Sqrt[2] y^2 ArcTan[(Sqrt[2] Sqrt[a^2 x^2] Sin[θ])/
    Sqrt[-a^2 x^2 + 2 y^2 + a^2 x^2 Cos[2 θ]]])/(
  Sqrt[a^2 x^2] Sqrt[-a^2 x^2 + 2 y^2 + a^2 x^2 Cos[2 θ]]) + 
  Sin[θ]) Tan[θ], θ -> π - α] - Limit[1/2 Sqrt[a^2 x^2 Cos[θ]^2 - 
 y^2 Cot[θ]^2] ((Sqrt[2] y^2 ArcTan[(Sqrt[2] Sqrt[a^2 x^2] Sin[θ])/
    Sqrt[-a^2 x^2 + 2 y^2 + a^2 x^2 Cos[2 θ]]])/(
  Sqrt[a^2 x^2]
    Sqrt[-a^2 x^2 + 2 y^2 + a^2 x^2 Cos[2 θ]]) + 
  Sin[θ]) Tan[θ], θ -> α]]

For the numerical integration the code is :

f2[a_, x_, y_, α_, ag_, pg_, wp_, mc_] := NIntegrate[ Sqrt[Cos[θ]^2 (a^2 x^2 - y^2/Sin[θ]^2)] Sin[θ], {θ, α, π - α}, AccuracyGoal -> ag, PrecisionGoal -> pg, WorkingPrecision -> wp, MaxRecursion -> mc]
$\endgroup$
  • 1
    $\begingroup$ Can you include the Mathematica code you used to evaluate your integral, so that potential helpers don't need to retype all that? $\endgroup$ – J. M. is away Mar 14 '18 at 15:40
  • $\begingroup$ I am sorry, added. $\endgroup$ – user91411 Mar 14 '18 at 15:49
  • $\begingroup$ @MariuszIwaniuk Yes, the code you wrote generates the same result with nintegrate but the analytic expression yields a different value. $\endgroup$ – user91411 Mar 14 '18 at 18:21
  • $\begingroup$ This looks to be yet another issue related to an antiderivative that is discontinuous within its interval of integration. Did you try plotting the antiderivative you got over the integration interval you have? $\endgroup$ – J. M. is away Mar 15 '18 at 1:31
  • $\begingroup$ @ J. M. You are right, it's deja vu. This one though is a slightly different variation on the theme because it has 6 discontinuity regions, as you see in my answer post below. $\endgroup$ – Vixillator Mar 15 '18 at 8:27
9
$\begingroup$

I apologize that the answer is a bit long because I have included a step-by-step explanation of what is going on, since a similar question was asked earlier and more might be asked on the same subject, so I thought it might be useful to post a reasonably complete answer.

This answer formulates a procedure for handling the symbolic indefinite integration of a multivariate improper integral over a domain in which the integrand contains one or more discontinuities.

(* Integrand *)
f1[\[Theta]_] :=  Sqrt[Cos[\[Theta]]^2 (a^2 x^2 - y^2/Sin[\[Theta]]^2)] Sin[\[Theta]];

Let us start by plotting the function f1:

Plot3D[f1[\[Theta]] /. {y -> 1, a -> 150}, {\[Theta], 0, 2 Pi}, {x, -10, 10}] (* Plot 1a *)

Plot 1a Plot 1a

DensityPlot[f1[θ]/.{y->1,a->150},{θ,0,2 Pi},{x,-10,10},GridLines->Automatic,ColorFunction->Function[{z},ColorData["RedGreenSplit"][ArcSin[z]]],ColorFunctionScaling->False] (*Plot 1b*)

Plot 1b Plot 1b

Plots 1a and 1b above show that there are 6 regions (w.r.t variable x) separated by discontinuities at [Theta]={Pi/2, 3 Pi/2} and x={0}. The discontinuity at [Theta]=Pi can be ignored here because the first derivative approaches the same value from both directions.

Plot3D[f1[\[Theta]] /. {x -> 10, a -> 150}, {\[Theta], 0, 2 Pi}, {y, -10,  10}] (* Plot 2a *)

Plot 2a Plot 2a

DensityPlot[f1[θ]/.{x->10,a->150},{θ,0,2 Pi},{y,-10,10},GridLines->Automatic,ColorFunction->Function[{z},ColorData["RedGreenSplit"][ArcSin[z]]],ColorFunctionScaling->False]  (*Plot 2b*)

Plot 2b Plot 2b

Plots 2a and 2b above show that there are 3 regions (w.r.t variable y) separated by discontinuities at [Theta]={Pi/2, 3 Pi/2}. No discontinuities orthogonal to the y axis. The discontinuity at [Theta]=Pi can be again ignored here because the first derivative approaches the same value from both directions.

Combining the regions from Plots 1a, 1b, 2a, and 2b, we conclude that there are 6 regions bounded by discontinuities in the domain of integration of the integrand, as shown in Plot 3 below. We will need to compute the indefinite integrals of each region separately, and then use the appropriate indefinite integral to compute the numerical value of integration for the region of interest, based on the values of the variables and constant parameters.

(* Logical structure of the 6 regions *)

R1 = {0 < \[Theta] < Pi/2 && x < 0};
R2 = {Pi/2 < \[Theta] < 3 Pi/2 && x < 0};
R3 = {3 Pi/2 < \[Theta] < 2 Pi && x <= 0};
R4 = {0 < \[Theta] < Pi/2 && x > 0}; R5 = {Pi/2 < \[Theta] < 3 Pi/2 && x > 0}; R6 = {3 Pi/2 < \[Theta] < 2 Pi && x > 0};

(* Let us show the numbering scheme of the 6 regions *)

legend = {"R1:{0<\[Theta]<Pi/2&&x<0}", "R2:{Pi/2<\[Theta]<3Pi/2&&x<0}","R3:{3 Pi/2<\[Theta]<2Pi&&x<0}", "R4:{0<\[Theta]<Pi/2&&x>0}","R5:{Pi/2<\[Theta]<3Pi/2&&x>0}", "R6:{3 Pi/2<\[Theta]<2Pi&&0x>0}"};

RegionPlot[{R1, R2, R3, R4, R5, R6}, {\[Theta], 0, 2 Pi}, {x, -10, 10}, GridLines -> {{Pi/2, 3 Pi/2}, {0}}, PlotLegends -> legend] (* Plot 3 *)

Plot 3

Plot 3

(* Indefinite integral for each region *)

IR1 = Integrate[f1[\[Theta]], \[Theta], Assumptions -> R1 && {y > 0, a > 0}];
IR2 = Integrate[f1[\[Theta]], \[Theta], Assumptions -> R2 && {y > 0, a > 0}];
IR3 = Integrate[f1[\[Theta]], \[Theta], Assumptions -> R3 && {y > 0, a > 0}];
IR4 = Integrate[f1[\[Theta]], \[Theta], Assumptions -> R4 && {y > 0, a > 0}];
IR5 = Integrate[f1[\[Theta]], \[Theta], Assumptions -> R5 && {y > 0, a > 0}];
IR6 = Integrate[f1[\[Theta]], \[Theta], Assumptions -> R6 && {y > 0, a > 0}];

(* Check to see if any are unique *)
IR1 == IR2 == IR3 == IR4 == IR5 == IR6 (* They are all the same! Note that for some other integrand functions each reagion may yield a different expressions for the indefinite integral *)

OUT: True

(* Since all expressions above are identical, let us use the expression of the indefinite integral for the entire domain *)
IR = Integrate[f1[\[Theta]], \[Theta]];

Let us now compute the indefinite integrals (in [Theta]) using the given parameters: f(x=10,y=1,a=150,[Alpha]=[Pi]/5) These parameters imply the integration limits of [Theta] from pi/5 to 4 Pi/5, which has a first part that is in Region 4 and a second part that is in Region 5. Therefore we need to evaluate the numerical results from the symbolic integrals in two parts, PIECEWISE WHILE HONORING THE DISCONTINUITY BOUNDARIES:

FirstPartInRegion4 = IR /. {x -> 10, y -> 1, a -> 150, \[Theta] -> {Pi/2 - Pi/5}} // N  (* Region 4 upper boundary is at \[Theta]=Pi/2 *)

OUT: {490.879 + 0.000523599 I}

SecondPartInRegion5 = IR /. {x -> 10, y -> 1, a -> 150, \[Theta] -> {4 Pi/5 - Pi/2}} // N(*  Region 5 lower boundary is at \[Theta]=Pi/2 *)

OUT: {490.879 + 0.000523599 I}

FullPart = FirstPartInRegion4 + SecondPartInRegion5 (* This sum covers the entire integration limits of the test parameters *)

OUT: {981.757 + 0.0010472 I} This is the correct value!

Let us now try an incorrect procedure which ignores the region boundaries, and try to compute the numerical value in one part:

(* Let us now try an incorrect procedure which ignores the region boundaries, and try to compute the numerical value in one part *)

IR /. {x -> 10, y -> 1, a -> 150, \[Theta] -> {4 Pi/5 - Pi/5}} // N (* This gives the wrong value *)

OUT: {-678.379 - 0.000523599 I} This is wrong!

(* Let us compare with the numerical definite integral *)

f2 = NIntegrate[f1[\[Theta]] /. {x -> 10, y -> 1, a -> 150}, {\[Theta], Pi/5, Pi - Pi/5}]

OUT: 981.762

(* We also get the same answer from the symbolic form of the definite integral, with assigned variables *)
f2 = Integrate[f1[\[Theta]] /. {x -> 10, y -> 1, a -> 150}, {\[Theta], Pi/5, Pi - Pi/5}] //N

OUT: 981.762 + 0. I

CONCLUSION: When using symbolic indefinite integration of an integrand function that is discontinuous, it is advisable to divide the domain of integration into regions that are free of discontinuities, and apply the indefinite integration strictly within the boundaries of each region.

$\endgroup$
  • $\begingroup$ Wow. I was not expecting such a detailed answer. Thank you. The antiderivative has branch cuts and that was the issue. Next time I will keep that in mind. $\endgroup$ – user91411 Mar 15 '18 at 8:34
  • $\begingroup$ You are welcome. Thank you for posting your question. Cheers! $\endgroup$ – Vixillator Mar 15 '18 at 8:41
6
$\begingroup$
  1. There's no point in using Limit since you're not actually extending to the singular boundary but integrating to regular interior points. /. would work as well.
  2. The integrand has a cusp singularity at Pi/2. The antiderivative need not be, and in fact is not, continuous. The Fundamental Theorem of Calculus does not apply.
  3. Even for nice functions the analytic expression may be discontinuous. See the first Possible Issues example under Definite Integrals on ref/Integrate.

The correct way to do the symbolic integral is as follows:

Block[{x = 10, y = 1, a = 150, \[Alpha] = Pi/5, ad} 
   ad = Integrate[Sqrt[Cos[θ]^2 (a^2 x^2 - y^2/Sin[θ]^2)] Sin[θ], θ];
   N[ (ad /. θ -> Pi - θ) 
      - Limit[ad, θ -> Pi/2, Direction -> "FromAbove"]
      + Limit[ad, θ -> Pi/2, Direction -> "FromBelow"]
      - (ad /. θ -> \[Alpha])
   ]
]
(* 981.762 + 2.1684*10^-19 I *)

To see point (2), evaluated the following

Block[{x = 10, y = 1, a = 150, \[Alpha] = Pi/5, ad, f},
       f = Sqrt[Cos[θ]^2 (a^2 x^2 - y^2/Sin[θ]^2)] Sin[θ];
       ad = Integrate[f, θ]; 
       Plot[{f, Re@ad}, {θ, 0, Pi}]
   ]
$\endgroup$
  • $\begingroup$ Thanks for the illustrative answer. $\endgroup$ – user91411 Mar 15 '18 at 8:36

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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