How to output Taylor formula in this format?

Question

I want to write a custom function that outputs the Taylor series as follows:

$$f(X)=f\left(X^{(0)}\right)+\nabla f\left(X^{(0)}\right)^{T} \Delta X+\frac{1}{2} \Delta X^{T} G\left(X^{(0)}\right) \Delta X+\cdots $$

$$G\left(X^{(0)}\right)=\left.\left(\begin{array}{cc} \frac{\partial^{2} f}{\partial x_{1}^{2}} & \frac{\partial^{2} f}{\partial x_{1} \partial x_{2}} \\ \frac{\partial^{2} f}{\partial x_{2} \partial x_{1}} & \frac{\partial^{2} f}{\partial x_{2}^{2}} \end{array}\right)\right|_{X^{(0)}}, \Delta X=\left(\begin{array}{c} \Delta x_{1} \\ \Delta x_{2} \end{array}\right)$$

That is to say, what I need is the following format that keeps the matrix from being computed:

$$f(X)=f\left(X^{(0)}\right)+\left(\frac{\partial f}{\partial x_{1}}, \frac{\partial f}{\partial x_{2}}\right)_{X^{(0)}}\left(\begin{array}{c} \Delta x_{1} \\ \Delta x_{2} \end{array}\right)+\left.\frac{1}{2}\left(\begin{array}{ccc} \frac{\partial^{2} f}{\partial x_{1}^{2}} & \frac{\partial^{2} f}{\partial x_{1} \partial x_{2}} \\ \frac{\partial^{2} f}{\partial x_{2} \partial x_{1}} & \frac{\partial^{2} f}{\partial x_{2}^{2}} \end{array}\right)\right|_{X^{(0)}}\left(\begin{array}{c} \Delta x_{1} \\ \Delta x_{2} \end{array}\right)+\cdots $$

(I need to keep the matrix $\left.\left(\begin{array}{cc} \frac{\partial^{2} f}{\partial x_{1}^{2}} & \frac{\partial^{2} f}{\partial x_{1} \partial x_{2}} \\ \frac{\partial^{2} f}{\partial x_{2} \partial x_{1}} & \frac{\partial^{2} f}{\partial x_{2}^{2}} \end{array}\right)\right|_{X^{(0)}}$ and vector $\left(\begin{array}{c} \Delta x_{1} \\ \Delta x_{2} \end{array}\right) $ in a non operational format).

The user-defined function should be universal and can be applied to the Taylor series of any n-variable function.

F[k_] := Sum[
  Binomial[k, r]*Δx^r*Δy^(k - r)*
   Derivative[r, k - r][f][x0, y0], {r, 0, k}]
Expand[Sum[F[i]/i!, {i, 0, 3}]]
grads = NestList[Grad[#1, {x1, x2, x3}] & , f[x1, x2, x3], 
    6] /. {x1 -> 0, x2 -> 0, x3 -> 0}; 
[DoubleStruckCapitalX] = {Δx1, Δx2, 

Δx3}; 
Expand[(1/0!)grads[[1]] + (1/1!)
   grads[[2]] . [DoubleStruckCapitalX] + (1/2!)*
   grads[[3]] . [DoubleStruckCapitalX] . [DoubleStruckCapitalX] + 
     (1/3!)*
   grads[[4]] . [DoubleStruckCapitalX] . [DoubleStruckCapitalX] . 

[DoubleStruckCapitalX]]

However, the results of the above codes can not keep the matrix and vector format. What can I do to get the output format I want?

Answer 1

Inactivate seems to be the function you are looking for

With

\[DoubleStruckCapitalX] = {\[CapitalDelta]x1, \[CapitalDelta]x2};
x = Transpose[{\[DoubleStruckCapitalX]}];
grads = NestList[Grad[#1,{x1, x2}] &, f[x1, x2], 2];
g2 = grads[[2]];
g3 = grads[[3]];

I get

f[x0] + Inactivate[g2.x] + 1/2 Inactivate[x. g3.x]

If you want to add the subscript X[0] you could try and modify this post

Answer 2

Here is not an answer but a review. Sinc for higher order,there are not simple matrix expression,we had to use tensor.Here is a test for Taylor expand for two variables and three orders.

Clear["`*"];
x0 = {0, 0};
x={p,q};
f[{u_, v_}] := Exp[u*v];
{1/0! (D[f[x], {x, 0}] /. Thread[x -> x0]), 
 1/1! (D[f[x], {x, 1}] /. Thread[x -> x0]).(x - x0), 
 1/2! (D[f[x], {x, 2}] /. Thread[x -> x0]).(x - x0).(x - x0), 
 1/3! (D[f[x], {x, 3}] /. Thread[x -> x0]).(x - x0).(x - x0).(x - x0)}
Total@%

Updated

taylorPoly[f_, x_, x0_, n_] := 
  Sum[(D[f, {x, i}] /. Thread[x -> x0]).Sequence @@ Table[x - x0, i]/
   i!, {i, 0, n}];
taylorPoly[Exp[x*y],{x,y},{0,0},3]