3
$\begingroup$

I have a function

f[pt_] := Log[1/(1 - Abs[x])] + Log[1/(1 - Abs[y])] /. {x -> pt[[1]], y -> pt[[2]]}

That looks like this:

enter image description here

I want to integrate it over a rotated rectangle:

Plot[NIntegrate[f[{x, y}], {x, y} \[Element] 
     TransformedRegion[Rectangle[{-1/2, -1/2}, {1/2, 1/2}], 
         RotationTransform[theta]]], {theta, 0, 2 Pi}]

But it's very slow, and I keep getting errors like:

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.

I don't think this function should be that difficult to integrate.

I tried to rotate it myself like this:

r[t_] := {{Cos[t], -Sin[t]}, {Sin[t], Cos[t]}}
Plot[NIntegrate[f[r[theta].{x, y}], {x, -1/2, 1/2}, {y, -1/2, 1/2}], {theta, 0, 2 Pi}]

But it is still very slow. As in I had to kill it after multiple minutes of no results. I'm using Mathematica 12.1.1.0. Is there anything I can do to speed this up?


Bonus problem: I also want to find the maximum of the function. So I take

FindMaximum[{NIntegrate[f[r[theta].{x, y}], {x, y} \[Element] Rectangle[{-1/2, -1/2}, {1/2, 1/2}]], 0 <= theta <= 2 Pi}, {theta, .1}]

However, this gives me a bunch of errors like

NIntegrate::inumr: The integrand g[x Cos[theta]-y Sin[theta],
y Cos[theta]+x Sin[theta]] has evaluated to non-numerical values
for all sampling points in the region with boundaries
{{-(1/2),0.},{-(1/2),1/2}}.

I tried to change the definition of r to

r[t_?NumericQ] := {{Cos[t], -Sin[t]}, {Sin[t], Cos[t]}}

But it doesn't seem to improve things. It only changes the error to

NIntegrate::inumr: The integrand g[r[theta],{x,y}] has evaluated to
non-numerical values for all sampling points in the region with
boundaries {{-(1/2),1/2},{-(1/2),1/2}}.

Note that the function is definitely defined at this location, and for all valid points in the region for any rotation.

Can anyone help me understand what I'm doing wrong?

$\endgroup$
2
  • $\begingroup$ Work fine in v13.1 $\endgroup$
    – cvgmt
    Oct 25, 2022 at 4:48
  • $\begingroup$ @cvgmt I just tried on Mathematica Online (which I presume is v13) and the plots are still very slow, while the Maximization still gives "non-numerical values" errors. (Though it doesn't find a result.) $\endgroup$ Oct 25, 2022 at 4:56

2 Answers 2

5
$\begingroup$

We can also fix the integral region but transform the function by RotationMatrix.

f[{x_, y_}] = Log[1/(1 - Abs[x])] + Log[1/(1 - Abs[y])];
Plot[NIntegrate[
  Evaluate[f[{x, y} . RotationMatrix[θ]]], {x, -1/2, 
   1/2}, {y, -1/2, 1/2}], {θ, 0, 2 π}]
int[θ_?NumericQ] := 
 NIntegrate[
  f[{x, y} . RotationMatrix[θ]] // Evaluate, {x, -1/2, 
   1/2}, {y, -1/2, 1/2}]
FindMaximum[{int[θ], 0 <= θ <= 2 π}, {θ, .1}]

enter image description here

{0.613705, {θ -> 0.00170899}}

$\endgroup$
2
  • $\begingroup$ Is the trick here the use of // Evaluate or the ?NumericQ in int? It's still not intuitive for me when to use what $\endgroup$ Oct 25, 2022 at 16:25
  • $\begingroup$ @ThomasAhle The ?NuericQ skill for FindMaximum, and the Evaluate is not necessary, the main ideas is using the integral transformation, avoid transform the domain. $\endgroup$
    – cvgmt
    Oct 25, 2022 at 23:19
4
$\begingroup$
$Version

(* "13.1.0 for Mac OS X x86 (64-bit) (June 16, 2022)" *)

Clear["Global`*"]

f[pt_] := Total[Log[1/(1 - Abs[#])] & /@ pt]

rgn[theta_] = 
  TransformedRegion[Rectangle[{-1/2, -1/2}, {1/2, 1/2}], 
    RotationTransform[theta]] // Simplify;

int[theta_?NumericQ] := 
 NIntegrate[f[{x, y}], {x, y} ∈ rgn[theta]]

Plot[int[theta], {theta, 0, 2 Pi}] // AbsoluteTiming

enter image description here

({max, arg} = FindMaximum[int[theta], {theta, 1.6}]) // AbsoluteTiming

(* {1.1356, {0.613706, {theta -> 1.57082}}} *)
$\endgroup$
3
  • $\begingroup$ What is the whole trick with Total and &/@? $\endgroup$ Oct 25, 2022 at 20:24
  • 1
    $\begingroup$ Each of the terms in your function is just the pure function Log[1/(1 - Abs[#])] & evaluated with the parts of pt. So the pure function is mapped ( /@ ) onto the elements of pt. Then the Total adds terms to produce the result. See the documentation for Function, Map, Total $\endgroup$
    – Bob Hanlon
    Oct 25, 2022 at 20:39
  • $\begingroup$ That makes sense. Thanks. Do you think there's anything I can do to speed up the function? As you've shown, currently plotting takes nearly a minute. Is there maybe a way to use NDSolve instead of NumericalIntegrate? $\endgroup$ Oct 25, 2022 at 21:19

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.