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

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.

• Off the top of my head, I would say that any Mathematicagician is likely to be listening to Mathemetallica. – march Sep 22 '15 at 6:11
• @march Oh God, you said it out loud. Now it's only a matter of time until Mathemetallica appears on YouTube. Enter Neumann... Call of Kolomogrov... Master of Publications... Simplify & Differentiate. It's only going to get worse from there. – Jerry Guern Sep 22 '15 at 6:58
• @JerryGuern and Fuel will keep its name, but it will be about coffee, instead! – rcollyer Sep 22 '15 at 10:42
• Now I'm afraid I'm going to get $2^{10}$ Mathemetallica jokes on this thread and no Answers to my Question. :-( Serves me right I guess. – Jerry Guern Sep 22 '15 at 11:24

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.

• Thanks! It's going to take me some time to follow this. – Jerry Guern Sep 22 '15 at 14:57
• Okay, I finally got to side down and work through this. I get it now, and it's hugely helpful. This isn't really a general expression in the sense of an analytic formula for the nth derive of g[x,v] wrt v, but I don't think that exists, and this does let me grind out as many terms of the expansion as I wish. Thanks for the assist! – Jerry Guern Sep 26 '15 at 8:29
• The Integrate[] fails if y[x] is higher order in x than x^2. – Jerry Guern Sep 26 '15 at 8:43
• @ Jerry Guern: of course, the general expression does not gurantee that MMA can solve the integrals symbolically. This must be done as usual with various methods depending on the case in question. Example: for y[z]=z^4 the Ft is computed by MMA and gives a combination of Gamma and Hypergeometric functions. In order to calculate h, an expansion into powers of t allows MMA to compute the iFt. – Dr. Wolfgang Hintze Sep 26 '15 at 8:59
• Dr. Hintze, I noticed you post good math solutions as well as MMa solutions. Would you care to take a crack at a bounty question I posted on the Math SE? math.stackexchange.com/questions/1455526/… – Jerry Guern Oct 1 '15 at 12:23

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;