integration of atan2
- Subject: integration of atan2
- From: Barton Willis
- Date: Thu, 6 Jun 2013 19:53:27 +0000
Comments are welcomed; my second waffle:
(1) delete the atan2 --> atan conversion in integrator (defined in sin.lisp)
(2) append $atan2 to the list of functions risch attempts to integrate (append to rischform)
(3) revised trigin1 (defined in risch.lisp)
Before and after:
Nounform :(
(%i1) integrate(exp(cos(x)) * cos(sin(x)),x,0,2*%pi);
(%o1) integrate(%e^cos(x)*cos(sin(x)),x,0,2*%pi)
Bogus (integrand is positive):
(%i2) expand(integrate(atan2(1,x),x,-1,1));
(%o2) -%pi
Invalid for x < 0:
(%i3) integrate(atan2(x,x),x);
(%o3) (%pi*x)/4
Experimental code (see below)
(%i4) load("F:/larry.lisp")$
(%i5) load("F:/sin.lisp")$
Correct, (I think) and not a nounform:
(%i6) integrate(exp(cos(x)) * cos(sin(x)),x,0,2*%pi);
(%o6) 2*%pi
Correct:
(%i7) expand(integrate(atan2(1,x),x,-1,1));
(%o7) %pi
OK on reals:
(%i8) integrate(atan2(x,x),x);
(%o8) x*atan2(x,x)
(%i9) diff(%,x);
(%o9) atan2(x,x)
--Barton
defun rischform (l)
(cond ((or (atom l) (alike1 intvar l) (freeof intvar l)) nil)
((polylogp l)
(if (and (integerp (car (subfunsubs l)))
(signp g (car (subfunsubs l))))
(rischform (car (subfunargs l)))
(setq operator t)))
((atom (caar l))
(case (caar l)
((%sin %cos %tan %cot %sec %csc)
(setq trigint t $exponentialize t)
(rischform (cadr l)))
((%asin %acos %atan %acot %asec %acsc $atan2) ; <--- append $atan2
(setq trigint t $logarc t)
(rischform (cadr l)))
((%sinh %cosh %tanh %coth %sech %csch)
(setq hypertrigint t $exponentialize t)
(rischform (cadr l)))
((%asinh %acosh %atanh %acoth %asech %acsch)
(setq hypertrigint t $logarc t)
(rischform (cadr l)))
((mtimes mplus mexpt rat %erf %log)
(mapc #'rischform (cdr l)))
(t (setq operator (caar l)))))
(t (setq operator (caar l)))))
(defun attempt-analytic-continuation (e)
(let (($opsubst t))
;; atan2(cos(x),sin(x)) --> x
(setq e ($substitute #'(lambda (a b)
(if (and (consp a) (consp (car a)) (eq (mop a) '%cos)
(consp b) (consp (car b)) (eq (mop b) '%sin)
(alike (second a) (second b)))
(second a) (take '($atan2) a b))) '$atan2 e))
;; atan(tan(x)) --> x
(setq e (let (($triginverses '$all) ($exponentialize nil) ($demoivre t)) ($trigreduce e))) ; trigreduce does sin(x)/cos(x) -> tan(x).
;; Convert incomplete_gamma(0,x) to exponential_ei functions; and convert
;; exponential_ei(z) --> %gamma + (log(z) - log(1/z))/2 + z * hypergeometric([1,1],[2,2],z).
(let (($gamma_expand t) (dosimp t) ($exponentialize t) ($logexpand '$all))
(setq e (sratsimp (resimplify e))) ; gotta resimplify to convert from incomplete_gamma to expintegral_ei.
(setq e ($substitute #'(lambda (z)
(add '%gamma
(div (sub (take '(%log) z) (take '(%log) (div 1 z))) 2)
(mul z (take '($hypergeometric) (list '(mlist) 1 1) (list '(mlist) 2 2) z)))) '%expintegral_ei e))
(sratsimp e))))
(defun trigin1 (*exp var)
(let ((rischp var) (rp-polylogp t) ($logarc nil) ($exponentialize nil)
(result (hypertrigint1 *exp var nil)))
;; When result has a integral nounform, set *in-risch-p* to true and resimplify.
;; Presumably setting *in-risch-p*to true avoids an infinite loop.
(if (isinop result '%integrate)
(let ((*in-risch-p* t)) (setq result (resimplify result))))
(setq result (let ((generate-atan2 t)) ($rectform result)))
(sratsimp (attempt-analytic-continuation (sratsimp result)))))
Replaces:
(defun trigin1 (*exp var)
(let ((yyy (hypertrigint1 *exp var nil)))
(setq yyy (div ($expand ($num yyy))
($expand ($denom yyy))))
(let ((rischp var) (rp-polylogp t) $logarc $exponentialize result)
(setq result (sratsimp (if (and (freeof '$%i *exp) (freeof '$li yyy))
($realpart yyy)
($rectform yyy))))
;; The result can contain solveable integrals. Look for this case.
(if (isinop result '%integrate)
;; Found an integral. Evaluate the result again.
;; Set the flag *in-risch-p* to make sure that we do not call
;; rischint again from the integrator. This avoids endless loops.
(let ((*in-risch-p* t))
(meval (list '($ev) result '$nouns)))
result))))
________________________________