src/airy.lisp

   1 ;;; Airy functions Ai(z) and Bi(z) - A&S 10.4
   2 ;;;
   3 ;;; airy_ai(z)   - Airy function Ai(z)
   4 ;;; airy_dai(z)  - Derivative of Airy function Ai(z)
   5 ;;; airy_bi(z)   - Airy function Bi(z)
   6 ;;; airy_dbi(z)  - Derivative of Airy function Bi(z)
   7
   8 ;;;; Copyright (C) 2005 David Billinghurst
   9
  10 ;;;; airy.lisp is free software; you can redistribute it
  11 ;;;; and/or modify it under the terms of the GNU General Public
  12 ;;;; License as published by the Free Software Foundation; either
  13 ;;;; version 2, or (at your option) any later version.
  14
  15 ;;;; airy.lisp is distributed in the hope that it will be
  16 ;;;; useful, but WITHOUT ANY WARRANTY; without even the implied
  17 ;;;; warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
  18 ;;;; See the GNU General Public License for more details.
  19
  20 ;;;; You should have received a copy of the GNU General Public License
  21 ;;;; along with command-line.lisp; see the file COPYING.  If not,
  22 ;;;; write to the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA
  23
  24 (in-package :maxima)
  25
  26 (declaim (special *flonum-op*))
  27
  28 ;; Airy Ai function
  29 (defmfun $airy_ai (z)
  30   "Airy function Ai(z)"
  31   (simplify (list '(%airy_ai) (resimplify z))))
  32
  33 (defprop %airy_ai simplim%airy_ai simplim%function)
  34 (defprop %airy_ai ((z) ((%airy_dai) z)) grad)
  35
  36 ;; airy_ai distributes over lists, matrices, and equations
  37 (defprop %airy_ai (mlist $matrix mequal) distribute_over)
  38
  39 ;; airy_ai has mirror symmetry
  40 (defprop %airy_ai t commutes-with-conjugate)
  41
  42 ;; Integral of Ai(z)
  43 ;; http://functions.wolfram.com/03.05.21.0002.01
  44 ;; (z/(3^(2/3)*gamma(2/3)))*hypergeometric([1/3],[2/3,4/3],z^3/9)
  45 ;; - (3^(1/6)/(4*%pi))*z^2*gamma(2/3)*hypergeometric([2/3],[4/3,5/3],z^3/9);
  46 (defprop %airy_ai
  47   ((z)
  48    ((mplus)
  49     ((mtimes)
  50      ((mexpt) 3 ((rat) -2 3))
  51      ((mexpt) ((%gamma) ((rat) 2 3)) -1)
  52      (($hypergeometric)
  53       ((mlist) ((rat) 1 3))
  54       ((mlist) ((rat) 2 3) ((rat) 4 3))
  55       ((mtimes) ((rat) 1 9) ((mexpt) z 3)))
  56      z)
  57    ((mtimes)
  58     ((rat) -1 4) ((mexpt) 3 ((rat) 1 6)) ((mexpt) $%pi -1) ((%gamma) ((rat) 2 3))
  59     (($hypergeometric)
  60      ((mlist) ((rat) 2 3))
  61      ((mlist) ((rat)  4 3) ((rat) 5 3))
  62      ((mtimes) ((rat) 1 9) ((mexpt) z 3)))
  63     ((mexpt) z 2))))
  64   integral)
  65
  66 (defun airy-ai (z)
  67   (cond ((floatp z) (airy-ai-real z))
  68         ((complexp z) (airy-ai-complex z))))
  69
  70 (setf (gethash '%airy_ai *flonum-op*) #'airy-ai)
  71
  72 (defun simplim%airy_ai (expr var val)
  73   ; Look for the limit of the argument
  74   (let ((z (limit (cadr expr) var val 'think)))
  75     (cond ((or (eq z '$inf)   ; A&S 10.4.59
  76                (eq z '$minf)) ; A&S 10.4.60
  77            0)
  78           (t
  79            ; Handle other cases with the function simplifier
  80            (simplify (list '(%airy_ai) z))))))
  81
  82 (def-simplifier airy_ai (z)
  83   (cond ((equal z 0)         ; A&S 10.4.4: Ai(0) = 3^(-2/3)/gamma(2/3)
  84          '((mtimes simp)
  85            ((mexpt simp) 3 ((rat simp) -2 3))
  86            ((mexpt simp) ((%gamma simp) ((rat simp) 2 3)) -1)))
  87         ((flonum-eval (mop form) z))
  88         (t (give-up))))
  89
  90
  91 ;; Derivative dAi/dz of Airy function Ai(z)
  92 (defmfun $airy_dai (z)
  93   "Derivative dAi/dz of Airy function Ai(z)"
  94   (simplify (list '(%airy_dai) (resimplify z))))
  95
  96 (defprop %airy_dai simplim%airy_dai simplim%function)
  97 (defprop %airy_dai ((z) ((mtimes) z ((%airy_ai) z))) grad)
  98 (defprop %airy_dai ((z) ((%airy_ai) z)) integral)
  99
 100 ;; airy_dai distributes over lists, matrices, and equations
 101 (defprop %airy_dai (mlist $matrix mequal) distribute_over)
 102
 103 ;; airy_dai has mirror symmetry
 104 (defprop %airy_dai t commutes-with-conjugate)
 105
 106 (defun airy-dai (z)
 107   (cond ((floatp z) (airy-dai-real z))
 108         ((complexp z) (airy-dai-complex z))))
 109
 110 (setf (gethash '%airy_dai *flonum-op*) #'airy-dai)
 111
 112 (defun simplim%airy_dai (expr var val)
 113   ; Look for the limit of the argument
 114   (let ((z (limit (cadr expr) var val 'think)))
 115     (cond ((eq z '$inf) ; A&S 10.4.61
 116            0)
 117           ((eq z '$minf) ; A&S 10.4.62
 118            '$und)
 119           (t
 120            ; Handle other cases with the function simplifier
 121            (simplify (list '(%airy_dai) z))))))
 122
 123 (def-simplifier airy_dai (z)
 124   (cond ((equal z 0)       ; A&S 10.4.5: Ai'(0) = -3^(-1/3)/gamma(1/3)
 125          '((mtimes simp) -1
 126            ((mexpt simp) 3 ((rat simp) -1 3))
 127            ((mexpt simp) ((%gamma simp) ((rat simp) 1 3)) -1)))
 128         ((flonum-eval (mop form) z))
 129         (t (give-up))))
 130
 131 ;; Airy Bi function
 132 (defmfun $airy_bi (z)
 133   "Airy function Bi(z)"
 134   (simplify (list '(%airy_bi) (resimplify z))))
 135
 136 (defprop %airy_bi simplim%airy_bi simplim%function)
 137 (defprop %airy_bi ((z) ((%airy_dbi) z)) grad)
 138
 139 ;; airy_bi distributes over lists, matrices, and equations
 140 (defprop %airy_bi (mlist $matrix mequal) distribute_over)
 141
 142 ;; airy_bi has mirror symmetry
 143 (defprop %airy_bi t commutes-with-conjugate)
 144
 145 ;; Integral of Bi(z)
 146 ;; http://functions.wolfram.com/03.06.21.0002.01
 147 ;; (z/(3^(1/6)*gamma(2/3)))*hypergeometric([1/3],[2/3,4/3],z^3/9)
 148 ;; + (3^(2/3)/(4*%pi))*z^2*gamma(2/3)*hypergeometric([2/3],[4/3,5/3],z^3/9);
 149 (defprop %airy_bi
 150   ((z)
 151    ((mplus)
 152     ((mtimes)
 153      ((mexpt) 3 ((rat) -1 6))
 154      ((mexpt) ((%gamma) ((rat) 2 3)) -1)
 155      (($hypergeometric)
 156       ((mlist) ((rat) 1 3))
 157       ((mlist) ((rat) 2 3) ((rat) 4 3))
 158       ((mtimes) ((rat) 1 9) ((mexpt) z 3)))
 159      z)
 160    ((mtimes)
 161     ((rat) 1 4) ((mexpt) 3 ((rat) 2 3)) ((mexpt) $%pi -1) ((%gamma) ((rat) 2 3))
 162     (($hypergeometric)
 163      ((mlist) ((rat) 2 3))
 164      ((mlist) ((rat)  4 3) ((rat) 5 3))
 165      ((mtimes) ((rat) 1 9) ((mexpt) z 3)))
 166     ((mexpt) z 2))))
 167   integral)
 168
 169 (defun airy-bi (z)
 170   (cond ((floatp z) (airy-bi-real z))
 171         ((complexp z) (airy-bi-complex z))))
 172
 173 (setf (gethash '%airy_bi *flonum-op*) #'airy-bi)
 174
 175 (defun simplim%airy_bi (expr var val)
 176   ; Look for the limit of the argument
 177   (let ((z (limit (cadr expr) var val 'think)))
 178     (cond ((eq z '$inf) ; A&S 10.4.63
 179            '$inf)
 180           ((eq z '$minf) ; A&S 10.4.64
 181            0)
 182           (t
 183            ; Handle other cases with the function simplifier
 184            (simplify (list '(%airy_bi) z))))))
 185
 186 (def-simplifier airy_bi (z)
 187   (cond ((equal z 0) ; A&S 10.4.4: Bi(0) = sqrt(3) 3^(-2/3)/gamma(2/3)
 188          '((mtimes simp)
 189            ((mexpt simp) 3 ((rat simp) -1 6))
 190            ((mexpt simp) ((%gamma simp) ((rat simp) 2 3)) -1)))
 191         ((flonum-eval (mop form) z))
 192         (t (give-up))))
 193
 194
 195 ;; Derivative dBi/dz of Airy function Bi(z)
 196 (defmfun $airy_dbi (z)
 197   "Derivative dBi/dz of Airy function Bi(z)"
 198   (simplify (list '(%airy_dbi) (resimplify z))))
 199
 200 (defprop %airy_dbi simplim%airy_dbi simplim%function)
 201 (defprop %airy_dbi ((z) ((mtimes) z ((%airy_bi) z))) grad)
 202 (defprop %airy_dbi ((z) ((%airy_bi) z)) integral)
 203
 204 ;; airy_dbi distributes over lists, matrices, and equations
 205 (defprop %airy_dbi (mlist $matrix mequal) distribute_over)
 206
 207 ;; airy_dbi has mirror symmetry
 208 (defprop %airy_dbi t commutes-with-conjugate)
 209
 210 (defun airy-dbi (z)
 211   (cond ((floatp z) (airy-dbi-real z))
 212         ((complexp z) (airy-dbi-complex z))))
 213
 214 (setf (gethash '%airy_dbi *flonum-op*) #'airy-dbi)
 215
 216 (defun simplim%airy_dbi (expr var val)
 217   ; Look for the limit of the argument
 218   (let ((z (limit (cadr expr) var val 'think)))
 219     (cond ((eq z '$inf) ; A&S 10.4.66
 220            '$inf)
 221           ((eq z '$minf) ; A&S 10.4.67
 222            '$und)
 223           (t
 224            ; Handle other cases with the function simplifier
 225            (simplify (list '(%airy_dbi) z))))))
 226
 227 (def-simplifier airy_dbi (z)
 228   (cond ((equal z 0) ; A&S 10.4.5: Bi'(0) = sqrt(3) 3^(-1/3)/gamma(1/3)
 229          '((mtimes simp)
 230            ((mexpt simp) 3 ((rat simp) 1 6))
 231            ((mexpt simp) ((%gamma simp) ((rat simp) 1 3)) -1)))
 232         ((flonum-eval (mop form) z))
 233         (t (give-up))))
 234
 235 ;; Numerical routines using slatec functions
 236
 237 (defun airy-ai-real (z)
 238   " Airy function Ai(z) for real z"
 239   (declare (type flonum z))
 240   ;; slatec:dai issues underflow warning for z > zmax.  See dai.{f,lisp}
 241   ;; This value is correct for IEEE double precision
 242   (let ((zmax 92.5747007268))
 243     (declare (type flonum zmax))
 244     (if (< z zmax) (slatec:dai z) 0.0)))
 245
 246 (defun airy-ai-complex (z)
 247   "Airy function Ai(z) for complex z"
 248   (declare (type (complex flonum) z))
 249   (multiple-value-bind (var-0 var-1 var-2 var-3 air aii nz ierr)
 250       (slatec:zairy (realpart z) (imagpart z) 0 1 0.0 0.0 0 0)
 251     (declare (type flonum air aii)
 252              (type f2cl-lib:integer4 nz ierr)
 253              (ignore var-0 var-1 var-2 var-3))
 254     ;; Check nz and ierr for errors
 255     (if (and (= nz 0) (= ierr 0)) (complex air aii) nil)))
 256
 257 (defun airy-dai-real (z)
 258   "Derivative dAi/dz of Airy function Ai(z) for real z"
 259   (declare (type flonum z))
 260   (let ((rz (sqrt (abs z)))
 261         (c (* 2/3 (expt (abs z) 3/2))))
 262     (declare (type flonum rz c))
 263     (multiple-value-bind (var-0 var-1 var-2 ai dai)
 264         (slatec:djairy z rz c 0.0 0.0)
 265       (declare (ignore var-0 var-1 var-2 ai))
 266       dai)))
 267
 268 (defun airy-dai-complex (z)
 269   "Derivative dAi/dz of Airy function Ai(z) for complex z"
 270   (declare (type (complex flonum) z))
 271   (multiple-value-bind (var-0 var-1 var-2 var-3 air aii nz ierr)
 272       (slatec:zairy (realpart z) (imagpart z) 1 1 0.0 0.0 0 0)
 273     (declare (type flonum air aii)
 274              (type f2cl-lib:integer4 nz ierr)
 275              (ignore var-0 var-1 var-2 var-3))
 276     ;; Check nz and ierr for errors
 277     (if (and (= nz 0) (= ierr 0)) (complex air aii) nil)))
 278
 279 (defun airy-bi-real (z)
 280   "Airy function Bi(z) for real z"
 281   (declare (type flonum z))
 282   ;; slatec:dbi issues overflows for z > zmax.  See dbi.{f,lisp}
 283   ;; This value is correct for IEEE double precision
 284   (let ((zmax 104.2179765192136))
 285     (declare (type flonum zmax))
 286     (if (< z zmax) (slatec:dbi z) nil)))
 287
 288 (defun airy-bi-complex (z)
 289   "Airy function Bi(z) for complex z"
 290   (declare (type (complex flonum) z))
 291   (multiple-value-bind (var-0 var-1 var-2 var-3 bir bii ierr)
 292       (slatec:zbiry (realpart z) (imagpart z) 0 1 0.0 0.0 0)
 293     (declare (type flonum bir bii)
 294              (type f2cl-lib:integer4 ierr)
 295              (ignore var-0 var-1 var-2 var-3))
 296     ;; Check ierr for errors
 297     (if (= ierr 0) (complex bir bii) nil)))
 298
 299 (defun airy-dbi-real (z)
 300   "Derivative dBi/dz of Airy function Bi(z) for real z"
 301   (declare (type flonum z))
 302   ;; Overflows for z > zmax.
 303   ;; This value is correct for IEEE double precision
 304   (let ((zmax 104.1525))
 305     (declare (type flonum zmax))
 306     (if (< z zmax)
 307         (let ((rz (sqrt (abs z)))
 308               (c (* 2/3 (expt (abs z) 3/2))))
 309         (declare (type flonum rz c))
 310         (multiple-value-bind (var-0 var-1 var-2 bi dbi)
 311             (slatec:dyairy z rz c 0.0 0.0)
 312           (declare (type flonum bi dbi)
 313                    (ignore var-0 var-1 var-2 bi))
 314           dbi))
 315       ;; Will overflow.  Return unevaluated.
 316       nil)))
 317
 318 (defun airy-dbi-complex (z)
 319   "Derivative dBi/dz of Airy function Bi(z) for complex z"
 320   (declare (type (complex flonum) z))
 321   (multiple-value-bind (var-0 var-1 var-2 var-3 bir bii ierr)
 322       (slatec:zbiry (realpart z) (imagpart z) 1 1 0.0 0.0 0)
 323     (declare (type flonum bir bii)
 324              (type f2cl-lib:integer4 ierr)
 325              (ignore var-0 var-1 var-2 var-3))
 326     ;; Check ierr for errors
 327     (if (= ierr 0) (complex bir bii) nil)))