4
$\begingroup$

I would like to analytically integrate a function in Mathematica. However, I have found a problem with the final result when I compare it with the numerically integration of the same function. An example follows:

Tex=1/(-1509.04 + 1/(-1.64387 - 1.05604 Cos[0.88 \[Theta]]) + 1/(
 0.000622346 + 0.0000998034 Cos[0.88 \[Theta]]) +
 397.918 Cos[0.88 \[Theta]])
p1 = Plot[NIntegrate[Tex, {\[Theta], 0, x}], {x, 0, 2*Pi*25}]
Tint = Integrate[Tex, \[Theta]]

When I integrate Tint, I obtain an expression with complex numbers. Even if I try to plot the real part of this expression, this does not reproduce the plot p1. How can I correctly compute this expression?

Improve this question
$\endgroup$
1
  • $\begingroup$ What do you mean by "integrateTint"? Do you mean "integrate Tint" (Tint is a variable name)? Or something else? Please respond by editing (changing) your question, not here in comments (without "Edit:", "Update:", or similar - the question should appear as if it was written today). $\endgroup$
    – Peter Mortensen
    Dec 6, 2021 at 0:03

6 Answers 6

Reset to default
5
$\begingroup$

Consider:

Plot[Tint // Chop, {\[Theta], 0, 10}, PlotRange -> All]

enter image description here

You see that the function has a jump around 3.56. Why is this? Look at the result from "Integrate":

Integrate[Tex, \[Theta]]

enter image description here

You see that you have an "ArcTan". This is actually a multivalued function. However, MMA folds this function into the range -Pi/ to Pi/2 (for real arguments). This creates the jumps. NIntegrate does not create this jumps and gives you a smooth result.

Improve this answer
$\endgroup$
3
  • $\begingroup$ thank you very much for your answer. Is there a way to compute this integral analytically? Or the only way to do it is by resorting to numerical procedures? $\endgroup$
    – VDF
    Dec 5, 2021 at 10:41
  • 3
    $\begingroup$ There is a problem with multivalued functions, which branch do you take? In all computer algebra system I know, one defines one branch as main-branch and folds the multivalued function into this range. This means that the computer-function occasionally jumps from one branch to another. To eliminate this and get a smooth function, you would have to follow the "computer-function" and add the corresponding value if it jumps to another branch. $\endgroup$
    – Daniel Huber
    Dec 5, 2021 at 14:20
  • 1
    $\begingroup$ @VDF: To follow up on Daniel Huber's comment, see this Wolfram blog post for more information on why all this "branch" business is necessary. $\endgroup$
    – Michael Seifert
    Dec 6, 2021 at 17:03
1
$\begingroup$

Even if I try to plot the real part of this expression, this does not reproduce the plot p1

Is there a way to compute this integral analytically? Or the only way to do it is by resorting to numerical procedures?

It works for me on V 12.3.1. But it helps when using Integrate to first convert the input to exact numbers as follows (when using exact solver, best to given them exact input)

Clear["Global`*"]
tex = 1/(-1509.04 + 1/(-1.64387 - 1.05604 Cos[0.88 θ]) + 
         1/(0.000622346 + 0.0000998034 Cos[0.88 θ]) + 
         397.918 Cos[0.88 θ]);
tex = SetAccuracy[tex, Infinity];
res = Integrate[tex, θ];
Plot[Re@res, {θ, 0, 10}, PlotRange -> {Automatic, {-0.15, 0.15}}]

Mathematica graphics

See Also ReImPlot

Improve this answer
$\endgroup$
2
  • $\begingroup$ Thanks @Nasser, but your result does not reproduce the one obtained with NIntegrate. $\endgroup$
    – VDF
    Dec 5, 2021 at 10:47
  • 3
    $\begingroup$ @VDF Well, you asked for analytical solution, and that is what Mathematica gives. If you think the result from Integrate is wrong, you could submit a bug report to WRI. I do not know myself which is correct. $\endgroup$
    – Nasser
    Dec 5, 2021 at 10:48
1
$\begingroup$

We can't use the antiderivative returned by Integrate because it contains multi-valued functions with default branch-cuts and in order to use a multi-valued antiderivative, the integration path needs to be in an analytic domain of the function. In order to accomplish this, replace the ArcTan function in the RootSum antiderivative with an analytically-continuous version. No need to replace the Log expression as the integration path does not cross a default branch cut.

First integrate a rational version of the integrand:

integrand = 
  1/(-1509.04 + 1/(-1.64387 - 1.05604 Cos[0.88 \[Theta]]) + 
     1/(0.000622346 + 0.0000998034 Cos[0.88 \[Theta]]) + 
     397.918 Cos[0.88 \[Theta]]);
ClearAll[currentArcVal, theCurrentRoot, currentZ];
newIntegrand = Rationalize[integrand, 10^-16];
antiD = Integrate[newIntegrand, \[Theta]]

6250/11 RootSum[
  189275993388171865138 + 1514214159714562294323 #1 + 
    4605782875704094238834 #1^2 + 6617810287079231882646 #1^3 + 
    4605782875704094238834 #1^4 + 1514214159714562294323 #1^5 + 
    189275993388171865138 #1^6 &, (1902663296339164 ArcTan[
        Sin[(22 \[Theta])/25]/(Cos[(22 \[Theta])/25] - #1)] - 
      951331648169582 I Log[1 - 2 Cos[(22 \[Theta])/25] #1 + #1^2] + 
      29652457743237434 ArcTan[Sin[(22 \[Theta])/25]/(
        Cos[(22 \[Theta])/25] - #1)] #1 - 
      14826228871618717 I Log[
        1 - 2 Cos[(22 \[Theta])/25] #1 + #1^2] #1 + 
      77680004276476128 ArcTan[Sin[(22 \[Theta])/25]/(
        Cos[(22 \[Theta])/25] - #1)] #1^2 - 
      38840002138238064 I Log[
        1 - 2 Cos[(22 \[Theta])/25] #1 + #1^2] #1^2 + 
      29652457743237434 ArcTan[Sin[(22 \[Theta])/25]/(
        Cos[(22 \[Theta])/25] - #1)] #1^3 - 
      14826228871618717 I Log[
        1 - 2 Cos[(22 \[Theta])/25] #1 + #1^2] #1^3 + 
      1902663296339164 ArcTan[Sin[(22 \[Theta])/25]/(
        Cos[(22 \[Theta])/25] - #1)] #1^4 - 
      951331648169582 I Log[
        1 - 2 Cos[(22 \[Theta])/
           25] #1 + #1^2] #1^4)/(1514214159714562294323 + 
      9211565751408188477668 #1 + 19853430861237695647938 #1^2 + 
      18423131502816376955336 #1^3 + 7571070798572811471615 #1^4 + 
      1135655960329031190828 #1^5) &]

It's messy for sure but the important part is the RootSum and the expression:

ArcTan[Sin[(22 [Theta])/25]/(Cos[(22 [Theta])/25] - #1)]

which is multivalued with the integration path crossing a default Mathematica branch-cut. First compute the six roots:

poly[x_] = 
  189275993388171865138 + 1514214159714562294323 x + 
   4605782875704094238834 x^2 + 6617810287079231882646 x^3 + 
   4605782875704094238834 x^4 + 1514214159714562294323 x^5 + 
   189275993388171865138 x^6;
theRoots = x /. NSolve[poly[x] == 0, x, WorkingPrecision -> 30]

In order to supply an analytically-continuous version of this function to the RootSum, we integrate six differential equations, one for each root of the RootSum across the integration path:

w0 = ArcTan[Sin[(22 z)/25]/(Cos[(22 z)/25] - theRoots[[4]])] /. 
  z -> myz[0]
wDerivTable = Table[
   w'[z] == (D[
      ArcTan[Sin[(22 z)/25]/(Cos[(22 z)/25] - theRoots[[i]])], z]),
   {i, 1, Length@theRoots}
   ];
maxZ = 50 Pi;
arcSol = NDSolveValue[{#, w[0] == 0}, w, {z, 0, maxZ}, 
    WorkingPrecision -> 25] & /@ wDerivTable

$\texttt{arcSol}$ now are the six analytically-continuous versions of the ArcTan expression in the RootSum.

In the RootSum, replace ArcTan with the six arcSol solutions, replace $\theta$ with z and # with currentRoot:

  newRootSumF[currentArcVal_, theCurrentRoot_, 
  currentZ_] := (6250/
    11 (1902663296339164 ArcTan[Sin[(22 \[Theta])/25]/(
         Cos[(22 \[Theta])/25] - #1)] - 
       951331648169582 I Log[1 - 2 Cos[(22 \[Theta])/25] #1 + #1^2] + 
       29652457743237434 ArcTan[Sin[(22 \[Theta])/25]/(
         Cos[(22 \[Theta])/25] - #1)] #1 - 
       14826228871618717 I Log[
         1 - 2 Cos[(22 \[Theta])/25] #1 + #1^2] #1 + 
       77680004276476128 ArcTan[Sin[(22 \[Theta])/25]/(
         Cos[(22 \[Theta])/25] - #1)] #1^2 - 
       38840002138238064 I Log[
         1 - 2 Cos[(22 \[Theta])/25] #1 + #1^2] #1^2 + 
       29652457743237434 ArcTan[Sin[(22 \[Theta])/25]/(
         Cos[(22 \[Theta])/25] - #1)] #1^3 - 
       14826228871618717 I Log[
         1 - 2 Cos[(22 \[Theta])/25] #1 + #1^2] #1^3 + 
       1902663296339164 ArcTan[Sin[(22 \[Theta])/25]/(
         Cos[(22 \[Theta])/25] - #1)] #1^4 - 
       951331648169582 I Log[
         1 - 2 Cos[(22 \[Theta])/
            25] #1 + #1^2] #1^4)/(1514214159714562294323 + 
       9211565751408188477668 #1 + 19853430861237695647938 #1^2 + 
       18423131502816376955336 #1^3 + 7571070798572811471615 #1^4 + 
       1135655960329031190828 #1^5)) /. {ArcTan[x__] -> 
    currentArcVal, #1 -> theCurrentRoot, \[Theta] -> currentZ}

Create the analytically-continuous antiderivative and plot it over the integration interval:

newAntiD[z_] := Module[{baseSum, endSum},
   baseSum = 0;
   For[i = 1, i <= 6, i++,
    currentArcVal = arcSol[[i]][0];
    theCurrentRoot = theRoots[[i]];
    baseSum += newRootSumF[currentArcVal,theCurrentRoot,z];
    ];
   endSum = 0;
   For[i = 1, i <= 6, i++,
    currentArcVal = arcSol[[i]][z];
    theCurrentRoot = theRoots[[i]];
    endSum += newRootSumF[currentArcVal,theCurrentRoot,z];
    ];
   (endSum - baseSum)
   ];
Plot[Re@newAntiD[z], {z, 0, 125}]

Used Re since will have small imaginary residue from numerical integration.

enter image description here

Improve this answer
$\endgroup$
8
  • 1
    $\begingroup$ Disagree. Both Integrate[Tex, {\[Theta], 0, t}, PrincipalValue -> True ] and Integrate[Tex, {\[Theta], 0, t}, PrincipalValue -> True, Assumptions -> t > 0] return the inputs. Integrate[Tex, {\[Theta], 0, t}, Assumptions -> t >= 0] performs (-2.41733 + 5.17109 I) ArcTan[(0.373211 + 0.00927184 I) Tan[0.44 t]] + 4.92469 ArcTan[ 0.387547 Tan[0.44 t]] + (2.58555 - 1.20866 I) Log[(1. - 3.46945*10^-18 I) - (0.00927184 + 0.373211 I) Tan[ 0.44 t]] - (2.58555 - 1.20866 I) Log[(1. - 3.46945*10^-18 I) + (0.00927184 + 0.373211 I) Tan[0.44 t]]. $\endgroup$
    – user64494
    Dec 6, 2021 at 17:11
  • 1
    $\begingroup$ In the above Tex = 1/(-1509.04 + 1/(-1.64387 - 1.05604 Cos[0.88 \[Theta]]) + 1/(0.000622346 + 0.0000998034 Cos[0.88 \[Theta]]) + 397.918 Cos[0.88 \[Theta]]);. $\endgroup$
    – user64494
    Dec 6, 2021 at 17:12
  • $\begingroup$ Ok thanks. I see that. And if we rationalize all the approx constants via SetPrecision[Tex,Infinity] which I suspect will only rationalize them to 10^-16, then antiderivative is even more complex with a RootSum. $\endgroup$
    – josh
    Dec 6, 2021 at 17:26
  • $\begingroup$ I have the feeling that this integral can be calculated as an analytical function in order to speed up the calculations, but I do not know how to do it. I also tried the above suggestions, but it did not work. If you can help me it would be really useful. $\endgroup$
    – VDF
    Dec 6, 2021 at 19:39
  • $\begingroup$ @VDF. I figured out how to do this by replacing the ArcTan in the RootSum antiderivative with its analytically-continuous version but you'll need to understand how to interpret a RootSum: take the roots in the first expression, substitute each root in the second expression and sum the results, and then how to analytically-continue the ArcTan function in the interval $0\leq \theta \leq 50\pi$ . I'll put something together and update my post above sometime later. $\endgroup$
    – josh
    Dec 7, 2021 at 15:59
1
$\begingroup$

Let me show how to simply add a correction term to continuate the ArcTan branchcuts.

Tex = 1/(-1509.04 + 1/(-1.64387 - 1.05604 Cos[0.88 \[Theta]]) + 
   1/(0.000622346 + 0.0000998034 Cos[0.88 \[Theta]]) + 
   397.918 Cos[0.88 \[Theta]]) /. \[Theta] -> th // Together // 
Expand

Plot[Tex, {th, 0, 20}, PlotRange -> All]

int = Integrate[Tex , th] // Expand

intRe = ComplexExpand[Re@int, TargetFunctions -> {Re, Im}] // 
     Simplify[#, 0 < th < 20] &

Plot[intRe, {th, 0, 20}]

ArcTan has to be corrected with a function of th/2 to get continuation. Have no time to show how i found this funtion, simply show it. Further you have to take signs into accout.

corr = ( 0.88` th/2 - 
  2 ArcTan[ Cot[  0.88` th/2] (-1 + Sqrt[Sec[  0.88` th/2]^2])])

Plot[corr, {th, 0, 20}, PlotRange -> All]

intfin = intRe /. 
ww_ ArcTan[uu_, vv_] :> 
ww (ArcTan[uu, vv] + 
   Sign[vv /. th -> 1.] Sign[uu /. th -> 1.]*corr);

Plot[intfin, {th, 0, 125}]

enter image description here

Improve this answer
$\endgroup$
1
$\begingroup$

I have found the solution in simpler and elegant way, which I post here, which may be useful for some other people, who will met my same problem.

Tex = 1/(-1509.0448660541342` + 
     1/(-1.6438723214688589` - 
      1.0560411247617332` Cos[0.88` \[Theta]]) + 1/(
     0.0006223456294207742` + 
      0.00009980344220763118` Cos[0.88` \[Theta]]) + 
     397.9181707318532` Cos[0.88` \[Theta]]);
pl1 = Plot[NIntegrate[Tex, {\[Theta], 0, x}], {x, 0, 50*Pi}, 
   PlotStyle -> Black];
Tint = Integrate[Tex, \[Theta]];
Pex = Pi/0.88`;
fint[t_] := Tint //. {\[Theta] -> t};
esp[x_] := 
 If[EvenQ[IntegerPart[x]], IntegerPart[(IntegerPart[x] + 2)/2], 
  IntegerPart[(IntegerPart[x] + 1)/2]]
FNint[t_] := esp[(t - Pex)/Pex];
fnint[t_] := If[0 <= t <= Pex, 0, 2*(FNint[t])*fint[Pex]]
Fint[t_] := fint[t] + fnint[t]
pl2 = Plot[Re[Fint[t]], {t, 0, 50*Pi}, PlotPoints -> 100, 
   PlotRange -> All, PlotStyle -> {Red, Dashed}];
Show[pl1, pl2]

enter image description here

Improve this answer
$\endgroup$
0
$\begingroup$

Just to be clearer, when I integrate the function I obtain this

Plot[Sum[Tex, {\[Theta], 0, x, 0.1}]/10, {x, 0, 2*Pi*25}, 
 PlotStyle -> Red, PlotPoints -> 50]

enter image description here

So I would like to obtain this function, but without resorting to integration. I would like to reproduce such a figure through an analytical function. This plot is the mixer of a periodic function and increasing function.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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