Menu

#4023 defint(,,0,inf) keyhole integration, patan issue

None
open
nobody
None
5
2022-08-27
2022-08-27
No
assume(sin(phi) > 0, cos(phi) > 0);
g_poles: exp(%i * [phi, -phi, %pi-phi, %pi+phi]);
g: 1/apply("*", makelist(w - pole, pole, g_poles));

defint(g, w, 0, inf);
fullratsimp(exponentialize(%));
(e*num(%))/factor(demoivre(expand((e*denom(%))))), e=exp(-2*%i*phi);

which simplifies to the incorrect

which_is:  %pi/4 * exp(%i*phi) / sin(2*phi);

but should be

gold: %pi / 4 / sin(phi);

Discussion regarding this can be found in thread
[Maxima-discuss] help with maxima-substitute on a plog(-w)

Barton Willis pointed to the function "patan" which
produces an incorrect result here.

patan has two arguments "r" and "i" which are meant to describe
a complex value "r + %i*i"
and patan should return either "nil" if it feels uncomfortable with the situation,
or the principal arg of the given value in the range of ]-%pi .. +pi].
This function is a private part of the defint business, and
was written in such a way as to return the well known pleasant angles
0 %pi/6 %pi/4 etc etc
or as a special case returns phiexpr when given
"r = expr * cos(phiexpr)" "i = expr * sin(phiexpr)"
where "phiexpr" can be some nontrivial expression itself.

The special case is done here in "patan" in csimp2.lisp

  (cond
      ;%...
       ((and (among '%cos r) (among '%sin i))
         (let ((var 'xz) genvar varlist)
            (numden (div* r i))
            (cond ((and (eq (caar nn*) '%cos) (eq (caar dn*) '%sin))
                 (return (cadr nn*)))))))

and will be off by one %pi if the common factor "expr" from above
is negative. Further it can be off by any integer amount of %pi
if phiexpr is not in [-%pi/2..%pi/2]

In the offending test case we have
"r = - cos(phi)" and
"i = - sin(phi)"
which will produce the incorrect return value "phi" instead of "phi - %pi"

A minimal sketch of a fix is to use "atan2" here like this

modified   src/csimp2.lisp
@@ -96,7 +96,7 @@
               (let ((var 'xz) genvar varlist)
        (numden (div* r i))
        (cond ((and (eq (caar nn*) '%cos) (eq (caar dn*) '%sin))

-              (return (cadr nn*)))))))
+              (return (simplify `(($atan2) ,i ,r))))))))
        (setq a ($sign r) b ($sign i))
        (cond ((eq a '$pos) (setq a 1))
         ((eq a '$neg) (setq a -1))

(this is a sketch only, at least use of `ftake' seems to be wanted here)

A more general fix by Barton Willis is is to dispose patan
all-together and to always use atan2 instead. This might lead
to defint not consider other strategies and thus yield results
using some unpleasant and potentially difficult atan2 based expressions.

The current maxima test suite seems to be ok with both fixes.

Discussion


Log in to post a comment.

MongoDB Logo MongoDB