← Back to branch summary

~ubuntu-branches/debian/squeeze/maxima/squeeze

« back to all changes in this revision

Viewing changes to src/defint.lisp

  • Committer: Bazaar Package Importer
  • Author(s): Camm Maguire
  • Date: 2010-04-30 13:30:33 UTC
  • mto: This revision was merged to the branch mainline in revision 12.
  • Revision ID: james.westby@ubuntu.com-20100430133033-wtewap0zdnmsix1y
Tags: upstream-5.21.1
Import upstream version 5.21.1

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/defint.lisp
779
779
  (setq *global-defint-assumptions*
780
780
        (cons (assume '((mlessp) epsilon 1.0e-8))
781
781
              *global-defint-assumptions*))
782
 
  ;; EPSILON is used in principal vaule code to denote the familiar
 
782
  ;; EPSILON is used in principal value code to denote the familiar
783
783
  ;; mathematical entity.
784
784
  (setq *global-defint-assumptions*
785
785
        (cons (assume '((mgreaterp) prin-inf 1.0e+8))
999
999
         (format t "Principal Value~%")
1000
1000
         (setq pcprntd t))))
1001
1001
 
 
1002
;; e is of form poly(x)*exp(m*%i*x)
 
1003
;; s is degree of denominator
 
1004
;; adds e to bptu or bptd according to sign of m
1002
1005
(defun rib (e s)
1003
1006
  (let (*updn c)
1004
1007
    (cond ((or (mnump e) (constant e))
1013
1016
                      (cond (*updn (setq bptu (cons e bptu)))
1014
1017
                            (t (setq bptd (cons e bptd))))))))))
1015
1018
 
 
1019
;; check term is of form poly(x)*exp(m*%i*x)
 
1020
;; n is degree of denominator
1016
1021
(defun ptimes%e (term n)
1017
 
  (cond ((polyinx term var nil) term)
1018
 
        ((and (mexptp term)
 
1022
  (cond ((and (mexptp term)
1019
1023
              (eq (cadr term) '$%e)
1020
1024
              (polyinx (caddr term) var nil)
1021
1025
              (eq ($sign (m+ (deg ($realpart (caddr term))) -1))
1065
1069
                                           (cadr n))
1066
1070
                                          (t 0.))))))))
1067
1071
 
1068
 
 
 
1072
;; exponentialize sin and cos
1069
1073
(defun sconvert (e)
1070
1074
  (cond ((atom e) e)
1071
1075
        ((polyinx e var nil) e)
1107
1111
           (t (return (simplify (cons (list (caar e))
1108
1112
                                      (mapcar #'esap (cdr e)))))))))
1109
1113
 
 
1114
;; computes integral from minf to inf for expressions of the form
 
1115
;; exp(%i*m*x)*r(x) by residues on either the upper half
 
1116
;;                plane or the lower half plane, depending on whether
 
1117
;;                m is positive or negative.  [wang p. 77]
 
1118
;;
 
1119
;; exponentializes sin and cos before applying residue method.
 
1120
;; can handle some expressions with poles on real line, such as
 
1121
;; sin(x)*cos(x)/x.
1110
1122
(defun mtosc (grand)
1111
1123
  (numden grand)
1112
1124
  (let ((n nn*)
1113
1125
        (d dn*)
 
1126
        ratterms ratans
1114
1127
        plf bptu bptd s upans downans)
1115
1128
    (cond ((not (or (polyinx d var nil)
1116
1129
                    (and (setq grand (%einvolve d))
1127
1140
                                       (t (sconvert n)))))
1128
1141
                (cond ((mplusp n)  (setq n (cdr n)))
1129
1142
                      (t (setq n (list n))))
1130
 
                (do ((do-var n (cdr do-var)))
1131
 
                    ((null do-var) t)
1132
 
                  (cond ((rib (car do-var) s))
 
1143
                (dolist (term n t)
 
1144
                  (cond ((polyinx term var nil)
 
1145
                         ;; call to $expand can create rational terms
 
1146
                         ;; with no exp(m*%i*x)
 
1147
                         (setq ratterms (cons term ratterms)))
 
1148
                        ((rib term s))
1133
1149
                        (t (return nil))))
1134
1150
;;;Function RIB sets up the values of BPTU and BPTD
1135
1151
                (cond ((car plf)
1136
1152
                       (setq bptu (subst (car plf) 'x* bptu))
1137
1153
                       (setq bptd (subst (car plf) 'x* bptd))
 
1154
                       (setq ratterms (subst (car plf) 'x* ratterms))
1138
1155
                       t)        ;CROCK, CROCK. This is TERRIBLE code.
1139
1156
                      (t t))
1140
1157
;;;If there is BPTU then CSEMIUP must succeed.
1141
1158
;;;Likewise for BPTD.
1142
 
                (cond (bptu (cond ((setq upans (csemiup (m+l bptu) d var)))
1143
 
                                  (t nil)))
 
1159
                (cond (ratterms (setq ratans 
 
1160
                                      (defint (m// (m+l ratterms) d) var '$minf '$inf)))
 
1161
                      (t (setq ratans 0)))
 
1162
                (cond (bptu (setq upans (csemiup (m+l bptu) d var)))
1144
1163
                      (t (setq upans 0)))
1145
 
                (cond (bptd (cond ((setq downans (csemidown (m+l bptd) d var)))
1146
 
                                  (t nil)))
 
1164
                (cond (bptd (setq downans (csemidown (m+l bptd) d var)))
1147
1165
                      (t (setq downans 0))))
1148
 
           (sratsimp (m* '$%pi (m+ upans (m- downans))))))))
 
1166
           
 
1167
           (sratsimp (m+ ratans
 
1168
                         (m* '$%pi (m+ upans (m- downans)))))))))
1149
1169
 
1150
1170
 
1151
1171
(defun evenfn (e var)
1576
1596
 
1577
1597
;; integrate(a*sc(r*x)^k/x^n,x,0,inf).
1578
1598
(defun ssp (exp)
1579
 
  (prog (u n c)
1580
 
     (when (notinvolve exp '(%sin %cos))
 
1599
  (prog (u n c arg)
 
1600
     ;; Get the argument of the involved trig function.
 
1601
     (when (null (setq arg (involve exp '(%sin %cos))))
1581
1602
       (return nil))
1582
1603
     ;; I don't think this needs to be special.
1583
1604
     #+nil
1584
1605
     (declare (special n))
1585
 
     ;; Replace (1-cos(x)^2) with sin(x)^2.
1586
 
     (setq exp ($substitute (m^t `((%sin) ,var) 2.)
1587
 
                            (m+t 1. (m- (m^t `((%cos) ,var) 2.)))
1588
 
                            exp))
 
1606
     ;; Replace (1-cos(arg)^2) with sin(arg)^2.
 
1607
     (setq exp ($substitute ;(m^t `((%sin) ,var) 2.)
 
1608
                            ;(m+t 1. (m- (m^t `((%cos) ,var) 2.)))
 
1609
                            ;; The code from above generates expressions with
 
1610
                            ;; a missing simp flag. Furthermore, the 
 
1611
                            ;; substitution has to be done for the complete
 
1612
                            ;; argument of the trig function. (DK 02/2010)
 
1613
                            `((mexpt simp) ((%sin simp) ,arg) 2)
 
1614
                            `((mplus) 1 ((mtimes) -1 ((mexpt) ((%cos) ,arg) 2)))
 
1615
                            exp))
1589
1616
     (numden exp)
1590
1617
     (setq u nn*)
1591
1618
     (cond ((and (setq n (findp dn*))
1639
1666
              (m^ r (m+ n -1))
1640
1667
              `((%signum) ,r)
1641
1668
              recursion)
1642
 
          ;; Recursion failed.  Return the integrand
1643
 
          (let ((integrand (div (pow `((%sin) ,(m* r var))
1644
 
                                     k)
1645
 
                                (pow var n))))
 
1669
          ;; Recursion failed.  Return the integrand
 
1670
          ;; The following code generates expressions with a missing simp flag 
 
1671
          ;; for the sin function. Use better simplifying code. (DK 02/2010)
 
1672
;         (let ((integrand (div (pow `((%sin) ,(m* r var))
 
1673
;                                    k)
 
1674
;                               (pow var n))))
 
1675
          (let ((integrand (div (power (take '(%sin) (mul r var))
 
1676
                                       k)
 
1677
                                (power var n))))
1646
1678
            (m* mult
1647
1679
                `((%integrate) ,integrand ,var ,ll ,ul)))))))
1648
1680
 
1995
2027
          #+nil
1996
2028
          ((not (equal ($asksign denom) '$zero))
1997
2029
           0)
1998
 
          ((equal ($asksign denom) '$zero)
 
2030
          ((equal ($csign denom) '$zero)
1999
2031
           '$undefined)
2000
2032
          (t (intsubs exp ll ul)))))
2001
2033
 
2272
2304
    (cond (result (sratsimp (m* (m- '$%i) result)))
2273
2305
          (t nil))))
2274
2306
 
 
2307
;; Evaluates the contour integral of GRAND around the unit circle
 
2308
;; using residues.
2275
2309
(defun unitcir (grand var)
2276
2310
  (numden grand)
2277
 
  (let ((result (princip (res nn* dn* #'(lambda (pt)
2278
 
                                          (ratgreaterp 1 (cabs pt)))
2279
 
                              #'(lambda (pt)
2280
 
                                  (alike1 1 (cabs pt)))))))
2281
 
    (cond (result (m* '$%pi result))
2282
 
          (t nil))))
 
2311
  (let* ((sgn nil)
 
2312
         (result (princip (res nn* dn* 
 
2313
                               #'(lambda (pt)
 
2314
                                   ;; Is pt stricly inside the unit circle?
 
2315
                                   (setq sgn (let ((limitp nil))
 
2316
                                               ($asksign (m+ -1 (cabs pt)))))
 
2317
                                   (eq sgn '$neg))
 
2318
                               #'(lambda (pt)
 
2319
                                   ;; Is pt on the unit circle?  (Use
 
2320
                                   ;; the cached value computed
 
2321
                                   ;; above.)
 
2322
                                   (prog1
 
2323
                                       (eq sgn '$zero)
 
2324
                                       (setq sgn nil)))))))
 
2325
    (when result
 
2326
      (m* '$%pi result))))
2283
2327
 
2284
2328
 
2285
2329
(defun logx1 (exp ll ul)
2802
2846
              ;; Make the substitution y=1/x.  If the integrand has
2803
2847
              ;; exactly the same form, the answer has to be 0.
2804
2848
              (return 0.))
 
2849
             ((setq ans (logx1 exp ll ul))
 
2850
              (return ans))
2805
2851
             ((setq ans (antideriv exp))
2806
2852
              ;; It's easy if we have the antiderivative.
2807
 
              (return (intsubs ans ll ul)))
2808
 
             ((setq ans (logx1 exp ll ul))
2809
 
              (return ans)))
 
2853
              ;; but intsubs sometimes gives results containing %limit
 
2854
              (return (intsubs ans ll ul))))
2810
2855
       ;; Ok, the easy cases didn't work.  We now try integration by
2811
2856
       ;; parts.  Set ANS to f(x).
2812
2857
       (setq ans (m// exp `((%log) ,arg)))