In a previous post, I wanted to compute the weighted average, with weight $\overline{r}(x,y)$, of the whole cartesian plane. This was not possible since the integrals did not converge.
Instead I'm computing the weighted average of points inside Region $\mathcal{C}$.
The definition of $\overline{r}(x,y)$ is the average ray function
Suppose we have curve $\mathcal{C}$ and a ray, from $\mathbf{p}\in\mathcal{C}$, rotating full circle at $\Delta \theta$ increments while intersecting $\mathcal{C}$. We will denote the intersections as $\left\{\mathbf{q}_i\right\}_{i=1}^{n}$ (for $n$ intersections) and the length of the ray as $r_i = \|{\mathbf{q}_i}\mathbf - \mathbf p\|$.
Then as $\Delta \theta\to0$, the average radius function is $\overline{r}(\mathbf{p})=\left(\sum\limits_{i=1}^{n}r_i\right)/{n}$
For star shaped curves, one can use an alternate and simpler definition
$$\overline{r}(\mathbf{p})=\frac1{2\pi}\oint_{\mathbf {q}_1\in\mathcal C}\|\mathbf {q}_1-\mathbf p\|\,\mathrm d\theta.$$
Thanks to @Rahul, in this post, I have a code that can implement $\overline{r}(x,y)$. For example, if $\mathcal{C}$ is defined by $x^2+y^2+\sin(4x)+\sin(4y)=4$
curve = DiscretizeRegion[
ImplicitRegion[
x^2+y^2+Sin[4*x]+Sin[4*y] == 4, {{x, -3, 3}, {y, -4, 4}}], {{-3, 3},
{-4, 4}}, AccuracyGoal -> 8]
q = MeshCoordinates[curve];
edges = First /@ MeshCells[curve, 1];
signedAngle[a_, b_] := Arg[(Complex @@ a)/(Complex @@ b)]
avgRadius[p_] :=
1/(2 \[Pi]) Abs[Sum[Module[{q1, q2, r, d\[Theta]}, q1 = q[[First@e]];
q2 = q[[Last@e]];
r = EuclideanDistance[p, (q1 + q2)/2];(*midpoint approximation*)
d\[Theta] = signedAngle[q1 - p, q2 - p];
r d\[Theta]], {e, edges}]]
avgRadius[{x, y}]
Edit: A faster code can be implemented using @HenrySchumacer's answer in the this post.
avgRadiusFunction[curve_MeshRegion] :=
Module[{q, e1, e2}, q = MeshCoordinates[curve];
{e1, e2} = Transpose[Developer`ToPackedArray[MeshCells[curve, 1][[All, 1]]]];
With[{q11 = q[[e1, 1]], q12 = q[[e1, 2]], q21 = q[[e2, 1]], q22 = q[[e2, 2]]},
Compile[{{p, _Real, 1}},
Module[{u11, u12, u21, u22, r, dθ, r2squared, r2, p1, p2},
p1 = Part[p, 1];
p2 = Part[p, 2];
u11 = q11 - p1;
u12 = q12 - p2;
u21 = q21 - p1;
u22 = q22 - p2;
r2squared = (u21^2 + u22^2);
dθ = ArcTan[(u11 u21 + u12 u22)/r2squared, (u12 u21 - u11 u22)/r2squared];
r2 = Sqrt[r2squared];
(RotateRight[r2] + r2).dθ/(4. Pi)
],
RuntimeAttributes -> {Listable},
Parallelization -> True]
]
];
curve = DiscretizeRegion[ ImplicitRegion[x^2+y^2+Sin[4*x]+Sin[4*y] == 4, {{x, -3, 3}, {y, -4, 4}}], {{-3, 3},{-4, 4}}, AccuracyGoal -> 8]
f = avgRadiusFunction[curve];
We can plot the contour of $\overline{r}(x,y)$ on $\mathcal{C}$.
k = ContourPlot[f[{x, y}], {x, -3, 3}, {y, -3, 3},
Exclusions -> None, Contours -> 100, PlotLegends -> Automatic]
Hence the weighted average in $\mathcal{C}$ is
$$\begin{array}{cc} {x_w=\frac{\int\int_{\mathcal{C}}x\overline{r}(x,y) \ dA }{\int\int_{\mathcal{C}}\overline{r}(x,y) \ dA}} & { y_w=\frac{\int\int_{\mathcal{C}} y\overline{r}(x,y) \ dA}{\int\int_{\mathcal{C}}\overline{r}(x,y) \ dA} }\end{array}$$
However, even for the simplest double integrals, I'm unable to approximate using NIntegrate
In less than an hour.
x_w=(NIntegrate[x*avgRadius[{x,y}],{x,-10,10},
{y,-10,10}])/(NIntegrate[avgRadius[{x,y}],{x,-10,10},{y,-10,10}])
y_w=(NIntegrate[y*avgRadius[{x,y}],{x,-10,10},{y,-10,10}])/(NIntegrate[avgRadius[{x,y}],{x,-10,10},{y,-10,10}])
How do I do accurately and efficiently approximate the weighted average of the $x$ and $y$-coordinates over curve $\mathcal{C}$ with weights $\overline{r}(x,y)$?
"Accurately" means the error should be less than .00001. "Efficently" means computing the integral under $15$ minutes? I am unable to do so using NIntegrate
avgRadius[{x, y}]
over-and-over. Cache that result. I'd suggest replacing theEuclideanDistance
with its appropriate symbolic definition, calculatingsymRad = avgRadius[{x, y}]
and trying toCompile
that result (usingCompile[ ... , Evaluate@symRad, ...]
). This will hopefully be fast enough to actually compute this integral. $\endgroup$Table
andListPlot3D
Mathematica clips it. Then there's a different issue in theNIntegrate
. All I know is it's unhappy with your function for one reason or another. $\endgroup$AccuracyGoal
is probably way overkill. Try dropping it to 3 or 4. $\endgroup$