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?
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$τ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$