quadde/quadde.lisp at main · incfbelgiannode/quadde · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
;;; Maxima interface to double exponential integration
;;
; (c) Dimiter Prodanov
; version 1.1.2 17 Sept 2023
; - new condition in quadde_argument_check
; version 1.1.1 6 Oct 2022
; - change in quad_intdeo
; version 1.1 13 March 2021
; - change in condition checking
; - compatibilty with Maxima 44
; version 1.0 23 Jan 2021
; updates in quadde_argument_check
; input from Michel Talon
;
; * M.Mori, Developments in the double exponential formula for numerical integration, Proceedings of the International Congress of Mathematicians, Kyoto 1990, 1991, Springer-Verlag, 1585-1594.
; * H.Takahasi and M.Mori, Double exponential formulas for numerical integration, Pub. RIMS Kyoto Univ. 9, 1974, 721-741
; * T.Ooura and M.Mori, The Double exponential formula for oscillatory functions over the half infinite interval, Journal of Computational and Applied Mathematics 38, 1991, 353-360
; * M.Mori and T.Ooura, Double exponential formulas for Fourier type integrals with a divergent integrand, Contributions in Numerical Mathematics, ed. R.P.Agarwal, World Scientific Series in Applicable Analysis, 2, 1993, 301-308
; * H.Toda and H.Ono, Some remarks for efficient usage of the double exponential formulas(in Japanese), Kokyuroku, RIMS, Kyoto Univ. 339, 1978, 74-109.
;
; 2021

(in-package :maxima)

(defvar *debug-quadde*
  NIL
  "Set to non-NIL to enable printing of the error object when the
  intde routines throw an error.")

(defun coerce-float-d (x)
	 (coerce (maxima::$float (maxima::meval* x)) 'double-float)
)

;; error checking similar to that done by $defint
;; we don't really need to simplify
(defun quadde_argument_check (expr var ll ul)
  (cond (($constantp var)
	 (merror "Variable of integration not a variable: ~M"
		 var)))
  (cond (($subvarp var)
	 (merror "Improper variable of integration: ~M"
		 var)))
  (cond ( (not (atom var))
	 (merror "Variable of integration: ~M  not an atom"
		 var)))
  (cond ( (and  ($freeof var expr) (not (= 0.0 expr ) ) )
	 (merror "Variable ~M not in ~M"
		 var expr)))
   (cond ((or (among var ul) (among var ll))
	  (merror "Terminal contains variable of integration: ~M" var)))
   (cond ((not (and ($numberp (coerce-float-d ul)) ($numberp (coerce-float-d ll))))
	  (merror "Terminal not a number ~M , ~M" ll ul)))
)


;; wrapper around intde
(defmfun $quad_intde (fun var a b &key
		 (epsrel 1e-8)
		 (epsabs 0.0))
  (quadde_argument_check fun var a b)
  (let
  ( (f (coerce-float-fun fun)))
	(handler-case
	(flet ((ff ( x)
          (funcall f x)))
	(multiple-value-bind ( i err )
		(quadde::intde  #'ff
			 (float-or-lose a)
			 (float-or-lose b)
			 (float-or-lose epsrel)
			 (float-or-lose epsabs)
		)
	  (list '(mlist) i err))
	  )
	   (error (e)
	 (when *debug-quadde*
	   (format t "~S" e))
	 `(($quad_intde) ,fun ,var ,a ,b
	   ((mequal) $epsrel ,epsrel)
	   ((mequal) $epsabs ,epsabs)
	  )
	  ))
	 )
)

;; wrapper around intdei
(defmfun $quad_intdei (fun var a  &key
		 (epsrel 1e-8)
		 (epsabs 0.0))
  (quadde_argument_check fun var a a)
  (let
  ( (f (coerce-float-fun fun)))
	(handler-case
	(flet ((ff ( x)
          (funcall f x)))
	(multiple-value-bind ( i err )
		(quadde::intdei  #'ff
			 (float-or-lose a)
			 (float-or-lose epsrel)
			 (float-or-lose epsabs)
		)
	  (list '(mlist) i err))
	  )
	   (error (e)
	 (when *debug-quadde*
	   (format t "~S" e))
	 `(($quad_intdei) ,fun ,var ,a
	   ((mequal) $epsrel ,epsrel)
	   ((mequal) $epsabs ,epsabs)
	  )
	  ))
	 )
)

;; wrapper around intdeo
(defmfun $quad_intdeo (fun var a omega &key
		 (epsrel 1e-8)
		 (epsabs 0.0))
  (quadde_argument_check fun var a omega)
  (let
  ( (f (coerce-float-fun fun)))
	(handler-case
	(flet ((ff ( x)
          (funcall f x)))
	(multiple-value-bind ( i err )
		(quadde::intdeo  #'ff
			 (float-or-lose a)
			 (coerce-float-d omega)
			 (float-or-lose epsrel)
			 (float-or-lose epsabs)
		)
	  (list '(mlist) i err))
	  )
	   (error (e)
	 (when *debug-quadde*
	   (format t "~S" e))
	 `(($quad_intdeo) ,fun ,var ,a ,omega
	   ((mequal) $epsrel ,epsrel)
	   ((mequal) $epsabs ,epsabs)
	  )
	  ))
	 )
)