2
$\begingroup$

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] 

enter image description here

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:

enter image description here

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?

$\endgroup$
6
  • 1
    $\begingroup$ Welcome to Mathematica.SE! 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2) Take the tour and check the faqs! 3) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! $\endgroup$
    – user9660
    Mar 22, 2016 at 8:50
  • 2
    $\begingroup$ I always try to plot the function I'm integrating to get an idea of why NIntegrate is having trouble, have you done this? You see the discontinuity at τ = 0.441549, is this the singularity you want to find? Also, do you need the answer to such high precision or are you just trying to get it to work without giving an error? If the latter, (NIntegrate[-First[curve02] D[Last@curve02, \[Tau]] + Last[curve02] D[First@curve02, \[Tau]], {\[Tau], 0, \[Pi]/6}, MaxRecursion -> 100, Method -> "LocalAdaptive"]) works very quickly and without error. $\endgroup$
    – Jason B.
    Mar 22, 2016 at 9:43
  • $\begingroup$ @user6043040 - see below, I've included 2 methods for the inverse problem, and one method for finding the point of discontinuity. $\endgroup$
    – Jason B.
    Mar 22, 2016 at 12:01
  • $\begingroup$ I used τp = (y - x) /. FindRoot[(curve02 /. τ -> x) == (curve02 /. τ -> y), {x, Pi/20}, {y, 2 Pi}]; to estimate the period of $\tau$. This value is also used to display the second figure in the post. $\endgroup$ Mar 22, 2016 at 12:34
  • $\begingroup$ This is related. $\endgroup$ Mar 22, 2016 at 22:57

2 Answers 2

4
$\begingroup$

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}]

enter image description here

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}} *)
$\endgroup$
3
  • $\begingroup$ Thank you this exactly solved all my problems on it. Why Evaluated -> False option is necessary for this case? I tried it without this option, but it seems no iteration occurs, and the return is only the initial $x$ I used. $\endgroup$ Mar 22, 2016 at 12:03
  • 1
    $\begingroup$ I googled "NIntegrate in FindRoot" and found this page, and used the solution presented there $\endgroup$
    – Jason B.
    Mar 22, 2016 at 12:05
  • $\begingroup$ That's it. The error message does indicate ...has evaluated to non-numerical values for all sampling points . I just thought it was weird, now clear! thank you very much! $\endgroup$ Mar 22, 2016 at 12:21
2
$\begingroup$

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] &

discretized

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 *)
$\endgroup$
1
  • $\begingroup$ Thank you MacroB. As a matter of fact, for calculation of the area, I also used a polygon to approximate it, then shoelace method(discrete version of Green's theorem,en.wikipedia.org/wiki/Shoelaces) can be applied. But such approximation may not be appropriate for the inverted problem. $\endgroup$ Mar 23, 2016 at 2:14

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.