[Maxima-discuss] Completing the Square in N dimensions

[Mike V. <mic...@gm...>]

Hello all,

I have a bit of code (included below) for completing the square in
N-dimensions. Three questions:

1) I am just unsure what I should have the method do in the case where the
matrix of second order terms is non-invertible (see code). In these cases,
factorsum actually appears to work well, so maybe I should have it return
that result instead?

2) What is the most elegant way I could extend to this to operate on
general expressions where only the second order terms get rearranged. One
solution, although not necessarily the best, is to take the original
expression, subtract my result, simplify, add my result. Maxima might have
a nicer way. Maybe I should spend more effort parsing the expression?

3) What is the most elegant way to get scanmap to operate on methods having
the headings like method(expression, [vars]). The most straightforward way
I can think of is to use a lambda function: scanmap(expression,
lambda([expression], mymethod(expression, [vars])).



The following, incomplete code is released under the third Gnu General
Public License (GPL3).
-----
/* This file is released under the GNU General Public License 3 */
/* http://www.gnu.org/licenses/gpl-3.0.en.html */

/* This is the basic outline for completing the square in N dimensions */
/* I would appreciate it if someone can do some better error checking */
/* known issue when c_mat is non-invertible:
 complete_square(expand(x^2+y^2), [x,y]) */

"Use: complete_square(Quadratic expression, [vars])";
"Option: complete_square_alt_form: 0, 1, or 2";
" 0) This returns 3 elements [v,C,a] where the completed square is 1/2 *
transpose(v) . C . v + a";
" 1) This mode carries out the matrix multiplication in 0";
" 2) This mode attempts to carry out the matrix multiplication in a
non-standard way to produce a nicer form";

complete_square_alt_form:0$
complete_square(expr, vars):=block([a_const:expr, b_vec, c_mat, c_mat_inv,
new_const, offset, x],
  /* Needs proper error checking */
  if(not listp(vars)) then error("vars must be a list of variables"),
  for i: 1 thru length(vars) do
    if(not atom(vars[i]) or subvarp(vars[i])) then error("The variables
must be atoms or indexed variables"),

  /* Allocate space */
  b_vec: genmatrix(lambda([i,j], 0), length(vars), 1),
  c_mat: genmatrix(lambda([i,j], 0), length(vars), length(vars)),

  /* Build the vector of coefficients of variables with degree 1, purge
cross-terms, and isolate the constant */
  for i: 1 thru length(vars) do block(
    b_vec[i,1]: coeff(expr, vars[i]),
    b_vec[i,1]: ev(b_vec[i,1], map(lambda([avar], avar=0), vars)),
    a_const: ratcoef(a_const, vars[i], 0)
  ),

  /* Build the matrix of coefficients for degree 2 */
  for i: 1 thru length(vars) do
    for j: 1 thru length(vars) do
      c_mat[i,j]: ratcoef(expr, vars[i]*vars[j])/(2-kron_delta(i,j)),
  /* c_mat should be symmetric, which should make this easier */

  /* Build the offset and new constant*/
  c_mat_inv:invert(c_mat),
  offset: ratsimp( -c_mat_inv . b_vec),
  new_const: ratsimp(a_const - transpose(b_vec) . c_mat_inv . b_vec / 2),
  x: ''(vars - offset),

  /* Output in the most pretty sense */
  if(length(vars)=1) then ((vars - offset)^2/2 + ratsimp(new_const))[1,1]
  elseif(complete_square_alt_form=2) then matrix_sum_elements((x .
transpose(x))*c_mat/2) + ratsimp(new_const)
  elseif(complete_square_alt_form=1) then transpose(x) . c_mat . x/2 +
ratsimp(new_const)
  else [x, c_mat, ratsimp(new_const)]
)$

matrix_sum_elements(amat):= block([res_sum:0],
  for i: 1 thru length(amat) do res_sum: res_sum + lsum(x,x,amat[i]),
  res_sum
)$