0
$\begingroup$

The following gives compile warnings and returns unevaluated, any idea why?

obj2[a_] := 
 NIntegrate[(
   E^(1/2 (-x[1]^2 - x[2]^2))
     Log[Abs[1 - a (x[1]^2 + x[2]^2)]] x[
     1]^2)/(\[Pi] (x[1]^2 + x[2]^2)), {x[
     1], -\[Infinity], \[Infinity]}, {x[
     2], -\[Infinity], \[Infinity]}] + Log[2]/2;
obj2[.915]   (* -0.00319086 + 5.42344*10^-18 I *)

FindRoot[obj2[a], {a, 0.92}]   (* compile errors + returns unevaluated *)

Also, I'm curious why the result of NIntegrate is complex, the only suspect part is Log[Abs[x]], but it should be real for all x.

(it's the objective from here)

$\endgroup$
1
  • $\begingroup$ Using obj2[a_?NumericQ] and forcing WorkingPrecision -> 15 with NIntegrate works fine for me. May be due to Log and Abs and the singularity at Abs[1 - 183/200 (x[1]^2 + x[2]^2)] == 0 there is loss of workingPrecision below machinePrecision. $\endgroup$
    – Akku14
    Aug 23, 2022 at 18:08

2 Answers 2

2
$\begingroup$

The two problems seem to be related, or at least can be cured together:

obj2[a_?NumericQ] := 
  Chop[NIntegrate[(E^(1/2 (-x[1]^2 - x[2]^2)) Log[
         Abs[1 - a (x[1]^2 + x[2]^2)]] x[1]^2)/(\[Pi] (x[1]^2 + 
         x[2]^2)), {x[1], -\[Infinity], \[Infinity]}, {x[
       2], -\[Infinity], \[Infinity]}] + Log[2]/2];
obj2[.915]  (* -0.00319086 *)
FindRoot[obj2[a], {a, 0.92}]  (* {a -> 0.918611} *)

The two modifications are:

  • Using a test to sanitize the inputs of obj2. Using the pattern a_?NumericQ inside a function definition means that the function will only be evaluated if NumericQ evaluates to True on an argument, and will be return unevaluated if it does not. In particular, this prevents FindRoot from trying to take symbolic derivatives of obj2 in its internal algorithms, which it will do by default.

  • Chop removes "small" complex parts from numerical outputs. Like you, I'm unsure why NIntegrate is returning a complex part for this integrand at all; but it apparently causes FindRoot to invoke some kind of algorithm that does not converge. Making sure that obj2 only ever returns a real output makes FindRoot much happier.

$\endgroup$
1
  • $\begingroup$ OP: Please don't accept this answer immediately, as the Mathematica internals wizards on this forum might be able to provide a better insight into why all of this happened in the first place. $\endgroup$ Aug 23, 2022 at 16:24
1
$\begingroup$

Use the settings by @Akku14, we can find two roots in 0<a<1.

Clear[obj2]; 
obj2[a_?NumericQ] := 
 NIntegrate[(E^(1/2 (-x[1]^2 - x[2]^2)) Log[
       Abs[1 - a (x[1]^2 + x[2]^2)]] x[1]^2)/(π (x[1]^2 + 
        x[2]^2)), {x[1], -∞, ∞}, {x[
     2], -∞, ∞}, WorkingPrecision -> 15] + 
  Log[2]/2;
FindRoot[obj2[a], {a, 0.2}]

FindRoot[obj2[a], {a, 0.8}]

ReImPlot[obj2[a], {a, 0, 1}]

{a -> 0.121312}

{a -> 0.918611}

enter image description here

$\endgroup$

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.