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 ])
1
1
#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)
--
Mario
Mario Rodriguez Riotorto
2010-08-05
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.
--
Mario
Dieter Kaiser
2010-08-06
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.
Remark:
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
2010-08-06
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.
Remark:
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
Mario Rodriguez Riotorto
2010-08-06
Mario Rodriguez Riotorto
2010-08-06
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.
--
Mario