I want to evaluate a term like this,
$$\left.\exp\left({1\over 2}\sum_{i,j=1}^{n}(A^{-1})_{ij}{\partial \over \partial x_i}{\partial \over \partial x_j}\right) f(\vec{x})\right|_{\vec{x}=0}$$
and I want to make a function that can evaluate this term given an arbitrary polynomial function $f$, with an arbitrary number of variables $x_i$.
The trouble I'm having, is compactly writing out all the derivatives that will survive the series expansion. Since the function $f$ is polynomial, we only need to take the series expansion out to a finite number of terms, equal to $\lfloor d/2\rfloor$ where $d$ is the degree of the polynomial.
$$f + \frac{1}{2} \sum_{i,j=1}^{n}(A^{-1})_{ij}{\partial \over \partial x_i}{\partial \over \partial x_j}f + \frac{1}{8}\left( \sum_{i,j=1}^{n}(A^{-1})_{ij}{\partial \over \partial x_i}{\partial \over \partial x_j}\right)^2f + ....$$
Now the mixed partial derivatives of $f$ will be symmetric and since the matrix $A$ is also symmetric, the $(A)^{-1})_{ij}=(A)^{-1})_{ji}$, but I haven't taken advantage of that in my code below.
expD[f_, vars_, ainverse_] :=
Module[{combinations, nvars, fdegree},
fdegree =
Exponent[# /. ((# -> \[FormalX] RandomReal[]) & /@
Variables[#]), \[FormalX]] &@f;
nvars = Length@vars;
(f + Sum[combinations = Tuples[Range@Length@vars, 2 m];
1/(2^m m!)
Sum[Times @@ (Extract[ainverse, #] & /@
Partition[n, 2]) (Fold[D, f, vars[[#]] & /@ n]), {n,
combinations}],
{m, 1, Floor[fdegree/2]}]) /.
Inner[Rule, vars, ConstantArray[0, nvars], List]
]
Thanks to klgr for the code to find the degree of the polynomial. You can verify that it works by using a test function and a test matrix,
testf = 4 + 8 x^4 + 12 x^3 y - y^2 x;
testainv = {{7/2, -1}, {-1, 1/2}};
and comparing the output of this function to what you get if you manually write it out,
expD[testf, {x, y}, testainv]
(testf +
1/2 (testainv[[1, 1]] D[testf, {x, 2}, {y, 0}] +
2 testainv[[1, 2]] D[testf, {x, 1}, {y, 1}] +
testainv[[2, 2]] D[testf, {x, 0}, {y, 2}]) +
1/8 (testainv[[1, 1]] testainv[[1, 1]] D[testf, {x, 4}, {y, 0}] +
4 testainv[[1, 1]] testainv[[1, 2]] D[
testf, {x, 3}, {y, 1}])) /. {x -> 0, y -> 0}
(* 172 *)
(* 172 *)
NestList
for example. $\endgroup$