1
$\begingroup$

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]

Center

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

Improve this question
$\endgroup$
5
  • 1
    $\begingroup$ Your main issue is that your function is too slow. Also you shouldn't keep calling avgRadius[{x, y}] over-and-over. Cache that result. I'd suggest replacing the EuclideanDistance with its appropriate symbolic definition, calculating symRad = avgRadius[{x, y}] and trying to Compile that result (using Compile[ ... , Evaluate@symRad, ...]). This will hopefully be fast enough to actually compute this integral. $\endgroup$
    – b3m2a1
    Jan 13, 2018 at 18:11
  • 1
    $\begingroup$ You should also point to this old post of yours. $\endgroup$
    – Henrik Schumacher
    Jan 13, 2018 at 18:25
  • $\begingroup$ @b3m2a1 Why is there a singularity? $\endgroup$
    – Arbuja
    Jan 13, 2018 at 18:56
  • 1
    $\begingroup$ @Arbuja never mind it's just a strong peak (which makes more sense). If you just do a Table and ListPlot3D Mathematica clips it. Then there's a different issue in the NIntegrate. All I know is it's unhappy with your function for one reason or another. $\endgroup$
    – b3m2a1
    Jan 13, 2018 at 19:02
  • $\begingroup$ Also your AccuracyGoal is probably way overkill. Try dropping it to 3 or 4. $\endgroup$
    – b3m2a1
    Jan 13, 2018 at 19:04

2 Answers 2

Reset to default
5
$\begingroup$

With the function avgRadiusFunction from this post (hardened version!) and using the mass matrix from the finite element method (and coarsening the mesh considerably, this can be done as follows:

implicitcurve =  ImplicitRegion[ x^2 + y^2 + Sin[4*x] + Sin[4*y] == 4, {{x, -3, 3}, {y, -4, 4}}];
ρ = 0.01;
ρfine = 0.000001;
curvecoarse = DiscretizeRegion[implicitcurve, {{-3, 3}, {-4, 4}},  MaxCellMeasure -> ρ]
curvefine = DiscretizeRegion[implicitcurve, {{-3, 3}, {-4, 4}},  MaxCellMeasure -> ρfine];

enter image description here

As integration domain, I use simply the bounded domain enclosed by curvecoarse.

Needs["NDSolve`FEM`"]
domain = ToElementMesh[curvecoarse, "MeshOrder" -> 2, MaxCellMeasure -> ρ];
domain["Wireframe"]

enter image description here

This computes the mass matrix:

mass = Module[{vd, sd, cdata, x, y, u},
  vd = NDSolve`VariableData[{"DependentVariables","Space"} -> {{u}, {x, y}}];
  sd = NDSolve`SolutionData[{"Space"} -> {domain}];
  cdata = InitializePDECoefficients[vd, sd, "MassCoefficients" -> {{1}}]; 
  DiscretizePDE[cdata, InitializePDEMethodData[vd, sd], sd]["MassMatrix"]
  ];

And this computes the integral $\int_{\operatorname{domain}} \bar r(x,y) \, \operatorname{d} x \, \operatorname{d} y$

{x,y}=Transpose[domain["Coordinates"]];
avgradii = avgRadiusFunction[curvefine][domain["Coordinates"]];
totalavgradius = Total[mass.avgradii]

21.4127

The $x_w$ and $y_w$ can be obtained like this:

xw = (x.mass.avgradii)/totalavgradius
yw = (y.mass.avgradii)/totalavgradius

-0.0689768

-0.0690877

Of course, you are free to decrease with the edge length parameter ρ which controls the accuracy. For example, with ρ = 0.001;ρfine = 0.000001; I get

xw = (x.mass.avgradii)/totalavgradius
yw = (y.mass.avgradii)/totalavgradius

-0.0690308

-0.0690308

Improve this answer
$\endgroup$
2
  • 1
    $\begingroup$ Nice. I clearly need to get comfortable with FEMs. $\endgroup$
    – b3m2a1
    Jan 13, 2018 at 20:19
  • $\begingroup$ How do I use this to solve $x_w$ and $y_w$? Do I take s=x*avgRadiusFunction[curve] and apply Total[mass.s[curvfine][domain["Coordinates"]]] $\endgroup$
    – Arbuja
    Jan 13, 2018 at 20:56
4
$\begingroup$

I mentioned this in the comments, but we'll tackle this via compilation and a more manageable AccuracyGoal:

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 -> 4];
q = MeshCoordinates[curve];
edges = First /@ MeshCells[curve, 1];
avgRadius[p_] :=

  1/(2 π) Abs[
    Sum[Module[{q1, q2, d1, d2, r, dθ}, q1 = q[[First@e]];
      q2 = q[[Last@e]];
      r = Norm[p - (q1 + q2)/2];(*midpoint approximation*)

      d1 = q1 - p;
      d2 = q2 - p;
      dθ =
       Arg[(d1[[1]] + I*d1[[2]])/(d2[[1]] + I*d2[[2]])];
      r dθ], {e, edges}]];
symbRad = avgRadius[{x, y}];
compSym =
  Compile[{{x, _Real}, {y, _Real}},
   Evaluate@symbRad,
   RuntimeOptions -> {"EvaluateSymbolically" -> False},
   CompilationTarget -> "C"
   ];

Then we can see how long it takes on a reasonably small grid:

AbsoluteTiming[
  heights =
    Table[
     compSym[x, y],
     {x, -10, 10, .1},
     {y, -10, 10, .1}
     ];
  ][[1]]

3.49649

Which is still mad slow. So the NIntegrate will probably take a looong time.

And so first we check to see if we messed up and introduced a "MainEvaluate":

<< CompiledFunctionTools`
pstr = CompilePrint[compSym];
StringContainsQ[pstr, "MainEvaluate"]

False

And it turns out we didn't.

Now plotting those heights:

ListPlot3D[heights, PlotRange -> All]

asdasd

This seems like the kind of NIntegrate should be able to handle. So:

NIntegrate[
 compSym[x, y],
 {x, -10, 10},
 {y, -10, 10}
 ]

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.

142.029

And so there's something in your domain making NIntegrate unhappy.

The other integrals can also be handled via compilation, e.g.:

compXSym =
  Compile[{{x, _Real}, {y, _Real}},
   Evaluate[x * symbRad],
   RuntimeOptions -> {"EvaluateSymbolically" -> False},
   CompilationTarget -> "C"
   ];

compYSym =
  Compile[{{x, _Real}, {y, _Real}},
   Evaluate[y * symbRad],
   RuntimeOptions -> {"EvaluateSymbolically" -> False},
   CompilationTarget -> "C"
   ];

But fundamentally NIntegrate is unhappy with something in your problem and it's worth trying to figure out what that is.

Improve this answer
$\endgroup$
1
  • $\begingroup$ Maybe the code from this answer could help? $\endgroup$
    – Arbuja
    Jan 13, 2018 at 19:44

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.