1

I am trying to evaluate an integral that has two symbolic parameters, $x$ and $y$ in it

$$-\int_0^1\frac{\log (u)}{\sqrt{1-u^2}} \exp \left(\frac{2 u}{u+1}xy-\frac{u^2 (x-y)^2}{1-u^2}\right) \; du$$

Of course an analytic formula would be desirable, but I can't find it neither by hand nor with Mathematica. Therefore, I've tried the following numerical evaluation

NIntegrate[- Log[u]/Sqrt[1 - u^2]Exp[(2 u)/(1 + u) xy - u^2/(1 - u^2) (x - y)^2], {u, 0, 1}]

This idea makes sense to me since NIntegrate can find the result for specific $x$ and $y$ values, so there should be a formula that I should find. I'm getting the error

NIntegrate::inumr: The integrand -((E^((2 u xy)/(1+u)-(u^2 (x-y)^2)/(1-Power[<<2>>])) Log[u])/Sqrt[1-u^2]) has evaluated to non-numerical values for all sampling points in the region with boundaries {{0,1}}.

Then, I tried to look up what might cause the issue, and I ended up at this 9 year old question, describing exactly the same problem. Now, it looks like the problem was solved there with comments, but as it often happens with very old questions, the link is broken, and the solution is not explained in comments.

Could someone point me to the right direction on why NIntegrate is not able to handle symbolic parameters, and also on how to solve this problem using a more recent version of Mathematica (maybe the solution had changed in the last 9 years)?

  • Is parameter a a typo? – Ulrich Neumann Feb 16 '22 at 09:35
  • Yes, indeed, Ive realized I had multiple typos (had too many failed attempts in my Notebook and somehow managed to copy wrong formulas) –  Feb 16 '22 at 09:48
  • 1
    This is an unusual question. You are asking why NIntegrate cannot handle symbolic parameters, despite this being documented and, to be blunt, fairly obvious. What method might it use to give a symbolic result that Integrate cannot provide? – Daniel Lichtblau Feb 16 '22 at 13:44
  • What is the domain for {x, y}? – Michael E2 Feb 16 '22 at 15:02
  • The link from the 9 year old question is to the article "How do I use ?NumericQ to affect order of evaluation?" and may be found here, together with links to Q&A that give the same answer as @Ulrich. – Michael E2 Feb 16 '22 at 15:23

2 Answers2

2

You might calculate as follows. Let us define the function:

J[x_, y_] :=NIntegrate[-Log[u]/Sqrt[1 - u^2] Exp[(2 u)/(1 + u) x*y - 
      u^2*(x - y)^2/(1 - u^2)], {u, 0, 1}];

Then let us calculate the list of integrals in various points x and y, within the intervals 0<x,y<1:

lst = Flatten[Table[{x, y, J[x, y]}, {x, 0, 1, 0.1}, {y, 0, 1, 0.1}],1]

The result is rather long, and I do not show it here, but its evaluation is straightforward and fast.

Now let us plot it:

ListPlot3D[lst, PlotRange -> All, 
 AxesLabel -> {Style["x", Italic, 16], Style["y", Italic, 16], 
   Style["J", Italic, 16]}]

enter image description here

If you need an analytic formula for J=J(x,y) you can look for a suitable fit for the list lst.

This approach has a born-in drawback: one can only do the calculation within a finite domain in the (x,y) plane.

Edit: To address your question. Let us fit the result:

model = a Exp[b*x*y + c (x - y)^2];
ff = FindFit[lst, model, {a, b, c}, {x, y}]

(* {a -> 1.07657, b -> 0.414825, c -> -0.210946}  *)

And now let us inspect the fitting visually:

Show[{
  ListPointPlot3D[lst, PlotRange -> All, 
   AxesLabel -> {Style["x", Italic, 16], Style["y", Italic, 16], 
     Style["J", Italic, 16]}, PlotStyle -> {Red, PointSize[0.01]}],
  Plot3D[model /. ff, {x, 0, 1}, {y, 0, 1}, 
   PlotStyle -> {Opacity[0.5], Blue}]
  }]

yielding this:

enter image description here

My feeling is that this fit is reasonably good.

A minor comment: The letter C is reserved in Mma. Its use in fitting is an error. In general, it is recommended to avoid capital letters and to always begin any custom variable name with lowercase ones.

Have fun!

Alexei Boulbitch
  • 39,397
  • 2
  • 47
  • 96
  • I like the idea of fitting the table. Unfortunately I am new to this, and I am running into the following issue:

    lst = Flatten[Table[{x, y, f[x, y]}, {x, -2, 2, 0.1}, {y, -2, 2, 0.1}]];

    model[A_, B_, C_, x_, y_] := A Exp[B x y - C (x - y)^2]

    FindFit[lst, model, {A, B, C}, {x, y}]

    returns the error: FindFit::fitc: Number of coordinates (1) is not equal to the number of variables (2).

    Do you know what I am missing?

     –  Feb 16 '22 at 10:00
  • You could also interpolate the data. The interpolation can be used in subsequent calculations as if it was a function. – Hugh Feb 16 '22 at 10:50
  • @La Belle Croissant Please see the edit. – Alexei Boulbitch Feb 16 '22 at 13:40
2

If I take the first expression ("without parameter" a), define a function depending only on numerical arguments x,y

int[x_?NumericQ, y_?NumericQ] := NIntegrate[-Log[u]/Sqrt[1 - u^2] Exp[(2 u)/(1 + u) x y - u^2 (x - y)^2/(1 - u^2)], {u, 0, 1}]

int[1,1] (1.66861)
Ulrich Neumann
  • 53,729
  • 2
  • 23
  • 55