5
$\begingroup$

Is there a fast method to get the coefficients of Taylor series expansion of function $f(x_1,x_2,...,x_d)$ with maximal summed partial derivative up to $n$, where $d,n$ can be relatively large (for example, $d=10,n=10$)?

For example, in mathematica, the series expansion

n = 1; Series[f[x, y], {x, 0, n}, {y, 0, n}]

$\left(f(0,0)+y f^{(0,1)}(0,0)+O\left(y^2\right)\right)+x \left(f^{(1,0)}(0,0)+f^{(1,1)}(0,0) y+O\left(y^2\right)\right)+O\left(x^2\right)$

I need the coefficient list $(1, x, y)$ for function list: ($f(0,0), f^{(1,0)}(0,0), f^{(0,1)}(0,0)$)

The conventional Series function is very slow, for example,

In[655]:=n = 1; Series[f[x1, x2, x3, x4, x5, x6, x7], {x1, 0, n}, {x2, 0, n}, {x3, 0, n}, {x4, 0, n}, {x5, 0, n}, {x6, 0, n}, {x7, 0, n}]; // AbsoluteTiming

Out[655]= {0.000259252, Null}

In[656]:= n = 2; Series[f[x1, x2, x3, x4, x5, x6, x7], {x1, 0, n}, {x2, 0, n}, {x3, 0, n}, {x4, 0, n}, {x5, 0, n}, {x6, 0, n}, {x7, 0, n}]; // AbsoluteTiming

Out[656]= {0.442849, Null}

In[657]:= n = 3; Series[f[x1, x2, x3, x4, x5, x6, x7], {x1, 0, n}, {x2, 0, n}, {x3, 0, n}, {x4, 0, n}, {x5, 0, n}, {x6, 0, n}, {x7, 0, n}]; // AbsoluteTiming

Out[657]= {2.95541, Null}

In[665]:= n = 5; Series[f[x1, x2, x3, x4, x5, x6, x7], {x1, 0, n}, {x2, 0, n}, {x3, 0, n}, {x4, 0, n}, {x5, 0, n}, {x6, 0, n}, {x7, 0, n}]; // AbsoluteTiming

Out[665]= {42.6434, Null}

The series function seems to be inadequate for obtaining the coefficient list and derivative list of function with large $(n,d)$.

I have found a solution in https://arxiv.org/abs/1905.02809

Mathematica code for Taylor series in several variables

Based on the multi-index, the Taylor series expansion of a multi-variable scalar function $ u(x_{1},...,x_{d})$ at $\mathbf {a} =(0,...,0)$ can be written as

$ u(x_{1},...,x_{d})=\sum _{(n_{1},...,n_{d})\in \alpha _{d}^{n}}{\frac {x_{1}^{n_{1}}...x_{d}^{n_{d}}}{n_{1}!...n_{d}!}}u_{,n_{1}...n_{d}}+O(r^{n+1})$ ,

with $ r={\sqrt {x_{1}^{2}+...+x_{d}^{2}}}$ , $ u_{,n_{1}...n_{d}}={\frac {\partial ^{n_{1}+...+n_{d}}u}{\partial x_{1}^{n_{1}}...\partial x_{d}^{n_{d}}}}|_{x_{1}=0,...,x_{d}=0}$ and $ \alpha _{d}^{n}$ is the set of multi-indexes in $d$ dimensions with maximal single index up to $n$ , i.e. $ \mathbf {\alpha } _{d}^{n}=\{(n_{1},...,n_{d})|0\leq \sum _{i=1}^{d}n_{i}\leq n,\,n_{i}\in \mathbb {N} ^{0},1\leq i\leq d\},$

which can be obtained by an efficient Mathematica code as

MultiIndexList[d_,n_]:=Block[{a,b,c},a=Subsets[Range[d+n],{d}];
Do[ c=a[[i]];b=c-1;b[[2;;]]-=c[[1;;-2]];a[[i]]=b,{i,Length[a]}];a];
(*note: d=number of spatial dimensions, n=maximal order of derivative*)

(*
(*Some Performance test*)
MultiIndexList[5, 20] // Length // AbsoluteTiming
MultiIndexList[5, 40] // Length // AbsoluteTiming
MultiIndexList[10, 10] // Length // AbsoluteTiming
MultiIndexList[200, 3] // Length // AbsoluteTiming
MultiIndexList[200, 1] // Length // AbsoluteTiming
(*
{0.5810 seconds,53130 terms}
{9.777 seconds,1221759 terms}
{1.583 seconds,184756 terms}
{89.013 seconds,1373701 terms}
{0.01451 seconds,201 terms}
*)
*)

For each multi-index $ (n_{1},...,n_{d})\in \alpha _{d}^{n}$ , the polynomial and partial derivative are

$ {\frac {x_{1}^{n_{1}}...x_{d}^{n_{d}}}{n_{1}!...n_{d}!}},\,\,u_{,n_{1}...n_{d}},\quad \forall (n_{1},...,n_{d})\in \alpha _{d}^{n}$

When the multi-index is written explicitly, the Taylor series expansion of the multi-variable function is straightforward.

$\endgroup$
  • $\begingroup$ Other than the function term, isn't the coefficient of a derivative term given by term /. Derivative[nn__][f][xx__] :> Multinomial[nn] Apply[Times, {xx}^{nn}]? $\endgroup$ – Michael E2 Jan 23 at 16:39
  • $\begingroup$ Yes, the coefficient is right. But this needs the list of nn at the first step $\endgroup$ – hlren Jan 23 at 21:34
  • 1
    $\begingroup$ But you know the list nn: i.stack.imgur.com/puvUO.png $\endgroup$ – Michael E2 Jan 26 at 21:36
  • 1
    $\begingroup$ Can you explain why Jens's answer in the other thread did not suit your needs? With[{n = 5}, Normal[Series[Apply[f, {x1, x2, x3, x4, x5, x6, x7} t], {t, 0, n}]] /. t -> 1] $\endgroup$ – J. M.'s technical difficulties Jan 28 at 13:22
  • 1
    $\begingroup$ Additionally, you can adapt the code from this answer to generate your required multi-indices: MultiIndexList[d_, n_] := Flatten[Table[Permutations[ip] - 1, {k, n}, {ip, IntegerPartitions[k + d, {d}]}], 2] $\endgroup$ – J. M.'s technical difficulties Jan 28 at 13:56
6
$\begingroup$

With some discussion with @J.M.willbebacksoon and @MichaelE2, I found the key to the Taylor series expansion of multi-variable function is the generation of the multi-index list. Two solutions for the multi-indice are

 MultiIndexList0[d_,n_]:=Block[{a,b,c},a=Subsets[Range[d+n],{d}];
 Do[ c=a[[i]];b=c-1;b[[2;;]]-=c[[1;;-2]];a[[i]]=b,{i,Length[a]}];a];

and

 MultiIndexList1[d_, n_] := Flatten[Table[Permutations[ip] - 1, {k, n}, {ip, 
 IntegerPartitions[k + d, {d}]}], 2]

MultiIndexList1 is proposed by @J.M.willbebacksoon. Both MultiIndexList0 and MultiIndexList1 are very concise and fast. MultiIndexList1 is more than 10 times faster than MultiIndexList0.

The mathematica functions based on two multi-indices functions for the Taylor series expansion of a one variable or multi-variable function can be written as

TaylorSeries0[f_, {X_, X0_, n_}] := 
  Block[{vars = Flatten[{X}], vars0 = Flatten[{X0}], d, alist, xlist, 
    dflist, xx}, d = Length[vars]; xx = vars - vars0; 
    alist = MultiIndexList0[d, n]; 
    xlist = Table[1/Apply[Times, Factorial[nn]] Apply[Times, xx^nn], {nn, alist}]; 
    dflist = Table[D[f, Sequence @@ Transpose[{vars, nn}]]/. (Rule @@@ ({vars, vars0}\[Transpose])), {nn,alist}];{xlist, dflist}];
TaylorSeries1[f_, {X_, X0_, n_}] := 
  Block[{vars = Flatten[{X}], vars0 = Flatten[{X0}], d, alist, xlist, 
   dflist, xx}, d = Length[vars]; xx = vars - vars0; 
   alist = MultiIndexList1[d, n]; 
   xlist = Table[ 1/Apply[Times, Factorial[nn]] Apply[Times, xx^nn], {nn, alist}]; 
   dflist =Table[D[f, Sequence @@ Transpose[{vars, nn}]]/. (Rule @@@ ({vars, vars0}\[Transpose])), {nn,alist}]; {xlist, dflist}];

The coefficients and derivatives are listed seperately.

Verfication:

Single variable function:

 In[84]:= TaylorSeries0[f[x, y], {x, 1, 5}] // TeXForm

 Out[84]//TeXForm=

$ \left( \begin{array}{cccccc} 1 & x-1 & \frac{1}{2} (x-1)^2 & \frac{1}{6} (x-1)^3 & \frac{1}{24} (x-1)^4 & \frac{1}{120} (x-1)^5 \\ f(1,y) & f^{(1,0)}(1,y) & f^{(2,0)}(1,y) & f^{(3,0)}(1,y) & f^{(4,0)}(1,y) & f^{(5,0)}(1,y) \\ \end{array} \right) $

Multi-variable function:

 In[82]:= TaylorSeries0[f[x, y], {{x, y}, {1, 1}, 2}] // TeXForm

 Out[82]//TeXForm=

$\left( \begin{array}{cccccc} 1 & y-1 & \frac{1}{2} (y-1)^2 & x-1 & (x-1) (y-1) & \frac{1}{2} (x-1)^2 \\ f(1,1) & f^{(0,1)}(1,1) & f^{(0,2)}(1,1) & f^{(1,0)}(1,1) & f^{(1,1)}(1,1) & f^{(2,0)}(1,1) \\ \end{array} \right)$

Some performance tests:

In[73]:= tf = TaylorSeries0[ f[x1, x2, x3, x4, x5, x6, x7], 
{{x1, x2, x3, x4, x5, x6, x7}, Table[0, {7}], 4}]; // AbsoluteTiming

Out[73]= {0.027772, Null}

In[76]:= tf =TaylorSeries0[f[x1, x2, x3, x4, x5, x6, x7], 
{{x1, x2, x3, x4, x5, x6, x7}, Table[0, {7}], 10}]; // AbsoluteTiming

Out[76]= {1.4329, Null}

In[77]:= tf =TaylorSeries1[f[x1, x2, x3, x4, x5, x6, x7],
{{x1, x2, x3, x4, x5, x6, x7},Table[0, {7}], 10}]; // AbsoluteTiming

Out[77]= {1.59547, Null}

In[75]:= tf =TaylorSeries1[f[x1, x2, x3, x4, x5, x6, x7],
 {{x1, x2, x3, x4, x5, x6, x7},Table[0, {7}], 4}]; // AbsoluteTiming

Out[75]= {0.0153569, Null}

In[55]:=tf = TaylorSeries0[f[x1, x2, x3, x4, x5, x6, x7], 
{{x1, x2, x3, x4, x5, x6, x7},Table[0, {7}], 20}]; // AbsoluteTiming

Out[55]= {84.3272, Null}

In[78]:= tf = TaylorSeries1[f[x1, x2, x3, x4, x5, x6, x7], 
{{x1, x2, x3, x4, x5, x6, x7},Table[0, {7}], 20}]; // AbsoluteTiming

Out[78]= {116.976, Null}

The TaylorSeries0 based on MultiIndexList0 seems to be slightly faster than TaylorSeries1.

| improve this answer | |
$\endgroup$
1
$\begingroup$

A simple answer:

(*ref/Grad  Maclaurin*)
orderedForm[poly_, var_List] := 
  HoldForm[+##] & @@ 
   MonomialList[poly, 
     var][[Ordering[-Total[#] & @@@ CoefficientRules[poly, var], All, 
      GreaterEqual]]];
(*f[x_,y_,z_] :=g[x,y,z]*)
grads = NestList[Grad[#, {x, y, z}] &, f[x, y, z], 6] /. {x -> 0, 
    y -> 0, z -> 0};
\[DoubleStruckCapitalX] = {x, y, z};

orderedForm[
 Sum[1/(k - 1)! Dot[grads[[k]], ##] & @@ 
   Table[\[DoubleStruckCapitalX], {k - 1}] , {k, 
   3}], \[DoubleStruckCapitalX]]
| 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.