4
$\begingroup$

I would like to approximate the integral $$ \int_0^\infty dy\,\frac{1}{\sqrt{2\pi y\sigma^2}}\exp\left(-\frac{(x-y)^2}{2y\sigma^2}\right)f(y), $$ as a series expansion in the limit $\sigma\rightarrow 0$, where $f(y)$ is an arbitrary function. I can make a change of variables $$ z=\frac{x-y}{\sigma\sqrt{y}} $$ and then expand the integrand as a power series in $\sigma$ but that is a lot of manual manipulation (particularly because the change of variables is not injective).

Is there an easier way to approximate the integral in Mathematica?

What I have so far

Define the integrand substitution rule

ruleZ=z -> (x - y)/(Sqrt[y] σ);

Obtain expressions for y in terms of z

slnY = Solve[z == (z /. ruleZ), y]
Out: {{y -> 1/2 (2 x + z^2 σ^2 - Sqrt[4 x z^2 σ^2 + z^4 σ^4])}, 
      {y -> 1/2 (2 x + z^2 σ^2 + Sqrt[4 x z^2 σ^2 + z^4 σ^4])}}

Define the integrand

ρ = σ Sqrt[y];
integrand = 1/Sqrt[2 Pi ρ^2] Exp[-(x - y)^2/(2 ρ^2)] f[y];

Change variables to $z$

integrandZ = FullSimplify[(integrand/Abs[D[z /. ruleZ, y]]) /. slnY];
Out: too nasty to reproduce

Expand the integrand as a series

approximateIntegrand=Series[integrandZ, {σ, 0, 2}];
Out: too nasty to reproduce

Integrate and obtain coefficients of $\sigma$

result = CoefficientList[Integrate[approximateIntegrand, {z, -∞, ∞}, σ]
Out: {{f[x], -(f[x]/(Sqrt[2 π] Sqrt[x])) - Sqrt[2/π] Sqrt[x] Derivative[1][f][x], 
       Derivative[1][f][x] + 1/2 x (f^′′)[x]}, 
      {f[x], f[x]/(Sqrt[2 π] Sqrt[x]) + Sqrt[2/π] Sqrt[x] Derivative[1][f][x], 
       Derivative[1][f][x] + 1/2 x (f^′′)[x]}}

In the limit $\sigma\rightarrow 0$ both solutions converge to $f(x)$ as expected. The first order terms are equal but opposite and the second order terms are equal.

What is still causing me trouble

I would expect the first order term to vanish. The full solution should be a sum over all branches of the change of variables (cf. http://www.math.utah.edu/~levin/M5080/transforms.pdf page 3) but I obtain twice the result I expect in said case.

Improve this question
$\endgroup$
5
  • $\begingroup$ Why not evaluate the integral directly? $\endgroup$
    – Julian
    Commented Feb 26, 2015 at 12:02
  • $\begingroup$ Because $f(y)$ is an arbitrary function. $\endgroup$
    – Till Hoffmann
    Commented Feb 26, 2015 at 12:03
  • $\begingroup$ (at) Julian: See my solution, which assumes that f[x] can be expanded into a power series. $\endgroup$
    – Dr. Wolfgang Hintze
    Commented Feb 26, 2015 at 12:45
  • $\begingroup$ Because you have replaced σ by ρ in the integrand, I think you should expand in ρ, not σ, which now is just a coefficient in the definition of z that has no connection to ρ. More generally, it seems to me that you would do better to move this question to Mathematics.SE. The issue is how to transform the integrand into a form that can be expanded in ρ, which is not a Mathematica.SE question. Once that is solved, doing the expansion with Mathematica should be straightforward. $\endgroup$
    – bbgodfrey
    Commented Feb 26, 2015 at 15:43
  • $\begingroup$ Apologies, I forgot to define $\rho=\sigma\sqrt{y}$ which is now corrected. $\endgroup$
    – Till Hoffmann
    Commented Feb 26, 2015 at 15:49

1 Answer 1

Reset to default
7
$\begingroup$

The expression in question is (we have replaced sigma by s)

g = Integrate[
  Exp[-(x - y)^2/(2 y s^2)] f[y] 1/Sqrt[2 \[Pi] y s^2], {y, 0, \[Infinity]}]

$\int_0^{\infty } \frac{e^{-\frac{(x-y)^2}{2 s^2 y}} f[y]}{\sqrt{2 \pi } \sqrt{s^2 y}} \, dy$

First of all we notice that for y>0

"Limit[Exp[-(x - y)^2/(2 y s^2)] 1/Sqrt[2 \[Pi] y s^2], s -> 0] = 
 DiracDelta[x - y]"

Formally

Limit[Exp[-(x - y)^2/(2 y s^2)] 1/Sqrt[2 \[Pi] y s^2], 
 s -> 0, Assumptions -> {y > 0, x > 0, y != x}]

(*
Out[37]= 0
*)

In[39]:= Limit[Exp[-(x - y)^2/(2 y s^2)] 1/Sqrt[2 \[Pi] y s^2], 
 s -> 0, Assumptions -> {y > 0, x > 0, y == x}]

(*
Out[39]= \[Infinity]
*)

so that

Limit[g, s -> 0] = f[x]

Assuming now that f[y] can be expanded into a power series about y = 0 and consider

g1 = Integrate[
  Exp[-(x - y)^2/(2 y s^2)] y^k/Sqrt[2 \[Pi] y s^2], {y, 0, \[Infinity]}, 
  Assumptions -> {s > 0, x^2 > 0}]

(*
Out[16]= (E^(x/s^2) Sqrt[2/\[Pi]] Abs[x]^(1/2 + k) BesselK[1/2 + k, Abs[x]/s^2])/s
*)

Now the series expansion about s = 0 becomes

sg1 = Series[g1, {s, 0, 2}, Assumptions -> s > 0] // Normal

(*
Out[27]= E^(x/s^2 - Abs[x]/s^2) (1/2 k (1 + k) s^2 Abs[x]^(-1 + k) + Abs[x]^k)
*)

Let us distinguish the region x<0 and x>0:

Simplify[sg1, x < 0]

(*
Out[28]= 1/2 E^((2 x)/s^2) (k (1 + k) s^2 - 2 x) (-x)^(-1 + k)
*)

Obviously, the exponential factor cannot be expanded about s=0.

Simplify[sg1, x > 0]

(*
Out[31]= 1/2 x^(-1 + k) (k s^2 + k^2 s^2 + 2 x)
*)

In the limit

Limit[%, s -> 0]

(*
Out[32]= x^k
*)

i.e. we rediscover the effect of the delta function.

Edit (Till Hoffmann): I thought I'd add to this answer to complete it.

Note that the expansion $\mathtt{sg1}=x^k + \sigma^2 \frac{k(1+k)}{2} x^{k-1}$ can be related to the Taylor expansion of $$ f(x) + \frac{\sigma^2}{2} \frac{\partial^2}{\partial x^2} \left(x f(x)\right), $$ which is the desired series expansion in terms of the original function.

EDIT (Dr. Wolfgang Hintze, 28.02.15)

Let me add the general expansion in the form of Till Hoffmann, which I call s-expansion in the following

$$fs = \sum _{j=0}^{\infty } \frac{s^{2 j} }{j!2^{j}}\frac{\partial ^{2 j}\left(f(x) x^j\right)}{\partial x^{2 j}}$$

It is very interesting that this expression seems to hold not only for the powers of y considered up to now but for more general functions f(y).

In order to verify this hypothesis we shall study examples for other functions than simple powers of the s-expansion fs of our integral.

The s-expansion with a function f is given by

fs[x_, s_, f_, n_] := 
 Sum[s^(2 j) 1/(j! 2^j) D[x^(j) f, {x, 2 j}], {j, 0, n}]

Example 1: Simple pole on the negative real y-axis

The integral is numerically

gg1[x_, s_] := 
 NIntegrate[
  Exp[-(x - y)^2/(2 y s^2)] 1/(1 + y) 1/Sqrt[2 \[Pi] y s^2], {y,0, \[Infinity]}]

The first 3 terms of the s-expansion are

fs1[x_, s_] = Table[fs[x, s, 1/(1 + x), k], {k, 0, 2}] // Simplify

(*
{1/(1 + x), (-s^2 + (1 + x)^2)/(1 + x)^3, (
     3 s^4 - s^2 (1 + x)^2 + (1 + x)^4)/(1 + x)^5} 
*)

Graphs up to order s^4 for two values of s are

With[{s = 1}, 
 Plot[{gg1[x, s], {1/(1 + x), (-s^2 + (1 + x)^2)/(1 + x)^3, (
    3 s^4 - s^2 (1 + x)^2 + (1 + x)^4)/(1 + x)^5}}, {x, -1, 2}, 
  PlotRange -> {0, 1.5}, ImageSize -> 400, 
  PlotLabel -> 
   Style["s-expansion of integral\nf(y) = 1/(1+y), s = " <> ToString[s] <> 
     "\n", 14], AxesLabel -> {"x", "gg1[x]"}, 
  Epilog -> {Text[
     Style["Legend of curves:\nblue = integral, numeric (gg1)\nred = f(x)\n\
brown = O(\!\(\*SuperscriptBox[\(s\), \(2\)]\))\ngreen = \
O(\!\(\*SuperscriptBox[\(s\), \(4\)]\))", Medium], {1, 1.2}]}]]
(* 150228_s-expansion_f1 _s1.jpg *)

enter image description here

With[{s = 0.3}, 
 Plot[{gg1[x, s], {1/(1 + x), (-s^2 + (1 + x)^2)/(1 + x)^3, (
    3 s^4 - s^2 (1 + x)^2 + (1 + x)^4)/(1 + x)^5}}, {x, -1, 2}, 
  PlotRange -> {0, 1.5}, ImageSize -> 400, 
  PlotLabel -> 
   Style["s-expansion of integral\nf(y) = 1/(1+y), s = " <> 
     ToString[s] <> "\n", 14], AxesLabel -> {"x", "gg1[x]"}, 
  Epilog -> {Text[
     Style["Legend of curves:\nblue = integral, numeric (gg1)\nred = \
f(x)\nbrown = O(\!\(\*SuperscriptBox[\(s\), \(2\)]\))\ngreen = O(\!\(\
\*SuperscriptBox[\(s\), \(4\)]\))", Medium], {1, 1.2}]}]]
(* 150228_s-expansion_f1 _s0-3.jpg *)

enter image description here

Example 2: Simple conjugate poles on the imaginary y-axis

The integral is numerically

gg2[x_, s_] := 
 NIntegrate[
  Exp[-(x - y)^2/(2 y s^2)] 1/(1 + y^2) 1/Sqrt[2 \[Pi] y s^2], {y, 0, \[Infinity]}]

The first 3 terms of the s-expansion are

fs2[x, s]

(*
{1/(1 + x^2), (s^2 x (-3 + x^2) + (1 + x^2)^2)/(1 + x^2)^3, (
 s^2 x (-3 + x^2) (1 + x^2)^2 + (1 + x^2)^4 - 
  3 s^4 (1 - 10 x^2 + 5 x^4))/(1 + x^2)^5} 
*)

Graphs up to order s^4 for two values of s are

With[{s = 1}, 
 Plot[{gg2[x, s], {1/(1 + x^2), (
    s^2 x (-3 + x^2) + (1 + x^2)^2)/(1 + x^2)^3, (
    s^2 x (-3 + x^2) (1 + x^2)^2 + (1 + x^2)^4 - 
     3 s^4 (1 - 10 x^2 + 5 x^4))/(1 + x^2)^5}}, {x, -1, 2}, 
  PlotRange -> {0, 1.5}, ImageSize -> 400, 
  PlotLabel -> 
   Style["s-expansion of integral\nf(y) = 1/(1+\!\(\*SuperscriptBox[\(y\), \
\(2\)]\)), s = " <> ToString[s] <> "\n", 14], AxesLabel -> {"x", "gg2[x]"}, 
  Epilog -> {Text[
     Style["Legend of curves:\nblue = integral, numeric (gg1)\nred = f(x)\n\
brown = O(\!\(\*SuperscriptBox[\(s\), \(2\)]\))\ngreen = \
O(\!\(\*SuperscriptBox[\(s\), \(4\)]\))", Medium], {1, 1.2}]}]]
(* 150228_s-expansion_f2 _s0-9.jpg *)

enter image description here

With[{s = 0.4}, 
 Plot[{gg2[x, s], {1/(1 + x^2), (
    s^2 x (-3 + x^2) + (1 + x^2)^2)/(1 + x^2)^3, (
    s^2 x (-3 + x^2) (1 + x^2)^2 + (1 + x^2)^4 - 
     3 s^4 (1 - 10 x^2 + 5 x^4))/(1 + x^2)^5}}, {x, -1, 2}, 
  PlotRange -> {0, 1.5}, ImageSize -> 400, 
  PlotLabel -> 
   Style["s-expansion of integral\nf(y) = \
1/(1+\!\(\*SuperscriptBox[\(y\), \(2\)]\)), s = " <> ToString[s] <> 
     "\n", 14], AxesLabel -> {"x", "gg2[x]"}, 
  Epilog -> {Text[
     Style["Legend of curves:\nblue = integral, numeric (gg1)\nred = \
f(x)\nbrown = O(\!\(\*SuperscriptBox[\(s\), \(2\)]\))\ngreen = O(\!\(\
\*SuperscriptBox[\(s\), \(4\)]\))", Medium], {1, 1.2}]}]]
(* 150228_s-expansion_f2 _s0-4.jpg *)

enter image description here

We see that the agreement is what is to expected in the sense of an asymptotic expansion.

Conclusion: the formula for $fs$ seems to be valid for a broader class of function.

Improve this answer
$\endgroup$
5
  • $\begingroup$ Hi, thanks for the suggestion. Using the expansion $\mathtt{sg1}=x^k + \sigma^2 \frac{k(1+k)}{2} x^{k-1}$, the first term reproduces the desired function $f(x)$ in a Taylor expansion. How can I express the second term in terms of the original function and its derivatives? It almost looks like a second derivative. I'll have a play with it. $\endgroup$
    – Till Hoffmann
    Commented Feb 26, 2015 at 13:58
  • $\begingroup$ Dr. Wolfgang Hintze, thank you very much for your help. I would like to mention you in the acknowledgements of the paper for which I required the above series expansion. Please drop me a line at t.hoffmann13(AT)imperial(DOT)ac(DOT)uk and I will send you a draft before publication. $\endgroup$
    – Till Hoffmann
    Commented Feb 26, 2015 at 17:19
  • $\begingroup$ Although this is an elegant solution, it may be valid only if f is an entire function, so that its series converges throughout the range of integration. If not, then some additional terms associated with the poles of f in the complex plane may arise. $\endgroup$
    – bbgodfrey
    Commented Feb 26, 2015 at 17:34
  • 1
    $\begingroup$ @bbgodfrey I agree completely. I tried e.g. two poles f[y] = 1/(1+y^2), without success. The y-integral must be split into one over (0..1) and the other over (1,oo). But the resulting integrals had no closed expressions. $\endgroup$
    – Dr. Wolfgang Hintze
    Commented Feb 26, 2015 at 21:03
  • $\begingroup$ @bbgodfrey Thank you for the hint which led me to consider other functions than powers successfully. See my edit of today. $\endgroup$
    – Dr. Wolfgang Hintze
    Commented Feb 28, 2015 at 21:02

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.