7
$\begingroup$

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.

$\endgroup$
  • 3
    $\begingroup$ Off the top of my head, I would say that any Mathematicagician is likely to be listening to Mathemetallica. $\endgroup$ – march Sep 22 '15 at 6:11
  • 3
    $\begingroup$ @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. $\endgroup$ – Jerry Guern Sep 22 '15 at 6:58
  • 1
    $\begingroup$ @JerryGuern and Fuel will keep its name, but it will be about coffee, instead! $\endgroup$ – rcollyer Sep 22 '15 at 10:42
  • $\begingroup$ 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. $\endgroup$ – Jerry Guern Sep 22 '15 at 11:24
7
$\begingroup$

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.

| improve this answer | |
$\endgroup$
  • $\begingroup$ Thanks! It's going to take me some time to follow this. $\endgroup$ – Jerry Guern Sep 22 '15 at 14:57
  • $\begingroup$ 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! $\endgroup$ – Jerry Guern Sep 26 '15 at 8:29
  • $\begingroup$ The Integrate[] fails if y[x] is higher order in x than x^2. $\endgroup$ – Jerry Guern Sep 26 '15 at 8:43
  • $\begingroup$ @ 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. $\endgroup$ – Dr. Wolfgang Hintze Sep 26 '15 at 8:59
  • $\begingroup$ 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/… $\endgroup$ – Jerry Guern Oct 1 '15 at 12:23
2
$\begingroup$

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;
| improve this answer | |
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.