5

Is there a fast method to get the coefficients of a Taylor series expansion of the 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.

hlren
  • 117
  • 8

2 Answers2

6

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.

hlren
  • 117
  • 8
1

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]]
  • 1
    Thanks. You are the expert in mathematica. But I am new to mathematica. It would be great if some examples are added. – hlren Jul 08 '20 at 14:27