#2042 mnewton fails with cvs maxima


Example taken from documentation:

(%i2) load("mnewton")$
(%i3) mnewton([2*a^a-5],[a],[1]);
Improper argument to coefmatrix:
(2.0 transpose(h ) + [ - 3.0 ])
#0: solve_by_lu(eqs=(2.0*'transpose(h[1])+matrix([-3.0]))[1],vars=[h[1]],field=complexfield)(mnewton.mac line 51)
#1: mnewton(funclist=[2*a^a-5],varlist=[a],guesslist=[1])(mnewton.mac line 113)
-- an error. To debug this try: debugmode(true);

But other examples from documentation work as expected:

(%i5) mnewton([2*3^u-v/u-5, u+2^v-4], [u, v], [2, 2]);
(%o5) [[u = 1.066618389595407, v = 1.552564766841786]]
(%i6) mnewton([x1+3*log(x1)-x2^2, 2*x1^2-x1*x2-5*x1+1],
[x1, x2], [5, 5]);
(%o6) [[x1 = 3.756834008012769, x2 = 2.779849592817897]]

(%i7) build_info();

Maxima version: 5.21post
Maxima build date: 21:2 7/26/2010
Host type: i686-pc-linux-gnu
Lisp implementation type: CLISP
Lisp implementation version: 2.44.1 (2008-02-23) (built 3476002983) (memory 3489159770)



  • I think the error is in simp.lisp.

    In Maxima 5.21 we had:
    (%i30) display2d:false$
    (%i31) s*r[1]+[-5];
    (%o31) [r[1]*s-5]

    Now, in Maxima 5.22:
    (%i38) display2d:false$
    (%i39) s*r[1]+[-5];
    (%o39) r[1]*s+[-5]

    This change produces an error in coefmatrix, when it is called from mnewton.


  • Dieter Kaiser
    Dieter Kaiser

    It is revision 1.108 of simp.lisp which causes the change. This revision is an attempt to implement the functions scalarp and nonscalarp more consistent. This is the expression which causes the bug:

    (2.0*transpose(h[1]) + [ - 3.0 ]);

    In Maxima 5.21 an expression transpose(h[1]) is handled like a scalar by default, because it is not a scalar and not a nonscalar. We get

    (%i3) scalarp(transpose(h[1]));
    (%o3) false
    (%i4) nonscalarp(transpose(h[1]));
    (%o4) false

    The expression transpose(h[1]) is handled like a scalar by default:

    (%i8) (2.0*transpose(h[1]) + [ - 3.0 ]);
    (%o8) [2.0*'transpose(h[1])-3.0]

    This has changed with revision 1.108 of simp.lisp. Now the expression transpose(h[1]) is a nonscalar:

    (%i4) scalarp(transpose(h[1]));
    (%o4) false
    (%i5) nonscalarp(transpose(h[1]));
    (%o5) true

    Because the expression transpose(h[1]) is not handled as a scalar, Maxima does not add elementwise:

    (%i8) (2.0*transpose(h[1]) + [ - 3.0 ]);
    (%o8) 2.0*'transpose(h[1])+[-3.0]

    I tend to say, that the new implementation is more correct and that the problem should be solved in mnewton. The algorithm of mnewton should not depend on the behavior that an expression transpose(h[1]) is handled as a scalar.


    The new behavior is to assume that a subscripted variable is a nonscalar:

    (%i13) scalarp(r[1]);
    (%o13) false
    (%i14) nonscalarp(r[1]);
    (%o14) true

    To get elementwise addition we have to declare r to be a scalar:

    (%i15) declare(r,scalar);
    (%o15) done
    (%i16) s*r[1]+[-5];
    (%o16) [r[1]*s-5]

    Dieter Kaiser

  • Dieter Kaiser
    Dieter Kaiser

    In addition:

    A workaround is to declare h to be a scalar in the function mnewton. With this declaration we get again the expected solutions. The testsuite and share_testsuite will have no problems.

    (%i2) mnewton([2*a^a-5],[a],[1]);
    (%o2) [[a = 1.70927556786144]]

    (%i3) mnewton([2*3^u-v/u-5, u+2^v-4], [u, v], [2, 2]);
    (%o3) [[u = 1.066618389595407, v = 1.552564766841787]]

    (%i4) mnewton([x1+3*log(x1)-x2^2, 2*x1^2-x1*x2-5*x1+1],[x1,x2],[5,5]);
    (%o4) [[x1 = 3.756834008012769, x2 = 2.779849592817897]]

    But it is clear that we have to thought through the behavior of scalarp and nonscalarp. The old behavior is to assume an expression to be a scalar by default. The more 'consistent ' but not perfect implementation has the problem, that an element h[1] of a matrix h is not a scalar by default.

    It is easy to change the implementation in a way that a subscripted symbol is a scalar by default. But we get problems with the tensor packages. These packages expect the behavior that a subscripted symbol is a matrix too.

    Dieter Kaiser

    • status: open --> closed
  • Thanks for the detailed clarifications. I didn't know this behavior was intended.

    I'll fix the problem in mnewton and close this bug report.