3

I'm trying to get maxima to perform some "abstract" Taylor series expansions, and I'm running into a simplification issue. A prototype of the problem might be the finite-difference analog of the gradient,

g(x1,dx1) := (f(x1+dx1) - f(x1))/dx1;  /* dx1 is small */
taylor(g(x1,dx1), [dx1], [0], 0);

for which maxima returns

enter image description here

So far so good. But now try the finite-difference analog of the second derivative (Hessian),

h(x1,dx1) := (f(x1+dx1) - 2*f(x1) + f(x1-dx1))/dx1^2;
taylor(h(x1,dx1), dx1, 0, 0);

for which I get

enter image description here

which is not nearly as helpful.

A prototype of the "real" problem I want to solve is to compute the low-order errors of the finite-difference approximation to ∂^2 f/(∂x1 ∂x2),

(f(x1+dx1, x2+dx2) - f(x1+dx1, x2) - f(x1, x2+dx2) + f(x1, x2))/(dx1*dx2)

and to collect the terms up to second order (which involves up to 4th derivatives of f). Without reasonably effective simplification I suspect it will be easier to do by hand than by computer algebra, so I am wondering what can be done to coax maxima into doing the simplification for me.

tholy
  • 11,882
  • 1
  • 29
  • 42
  • This is a good question, but my guess is that it's more suitable for the Maxima discussion list. See: https://sourceforge.net/projects/maxima/lists/maxima-discuss Also, can you clarify what you mean by "not nearly as helpful"? What are the simplifications which you are hoping to see which are not being carried out? – Robert Dodier Jun 04 '18 at 15:01
  • Those two first derivatives cancel one another because of the opposite signs of the `dx1` terms. Likewise, the two second derivatives should add (because `-1 * -1 = 1`), resulting in just `d^2 f/dx^2`. – tholy Jun 07 '18 at 13:00

1 Answers1

1

Consider this example. It uses Barton Willis' pdiff package. I simplified notation a bit: moved center to [0, 0] and introduced notation for partial derivatives.

(%i1) load("pdiff") $
(%i2) matchdeclare([n, m], integerp) $
(%i3) tellsimpafter(f(0, 0), 'f00) $
(%i4) tellsimpafter(pderivop(f,n,m)(0,0), concat('f, n, m)) $
(%i5) e: (f(dx, dy) - f(dx, -dy) - f(-dx, dy) + f(-dx, -dy))/(4*dx*dy)$
(%i6) taylor(e, [dx, dy], [0, 0], 3);
                                    2         2
                              f31 dx  + f13 dy
(%o6)/T/                f11 + ----------------- + . . .
                                      6
slitvinov
  • 5,693
  • 20
  • 31