A short time ago, there was some e-mail about doing more numerical
stuff in maxima. (There's also an e-mail from around mid June, 2005
about this).
Here is a simple proof-of-concept for this idea for Minpack. The
example is taken from the documentation for lmdif. The problem is to
determine the values of x(1), x(2), and x(3) which provides the best
least-squares fit of
x(1) + u(i)/(v(i)*x(2) + w(i)*x(3)), i = 1 to 15
to the data
y = (0.14,0.18,0.22,0.25,0.29,0.32,0.35,0.39,
0.37,0.58,0.73,0.96,1.34,2.10,4.39),
where u(i) = i, v(i) = 16 - i, and w(i) = min(u(i),v(i)).
Assuming you've loaded the minpack package (obtained from f2cl), you
can compile this Lisp function to give access to lmdif:
(in-package :maxima)
(defun $lmdif (fcn m n x tol)
(let* ((lwa (+ (* m n)
(* 5 n)
m))
(iwa (make-array m :element-type '(signed-byte 32)))
(wa (make-array lwa :element-type 'double-float))
(x (make-array n :element-type 'double-float
:initial-contents
(mapcar #'(lambda (r)
(float r 1d0))
(cdr x))))
(fvec (make-array m :element-type 'double-float))
(info 0))
(flet ((f (m n x fvec iflag)
#+nil
(progn
(format t "m,n = ~A ~A~%" m n)
(format t "x = ~A~%" x)
(format t "x = ~A~%" (list* '(mlist) (coerce x 'list))))
(let ((res (mfuncall fcn (list* '(mlist) (coerce x 'list)))))
#+nil
(format t "res = ~A~%" res)
(map-into fvec #'(lambda (r)
($float r))
(cdr res))
(values m n x fvec iflag))))
(multiple-value-bind (z-ret z-m z-n z-x z-fvec z-tol z-info
z-iwa z-wa z-lwa)
(minpack::lmdif1 #'f m n x fvec tol info iwa wa lwa)
(declare (ignore z-ret z-m z-n z-x z-fvec z-tol z-iwa z-wa z-lwa))
(list '(mlist)
(list* '(mlist) (coerce x 'list))
(list* '(mlist) (coerce fvec 'list))
z-info)))))
Then you can do the following:
(%i54) fcn(x) :=
maplist(lambda([y],y[1]-(x[1]+y[2]/(y[3]*x[2]+min(y[2],y[3])*x[3]))),
[[0.14,1,15],[0.18,2,14],[0.22,3,13],[0.25,4,12],[0.29,5,11],[0.32,6,10],[0.35,7,9],[0.39,8,8],[0.37,9,7],[0.58,10,6],[0.73,11,5],[0.96,12,4],[1.34,13,3],[2.10,14,2],[4.39,15,1]])$
(%l55) lmdif(fcn, 15, 3, [1.0,1.0,1.0], 1d-8);
(%o55) [[.08241057586436833, 1.133036626537844, 2.343694664907586],
[.005881094857421248, 2.6535951786882395e-4, - 2.746740648905399e-4,
- .006541523948618011, 8.230030670372535e-4, .001299499711875207,
0.00446310745705325, .01996293932200199, - .08221605633903506,
.01821194944624882, .01481115754564644, 0.0147099696947427,
.01120798994323646, .004204030440223772, - .006807848068813627],
1]
A list of results is returned. The first element is the estimated
solution: [.08241057586436833, 1.133036626537844, 2.343694664907586].
The expected result from minpack is 0.8241057D-01, 0.1133037D+01,
0.2343695D+01.
The second element is the value of fcn at the estimated solution.
The third element is 1 to indicate success.
Ray