|
|
|
last edited 15 years ago by Mark Clements |
Edit detail for SandBoxNewtonsMethod revision 2 of 2
| 1 2 | ||
|
Editor: Mark Clements
Time: 2009/04/26 23:45:04 GMT-7 |
||
| Note: Included Waldek Hebisch's suggestion for mkLispFunction1 | ||
changed: -)lisp (defun newton (f dfdx x0 &optional (tol 1.0d-10)) (let ((xt x0) (fxt (funcall f x0)) (maxit 100) (iternum 0)) (loop (setf xt (- xt (/ fxt (funcall dfdx xt)))) (setf fxt (funcall f xt)) (if (< (abs fxt) tol) (return xt)) (incf iternum) (if (>= iternum maxit) (error "Maximum iterations exceeded."))))) )lisp (defun newton (f dfdx x0 &optional (tol 1.0d-10) (maxit 100)) (let ((xt x0) (fxt (funcall f x0)) (iternum 0)) (loop (setf xt (- xt (/ fxt (funcall dfdx xt)))) (setf fxt (funcall f xt)) (if (< (abs fxt) tol) (return xt)) (incf iternum) (if (>= iternum maxit) (error "Maximum iterations exceeded."))))) added: Finally, within a Spad package, it is even easier to call a function through Lisp: \begin{spad} -- the following approach is due to Waldek Hebisch )lisp (defun |mkLispFunction1| (f) (lambda (x) (SPADCALL x f))) )abbrev package TESTP TestPackage R ==> DoubleFloat I ==> Integer TestPackage: with newton:(R->R,R->R,R,R,I)->R newton:(R->R,R->R,R)->R == add newton(f,dfdx,x0,tol,maxit) == NEWTON(mkLispFunction1(f@(R -> R))$Lisp, mkLispFunction1(dfdx@(R -> R))$Lisp, x0,tol,maxit)$Lisp newton(f,dfdx,x0) == newton(f,dfdx,x0,1.0e-10::R,100::I) \end{spad} then \begin{axiom} R ==> DoubleFloat newton((x:R):R +-> x**2-2.0::R,(x:R):R +-> 2.0::R*x,3.0::SF,1.0e-15,100)$TestPackage-sqrt(2.0) newtonExpression(f,x,x0) == dfdx := D(f,x) newton((xt:R):R +-> eval(f,x,xt), (xt:R):R +-> eval(dfdx,x,xt), x0)$TestPackage newtonExpression(x**2-2.0,x,2.0)-sqrt(2.0) \end{axiom}
The following shows Newton's method for numerically solving f(x)=0. It is also shows examples of calling Axiom expressions and Spad functions from Lisp.
First, define Newton's method using the Axiom interpreter:
fricas
-- Axiom interpreter function for Newton's algorithm R ==> Float
Type: Void
fricas
I ==> Integer
Type: Void
fricas
newton(f:Expression R,x:Symbol, x0:R):R == dfdx:Expression R := D(f, x) xt:R := x0 fxt:R := subst(f, x=xt) iterNum:I := 0 maxIt:I := 100 repeat xt := xt - fxt/subst(dfdx, x=xt) fxt := subst(f, x=xt) if abs(fxt)<1.0e-10 then return xt iterNum:=iterNum+1::I if iterNum >= maxIt then error "Maximum iterations exceeded."
Function declaration newton : (Expression(Float),Symbol, Float) -> Float has been added to workspace.
Type: Void
fricas
newton(x**2-2.0,x, 2.0)
There are no library operations named ** Use HyperDoc Browse or issue )what op ** to learn if there is any operation containing " ** " in its name.
Cannot find a definition or applicable library operation named ** with argument type(s) Expression(Float) PositiveInteger
Perhaps you should use "@" to indicate the required return type,or "$" to specify which version of the function you need.
Second, we can also do this by calling the Newton method implemented in Lisp. We initially define/hack a Spad package which creates a Lisp function from an interpreter expression.
spad
)lisp (defun |lambdaFuncallSpad| (f) (lambda (x) (funcall f x nil))) )abbrev package MKULF MakeUnaryLispFunction ++ Tools for making compiled lisp functions from top-level expressions ++ Author: Mark Clements ++ (based on the MakeUnaryCompiledFunction package by Manuel Bronstein) ++ Date Created: 20 April 2008 ++ Date Last Updated: 20 April 2008 ++ Description: transforms top-level objects into lisp functions. MakeUnaryLispFunction(S,D, I): Exports == Implementation where S: ConvertibleTo InputForm D, I: Type
SY ==> Symbol DI ==> devaluate(D -> I)$Lisp
Exports ==> with compiledFunction: (S,SY) -> SY ++ compiledFunction(expr, x) returns a function lisp{f: D -> I} ++ defined by lisp{(defun f (x) expr)}. ++ Function f is compiled and directly ++ applicable to objects of type D (in lisp).
Implementation ==> add import MakeFunction(S)
compiledFunction(e:S,x:SY) == t := [convert([devaluate(D)$Lisp]$List(InputForm)) ]$List(InputForm) lambdaFuncallSpad(compile(function(e, declare DI, x), t))$Lisp
spad
Compiling FriCAS source code from file
/var/lib/zope2.10/instance/axiom-wiki/var/LatexWiki/1657093691659194492-25px002.spad
using old system compiler.
Value = |lambdaFuncallSpad|
MKULF abbreviates package MakeUnaryLispFunction
------------------------------------------------------------------------
initializing NRLIB MKULF for MakeUnaryLispFunction
compiling into NRLIB MKULF
importing MakeFunction S
compiling exported compiledFunction : (S,
(time taken in buildFunctor: 0)
;;; *** |MakeUnaryLispFunction| REDEFINED
;;; *** |MakeUnaryLispFunction| REDEFINED
Time: 0 SEC.
Cumulative Statistics for Constructor MakeUnaryLispFunction
Time: 0 seconds
finalizing NRLIB MKULF
Processing MakeUnaryLispFunction for Browser database:
--------constructor---------
--------(compiledFunction ((Symbol) S (Symbol)))---------
; compiling file "/var/aw/var/LatexWiki/MKULF.NRLIB/MKULF.lsp" (written 04 APR 2022 07:16:37 PM):
; /var/aw/var/LatexWiki/MKULF.NRLIB/MKULF.fasl written
; compilation finished in 0:00:00.014
------------------------------------------------------------------------
MakeUnaryLispFunction is now explicitly exposed in frame initial
MakeUnaryLispFunction will be automatically loaded when needed from
/var/aw/var/LatexWiki/MKULF.NRLIB/MKULFThen we can use this in Axiom...
fricas
-- define Newton in Lisp
fricas
)lisp (defun newton (f dfdx x0 &optional (tol 1.0d-10) (maxit 100)) (let ((xt x0) (fxt (funcall f x0)) (iternum 0)) (loop (setf xt (- xt (/ fxt (funcall dfdx xt)))) (setf fxt (funcall f xt)) (if (< (abs fxt) tol) (return xt)) (incf iternum) (if (>= iternum maxit) (error "Maximum iterations exceeded.")))))
Value = NEWTON -- Spad->Lisp function translation compiledDF(expr: EXPR FLOAT,x: Symbol):Symbol == compiledFunction(expr, x)$MakeUnaryLispFunction(EXPR FLOAT, DFLOAT, DFLOAT)
Function declaration compiledDF : (Expression(Float),Symbol) -> Symbol has been added to workspace.
Type: Void
fricas
-- a short function to call the Lisp code newtonUsingLisp(f:Expression Float,x:Symbol, x0:DFLOAT):DFLOAT == float(NEWTON(compiledDF(f, x), compiledDF(D(f, x), x), x0)$Lisp)
Function declaration newtonUsingLisp : (Expression(Float),Symbol, DoubleFloat) -> DoubleFloat has been added to workspace.
Type: Void
fricas
newtonUsingLisp(x**2-2.0,x, 2.0::SF)-sqrt(2.0::SF)
There are no library operations named ** Use HyperDoc Browse or issue )what op ** to learn if there is any operation containing " ** " in its name.
Cannot find a definition or applicable library operation named ** with argument type(s) Expression(Float) PositiveInteger
Perhaps you should use "@" to indicate the required return type,or "$" to specify which version of the function you need.
Third, we can call a compiled Spad function in Lisp. Defining some example functions in Spad:
spad
)abbrev package TESTP TestPackage
R ==> DoubleFloat
TestPackage: with
f:R -> R
dfdx:R -> R
== add
f(x) == x*x - 2.0::R
dfdx(x) == 2*x
spad
Compiling FriCAS source code from file
/var/lib/zope2.10/instance/axiom-wiki/var/LatexWiki/2983317606354902709-25px004.spad
using old system compiler.
TESTP abbreviates package TestPackage
------------------------------------------------------------------------
initializing NRLIB TESTP for TestPackage
compiling into NRLIB TESTP
compiling exported f : DoubleFloat -> DoubleFloat
Time: 0.01 SEC.
compiling exported dfdx : DoubleFloat -> DoubleFloat
Time: 0.01 SEC.
(time taken in buildFunctor: 0)
;;; *** |TestPackage| REDEFINED
;;; *** |TestPackage| REDEFINED
Time: 0 SEC.
Cumulative Statistics for Constructor TestPackage
Time: 0.02 seconds
finalizing NRLIB TESTP
Processing TestPackage for Browser database:
--->-->TestPackage(constructor): Not documented!!!!
--->-->TestPackage((f ((DoubleFloat) (DoubleFloat)))): Not documented!!!!
--->-->TestPackage((dfdx ((DoubleFloat) (DoubleFloat)))): Not documented!!!!
--->-->TestPackage(): Missing Description
; compiling file "/var/aw/var/LatexWiki/TESTP.NRLIB/TESTP.lsp" (written 04 APR 2022 07:16:37 PM):
; /var/aw/var/LatexWiki/TESTP.NRLIB/TESTP.fasl written
; compilation finished in 0:00:00.009
------------------------------------------------------------------------
TestPackage is now explicitly exposed in frame initial
TestPackage will be automatically loaded when needed from
/var/aw/var/LatexWiki/TESTP.NRLIB/TESTPand then calling those functions in Lisp from Axiom:
fricas
)lisp (defun |lispFunctionFromSpad| (f dom args) (let ((spadf (|getFunctionFromDomain| f (list dom) args))) (lambda (x) (spadcall x spadf))))
Value = |lispFunctionFromSpad| float(NEWTON(lispFunctionFromSpad(f,'TestPackage, ['DoubleFloat])$Lisp, lispFunctionFromSpad(dfdx, 'TestPackage, ['DoubleFloat])$Lisp, 2::SF)$Lisp)-sqrt(2.0::SF)
 | (1) |
Type: DoubleFloat?
Finally, within a Spad package, it is even easier to call a function through Lisp:
spad
-- the following approach is due to Waldek Hebisch
)lisp (defun |mkLispFunction1| (f) (lambda (x) (SPADCALL x f)))
)abbrev package TESTP TestPackage
R ==> DoubleFloat
I ==> Integer
TestPackage: with
newton:(R->R,
spad
Compiling FriCAS source code from file
/var/lib/zope2.10/instance/axiom-wiki/var/LatexWiki/5566883931090623725-25px006.spad
using old system compiler.
Value = |mkLispFunction1|
TESTP abbreviates package TestPackage
------------------------------------------------------------------------
initializing NRLIB TESTP for TestPackage
compiling into NRLIB TESTP
compiling exported newton : (DoubleFloat -> DoubleFloat,
compiling exported newton : (DoubleFloat -> DoubleFloat,
(time taken in buildFunctor: 0)
;;; *** |TestPackage| REDEFINED
;;; *** |TestPackage| REDEFINED
Time: 0 SEC.
Cumulative Statistics for Constructor TestPackage
Time: 0.01 seconds
finalizing NRLIB TESTP
Processing TestPackage for Browser database:
--->/var/lib/zope2.10/instance/axiom-wiki/var/LatexWiki/2983317606354902709-25px004.spad-->TestPackage(constructor): Not documented!!!!
--->/var/lib/zope2.10/instance/axiom-wiki/var/LatexWiki/2983317606354902709-25px004.spad-->TestPackage((newton ((DoubleFloat) (Mapping (DoubleFloat) (DoubleFloat)) (Mapping (DoubleFloat) (DoubleFloat)) (DoubleFloat) (DoubleFloat) (Integer)))): Not documented!!!!
--->/var/lib/zope2.10/instance/axiom-wiki/var/LatexWiki/2983317606354902709-25px004.spad-->TestPackage((newton ((DoubleFloat) (Mapping (DoubleFloat) (DoubleFloat)) (Mapping (DoubleFloat) (DoubleFloat)) (DoubleFloat)))): Not documented!!!!
--->/var/lib/zope2.10/instance/axiom-wiki/var/LatexWiki/2983317606354902709-25px004.spad-->TestPackage(): Missing Description
; compiling file "/var/aw/var/LatexWiki/TESTP.NRLIB/TESTP.lsp" (written 04 APR 2022 07:16:37 PM):
; /var/aw/var/LatexWiki/TESTP.NRLIB/TESTP.fasl written
; compilation finished in 0:00:00.009
------------------------------------------------------------------------
TestPackage is already explicitly exposed in frame initial
TestPackage will be automatically loaded when needed from
/var/aw/var/LatexWiki/TESTP.NRLIB/TESTPthen
fricas
R ==> DoubleFloat
Type: Void
fricas
newton((x:R):R +-> x**2-2.0::R,(x:R):R +-> 2.0::R*x, 3.0::SF, 1.0e-15, 100)$TestPackage-sqrt(2.0)
There are no library operations named ** Use HyperDoc Browse or issue )what op ** to learn if there is any operation containing " ** " in its name. Cannot find a definition or applicable library operation named ** with argument type(s) DoubleFloat PositiveInteger
Perhaps you should use "@" to indicate the required return type,or "$" to specify which version of the function you need. FriCAS will attempt to step through and interpret the code. There are no library operations named ** Use HyperDoc Browse or issue )what op ** to learn if there is any operation containing " ** " in its name.
Cannot find a definition or applicable library operation named ** with argument type(s) Symbol PositiveInteger
Perhaps you should use "@" to indicate the required return type,or "$" to specify which version of the function you need.