Re: [Maxima-discuss] Completing the Square in N dimensions

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

Oppps, also the answer is not quite correct - fixing that now too

On Sun, May 31, 2015 at 5:03 PM, Mike Valenzuela <mic...@gm...>
wrote:

> 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
> )$
>
>