Good coding for a tricky expansion with convolutions, limits, substitutions

Question

I am attempting to generate a power series in $v$ as $v\rightarrow0$ for

$$g(x,v)\equiv \ln\int_{-\infty}^{\infty}e^{-y\left(z\right)}N\left(x-z,v\right)\,\mathrm{d}z $$

where $y(x)\equiv\sum_{i=1}^{n}c_{i}x^{i}$, $n$ is even, and $c_{n}>0$ so the integral converges. $N(x,v)$ is the PDF of the Normal distribution centered in $x$ with variance $v$.

The key mathematical step that I can't reasonably expect Mathematica to figure out is:

$$\frac{d}{dv}\int f(z)N(z-x,v)dz=\int f(z)\frac{d}{dv}N(z-x,v)dz=\int f(z)\frac{1}{2}\frac{d^{2}}{d^{2}x}N(z-x,v)dz=\frac{1}{2}\frac{d^{2}}{d^{2}x}\int f(z)N(z-x,v)dz\rightarrow\frac{1}{2}\frac{d^{2}}{d^{2}x}\int f(z)\delta(z-x)=\frac{1}{2}\frac{d^{2}}{d^{2}x}f(x)$$

I've developed some truly ugly code that does the expansion of $g[x,v]$ in $v$, but it's an unreadable mess that relies on my doing most of the manipulations on paper beforehand, and I have to rewrite most of it every time I want to change the calculation a bit. How would a Mathematicagician more experienced than me do this calculation?

I only need the first several terms in the expansion in $v$, not a general expression. But I'll be impressed all day if someone comes up with a general expression.

Answer 1

We consider h = Exp[g].

The main idea is to calculate the Fourier transform of h, expand it into a series in powers of v, and then transform back

Here we go.

The normal distribution is explicitly

f = 1/(v Sqrt[2 \[Pi]]) Exp[-(x - z)^2/(2 v^2)];

Now the Fourier transformation of h acts only on f, with the result

ft = FourierTransform[f, x, t, Assumptions -> v > 0]

(*
Out[49]= E^(-(1/2) t^2 v^2 + I t z)/Sqrt[2 \[Pi]]
*)

For simplicity let (with a > 0)

y[z_] := a z^2

so that the convolution becomes

cft = Integrate[
  Exp[-y[z]] E^(-(1/2) t^2 v^2 + I t z)/Sqrt[
   2 \[Pi]], {z, -\[Infinity], \[Infinity]}, Assumptions -> a > 0]

(*
Out[51]= E^(-((t^2 (1 + 2 a v^2))/(4 a)))/(Sqrt[2] Sqrt[a])
*)

Now we expand in powers of v

scft = (Series[cft, {v, 0, 4}] // Normal)

(*
Out[53]= E^(-(t^2/(4 a)))/(Sqrt[2] Sqrt[a]) - (E^(-(t^2/(4 a))) t^2 v^2)/(
 2 Sqrt[2] Sqrt[a]) + (E^(-(t^2/(4 a))) t^4 v^4)/(8 Sqrt[2] Sqrt[a])
*)

and transform back

h = InverseFourierTransform[scft, t, x, Assumptions -> a > 0]

(*
Out[54]= 1/2 E^(-a x^2) (2 - 2 a v^2 - 12 a^3 v^4 x^2 + 4 a^4 v^4 x^4 + 
   a^2 (3 v^4 + 4 v^2 x^2))
*)

For general y[z] we don't calculate the z integral in the first place but write for the Fourier transform of h

hf = Integrate[
  Exp[-y[z]] Exp[-(1/2) t^2 v^2 + I t z], {z, -\[Infinity], \[Infinity]}]

Now we see that the z-integral just gives the Fourier transform (Ft) of Exp[g] giving

hf = Ft[Exp[-y[z]] , z, t] Exp[-(1/2) t^2 v^2]

This must be transformed back (iFt), giving

h = iFt[ Ft[Exp[-y[z]] , z, t] Exp[-(1/2) t^2 v^2], t, x]

The expansion into powers of v is obvious, giving powers in t^2. The final g = Log[h] can then be expanded into powers of v, correspondingly.

I hope this is the general expression you were looking for.

Answer 2

I am the OP.

Below is the code I came up with that grinds out the result. Basically, because I couldn't figure out how to get MMa to do the substitutions and limits, I just did as much of the problem as I could on paper and found this shortcut:

$$\frac{dg[x,v]}{dv}=\frac{1}{2}\frac{d^{2}g}{dx^{2}}+\frac{1}{2}\left(\frac{dg}{dx}\right)^{2}$$

So this gets the job done but I had to do a lot of trying stuff on paper to get it. This was my hack "around" not knowing how to get MMa to do all that for me.

y[x_] = Sum[c[j]*x^j, {j, 0, 4}];
ord = 5;
h = Table[0, {i, 1, ord}];
h[[1]] = g[x, v];
Do[
  h[[n]] = Expand[(D[h[[n - 1]], v]) /. Derivative[q_, 1][g][x, v] -> D[(D[g[x, v], {x, 2}] + D[g[x, v], x]^2)/2, {x, q}]], {n, 2,ord}];
hh = h /. Derivative[q_, 0][g][x, v] -> D[y[x], {x, q}] /. g[x, v] -> y[x] // Expand;