How to find the integrand singularity points when having NIntegrate::slwcon:

Question

In order to calculate the closed area of the curve below defined by parametric equation curve02,

ClearAll["Global`*"];
R = 4800/100;
z1 = 6;
r = R/z1;
z2 = z1 - 1;
e = 705/100;
f = r/e;
re = 126/10;
θ = ArcTan[Sin[z1 τ]/(f + Cos[z1 τ])] - τ;
φ = ArcSin[f Sin[θ + τ]] - θ;
ψ = z1/(z1 - 1) φ;

(* definition of parametric equations *)
curve01 = {(R - r) Sin[τ] + e Sin[z2 τ] - re Sin[θ], (R - r) Cos[τ] - e 
Cos[z2 τ] + re Cos[θ]}
// Simplify;
(*show stator curve*)
(* ParametricPlot[curve01, {τ, 0, 2 π}, Exclusions -> None, 
 MaxRecursion -> 15, PlotPoints -> 500, PlotStyle -> Red] *)


(*rotator curve*)
curve02 = {curve01[[1]] Cos[φ - ψ] - curve01[[2]] Sin[φ - ψ] - e Sin[ψ], 
    curve01[[1]] Sin[φ - ψ] + curve01[[2]] Cos[φ - ψ] - e Cos[ψ]} 
 //Simplify;

base = ParametricPlot[curve02, {τ, 0, 2 π}, 
  Exclusions -> None, MaxRecursion -> 15, PlotPoints -> 500, 
  PlotStyle -> Blue]; Show[base] 

I use NIntegrate per Green's theorem as below since the period of $\tau$ can be easily obtained as $5\pi/3$ considering the symmetry of the curve:

5 (NIntegrate[-First[curve02] D[Last@curve02, τ] + 
    Last[curve02] D[First@curve02, τ], {τ, 0, π/6}, 
   WorkingPrecision -> 50, PrecisionGoal -> 12, MaxRecursion -> 100])

Though the code results:

6385.1150304636387606396132810257250079650388434246

There is always warning message:

NIntegrate::slwcon: Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small.

Though I understand from the documents that if I know the singularity points of the integrand and specify them such warning message would be removed. But How can I easily find such singularity points especially for such a complicate case?

Additionally, I also want to solve the inverse problem of the above issue, i.e., given the area value, say, 6557.23, find the corresponding $r$ by using Newton's method in FindRoot.

I used the following function defined in Mathematica:

curveArea[r_]:=Module[{z1,z2,e,f,re,R,θ,φ,ψ,τ,curve01,curve02},
z1=6;
R=r z1;
z2=z1-1;
e=705/100;
f=r/e;
re=126/10;
θ=ArcTan[Sin[z1 τ]/(f+Cos[z1 τ])]-τ;
φ=ArcSin[f Sin[θ+τ]]-θ;
ψ=z1/(z1-1) φ;
curve01={(R-r) Sin[τ]+e Sin[z2 τ]-re Sin[θ],(R-r) Cos[τ]-e Cos[z2 τ]+re Cos[θ]};curve02={curve01[[1]] Cos[φ-ψ]-curve01[[2]] Sin[φ-ψ]-e Sin[ψ],curve01[[1]] Sin[φ-ψ]+curve01[[2]] Cos[φ-ψ]-e Cos[ψ]};
Quiet@(5 (NIntegrate[-First[curve02] D[Last@curve02,τ]+Last[curve02] D[First@curve02,τ],{τ,0,τp/20,τp/10},PrecisionGoal->8,MaxRecursion->100]))
]

but FindRoot just gives more error or warning messages and does not work:

FindRoot[curveArea[x] == 6557.23, {x, 7.5}]

How can I solve such a problem elegantly?

Answer 1

Firstly, you should try the different methods available for NIntegrate as some of them may be able to help with a given problem, for example,

5 (NIntegrate[-First[curve02] D[Last@curve02, τ] + 
    Last[curve02] D[First@curve02, τ], {τ, 0, π/6}, 
   MaxRecursion -> 100, Method -> "LocalAdaptive"])
(* 6385.116632516242 *)

works without giving an error. But if you try to increase the PrecisionGoal it will slow down considerably.

Your more general question of "How to find the integrand singularity points" is interesting here. I always like to plot the function being integrated to see what is going on,

func[τ_] = -First[curve02] D[Last@curve02, τ] + 
   Last[curve02] D[First@curve02, τ];
Plot[func[t], {t, 0, π/6}]

So you have what looks like a jump discontinuity around t=0.44, but if you want to feed this in to NIntegrate you need that value to higher precision. This is kind of a hack method, but it basically iteratively finds the discontinuity by building a list of the function, then applying Differences to that list to find where the jump is.

Module[{min = 0, max = π/6., vals, data},
 Table[
  vals = Range[min, max, 10.^-n];
  data = f1[#] & /@ vals;
  discontinuity = 
   Part[vals, Last@Ordering[Abs@Differences@data] + 1];
  {min, max} = discontinuity + 10.^-n {-1, 1};
  discontinuity
  , {n, 3, 8}]
 ]

(* {0.442`, 0.4416`, 0.44155`, 0.44155`, 0.44154950000000004`, 0.44154945`} *)

Now you can do your original integral,

5 (NIntegrate[-First[curve02] D[Last@curve02, τ] + 
    Last[curve02] D[First@curve02, τ], {τ, 0, discontinuity, π/6},
    WorkingPrecision -> 50, PrecisionGoal -> 12, MaxRecursion -> 100])

(* 6385.1150254108130393191308125092576916637822038014 *)

Your second question, about how to invert the procedure to find the value of r which gives the corresponding area, is also interesting. Your code contained calls to τp which is undefined, so I'll keep the integration limits the same as above. I will also use the LocalAdaptive method since it works without error,

curveArea[r_] := 
 Module[{z1, z2, e, f, re, R, θ, φ, ψ, τ, 
   curve01, curve02}, z1 = 6;
  R = r z1;
  z2 = z1 - 1;
  e = 705/100;
  f = r/e;
  re = 126/10;
  θ = ArcTan[Sin[z1 τ]/(f + Cos[z1 τ])] - τ;
  φ = ArcSin[f Sin[θ + τ]] - θ;
  ψ = z1/(z1 - 1) φ;
  curve01 = {(R - r) Sin[τ] + e Sin[z2 τ] - 
     re Sin[θ], (R - r) Cos[τ] - e Cos[z2 τ] + 
     re Cos[θ]}; 
  curve02 = {curve01[[1]] Cos[φ - ψ] - 
     curve01[[2]] Sin[φ - ψ] - e Sin[ψ], 
    curve01[[1]] Sin[φ - ψ] + 
     curve01[[2]] Cos[φ - ψ] - e Cos[ψ]};
  Quiet@(5 (NIntegrate[-First[curve02] D[Last@curve02, τ] + 
        Last[curve02] D[First@curve02, τ], {τ, 0, π/6}, 
       MaxRecursion -> 100, Method -> "LocalAdaptive"]))]

Using FindRoot works in this case, but is slow.

FindRoot[curveArea[x] == 6557.23, {x, 8.0}, 
  Evaluated -> False] // AbsoluteTiming

During evaluation of FindRoot::lstol: The line search decreased the step size to within tolerance specified by AccuracyGoal and PrecisionGoal but was unable to find a sufficient decrease in the merit function. You may need more than MachinePrecision digits of working precision to meet these tolerances. >>

(* {11.2733, {x -> 8.12944}} *)

If you want to repeatedly apply this, getting r for many different areas, you can build an interpolating function and apply FindRoot to it,

interp = Interpolation@Table[{r, curveArea[r]}, {r, 7, 10, .05}];

FindRoot[interp[x] == 6557.23, {x, 8.}, 
  Evaluated -> False] // AbsoluteTiming
(* {0.000659, {x -> 8.12944}} *)

Answer 2

Here is an approach to calculating the area by the curve using geometric computation and region functionality, rather than explicitly through Green's theorem:

base = ParametricPlot[                  
    curve02, {τ, 0, 5 Pi/3},                    (* note: 5 Pi/3 is a more natural endpoint *)
    Exclusions -> None, MaxRecursion -> 15,
    PlotStyle -> Blue, Axes -> False
  ];

region = Cases[base, Line[a__], Infinity] //
     First //
     FilledCurve //
     Graphics //
     DiscretizeGraphics[#, MaxCellMeasure -> 0.5] &

We can then use NIntegrate directly to calculate the area, or alternatively Area or RegionMeasure:

NIntegrate[1, {x, y} ∈ region]   (* Out: 6384.96 *)
Area[region]                     (* Out: 6384.96 *)
RegionMeasure[region]            (* Out: 6384.96 *)