share/numeric/romberg.lisp

   1 ;;; -*-  Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10 -*- ;;;;
   2 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
   3 ;;;     The data in this file contains enhancments.                    ;;;;;
   4 ;;;                                                                    ;;;;;
   5 ;;;  Copyright (c) 1984,1987 by William Schelter,University of Texas   ;;;;;
   6 ;;;     All rights reserved                                            ;;;;;
   7 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
   8 ;;;     (c) Copyright 1982 Massachusetts Institute of Technology         ;;;
   9 ;;; Original code by CFFK.  Modified to interface correctly with TRANSL  ;;;
  10 ;;; and the rest of macsyma by GJC                                       ;;;
  11 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  12
  13 (in-package :maxima)
  14
  15 (macsyma-module rombrg)
  16
  17 (load-macsyma-macros transm)
  18
  19 (defmvar  $rombergit 11 "the maximum number of iterations")
  20 (defmvar  $rombergmin 0 "the minimum number of iterations")
  21 (defmvar  $rombergtol 1e-4 "the relative tolerance of error")
  22 (defmvar  $rombergabs 0.0 "the absolute tolerance of error")
  23 (defmvar  $rombergit_used 0 "the number of iterations actually used.")
  24
  25 (defun romberg-subr (f left right)
  26    (let (a b
  27          (x 0.0)
  28          (tt (make-array  $rombergit :element-type 'flonum))
  29          (rr  (make-array  $rombergit :element-type 'flonum)))
  30      (let (($numer t) ($%enumer t))
  31        (setq
  32          a ($float left)
  33          b ($float right)))
  34      (if (not (and (numberp a) (numberp b)))
  35        (return-from romberg-subr (values nil a b)))
  36      (setq x (- b a))
  37      (let
  38        ((fa (funcall f a))
  39         (fb (funcall f b))
  40         (fab (funcall f (* (+ b a) 0.5))))
  41        (if
  42          (or
  43            (not (numberp fa))
  44            (not (numberp fb))
  45            (not (numberp fab)))
  46          (return-from romberg-subr (values nil a b))
  47          (progn
  48            (setf (aref tt 0) (* x (+ fb fa) 0.5))
  49            (setf (aref rr 0) (* x fab)))))
  50
  51      (do ((l 1 (1+ l))
  52           (m 4 (* m 2))
  53           (y 0.0)
  54           (z 0.0)
  55           (cerr 0.0))
  56          ((= l $rombergit)
  57           (merror "`romberg' failed to converge"))
  58        (declare (fixnum l m))
  59        (setq y (float m) z (/ x y))
  60        (setf (aref tt l) (* (+ (aref tt (1- l))
  61                                (aref rr (1- l))) 0.5))
  62        (setf (aref rr l) 0.0)
  63        (do ((i 1 (+ i 2)))
  64            ((> i m))
  65        (let
  66          ((fzi+a (funcall f (+ (* z (float i)) a))))
  67          (if (not (numberp fzi+a))
  68            (return-from romberg-subr (values nil a b))
  69            (setf (aref rr l) (+ fzi+a (aref rr l))))))
  70
  71        (setf (aref rr l) (* z (aref rr l) 2))
  72        (setq y 0.0)
  73        (do ((k l (1- k))) ((= k 0))
  74          (declare (fixnum k))
  75          (setq y (+ (* y 4) 3))
  76          (setf (aref  tt (1- k))
  77                (+ (/ (- (aref tt k)
  78                         (aref tt (1- k))) y)
  79                   (aref tt k)))
  80          (setf (aref rr (1- k))
  81                (+ (/ (- (aref rr k) (aref rr (1- k))) y)
  82                   (aref rr k))))
  83        (setq y (* (+ (aref tt 0) (aref rr 0)) 0.5))
  84 ;;; this is the WIN condition test.
  85        (cond ((and (or (not (< $rombergabs
  86                                (setq cerr (abs (- (aref tt 0) (aref rr 0))))))
  87                        (not (< $rombergtol
  88                                ;; cerr = "calculated error"; used for ^]
  89                                (setq cerr (/ cerr (if (zerop y) 1.0 (abs y)))))))
  90                    (> l $rombergmin))
  91               (setq $rombergit_used l)
  92               (return-from romberg-subr (values y nil nil)))))))
  93
  94 (defun $romberg (&rest args)
  95   (cond
  96     ((= (length args) 3)
  97      (multiple-value-bind
  98        (result left right)
  99        ;; BIND EVIL SPECIAL VARIABLE *PLOT-REALPART* HERE ...
 100        (let (($numer t) ($%enumer t) (*plot-realpart* nil))
 101          (romberg-subr
 102            (coerce-float-fun (first args))
 103            (second args)
 104            (third args)))
 105        (if (numberp result)
 106          result
 107          `(($romberg) ,(first args) ,left ,right))))
 108     ((= (length args) 4)
 109      (when ($constantp (second args))
 110        (merror
 111          "romberg: variable of integration cannot be a constant; found ~M"
 112          (second args)))
 113      (multiple-value-bind
 114        (result left right)
 115        ;; BIND EVIL SPECIAL VARIABLE *PLOT-REALPART* HERE ...
 116        (let (($numer t) ($%enumer t) (*plot-realpart* nil))
 117          (romberg-subr
 118            (coerce-float-fun (first args) `((mlist) ,(second args)))
 119            (third args)
 120            (fourth args)))
 121        (if (numberp result)
 122          result
 123          `(($romberg) ,(first args) ,(second args) ,left ,right))))
 124     (t
 125       (wna-err '$romberg))))