#4660 Incorrect integration leads to incorrect series in fourie package

[Devin Greene]

-------------------------------------------------------------
Maxima version: "5.45.1"
Maxima build date: "2021-12-06 23:50:08"
Host type: "x86_64-pc-linux-gnu"
Lisp implementation type: "GNU Common Lisp (GCL)"
Lisp implementation version: "GCL 2.6.12"
User dir: "/home/devingreene/.maxima"
Temp dir: "/tmp"
Object dir: "/home/devingreene/.maxima/binary/5_45_1/gcl/GCL_2_6_12"
Frontend: false
-------------------------------------------------------------

The integral over the interval [-pi,pi] of cos(n*x), where n is an integer is zero when n != 0 and 2*pi otherwise. However:

(%i1) declare(n,integer)$

(%i2) integrate(cos(n*x),x,-%pi,%pi); /* This should give a case prompt regarding n. */
(%o2)                                  0

I don't know the internals of the fourie package, but consider this:

(%i3) load(fourie)$

(%i4) fourier(cos(x),x,%pi); /* Incorrect a_1 */
(%t4)                               a  = 0
                                     0

(%t5)                               a  = 0
                                     n

(%t6)                               b  = 0
                                     n

(%o6)                           [%t4, %t5, %t6]

Obviously a_1 should be 1, not 0. I believe that the integration problem mentioned above is what is causing this.

[Barton Willis]

I think this has been suggested before: make the antiderivative of cos to be (x -> x sinc(x)). Doing so gives

(%i1) :lisp(setf (get '%cos 'integral) '((x) ((mtimes) x ((%sinc) x)))))

((X) ((MTIMES) X ((%SINC) X)))
(%i1) :lisp(def-simplifier sinc (x) (if (eql 0 x) 1 (give-up))))

SIMP-%SINC
(%i1) xxx : integrate(cos(n*x),x,-%pi,%pi);
(%o1)                 %pi sinc(%pi n) + %pi sinc(- %pi n)
(%i2) subst(n=0,xxx);
(%o2)                                2 %pi
(%i3) subst(n=3,xxx);
(%o3)                 %pi sinc(3 %pi) + %pi sinc(- 3 %pi)

Boosting up the sinc simp code to recognize that sinc(integer x %pi) = 0 is easy to do.

Switching to a sinc antiderivative might be OK if it were localized to the definite integral code, but doing this in general would almost surely cause problems with users and problems with vanishing expressions that don't simplify to zero.

[Emmanuel Roque]

Dear Devin Greene,

In case you are interested, our package pdefourier solves this problem and many other potential problems when the input contains trigonometric expressions or a product of a polynomial and a trigonometric expression

https://github.com/emmanuelroque/pdefourier

I hope you find it helpful.

There is also an arxiv preprint where we explain why even Mathematica has many issues regarding Fourier series.

https://arxiv.org/abs/2004.02604

Best,
Emmanuel

[Barton Willis]

Possibility: Define

(def-simplifier sinc (x) 
   (cond ((eql x 0) 1)
         (t 
           (let* ((npi (multiplep x '$%pi))
                  (n (if npi (div x '$%pi) nil)))
                  (cond ((and ($featurep n '$integer) (eq t (mnqp n 0))) 0)
                       (t (give-up)))))))

(defmvar *in-defint* nil)

(defun cos-integral (x)  
   (if *in-defint* 
    (mul x (ftake '%sinc x))
    (ftake '%sin x)))

(putprop '%cos  `((x) ,'cos-integral) 'integral)

And at the top-level of defint, set *in-defint* to true. Doing so, I get:

(%i3) declare(n,integer);
(%o3)                                done
(%i4) xxx : integrate(cos(n*x),x,-%pi,%pi);
(%o4)                 %pi sinc(%pi n) + %pi sinc(- %pi n)
(%i5) subst(n=0,%);
(%o5)                                2 %pi
(%i6) assume(n > 0)$

(%i7) xxx : integrate(cos(n*x),x,-%pi,%pi);
(%o7)                                  0

[Barton Willis]

Here is a stand alone workaround for this bug:

(in-package :maxima)

(defmvar $use_sinc_integration t)

(defun cos-integral (x)  
   (if $use_sinc_integration
    (mul x (ftake '%sinc x))
    (ftake '%sin x)))
(putprop '%cos `((x) ,'cos-integral) 'integral)

(defun sin-integral (x)  
   (if $use_sinc_integration
    (div (mul x x (ftake 'mexpt (ftake '%sinc (div x 2)) 2)) 2)
    (neg (ftake '%cos x))))
(putprop '%sin `((x) ,'sin-integral) 'integral)

(def-simplifier sinc (x) 
   (cond ((or (zerop1 x) (eq t (meqp x 0))) (add 1 x)) ; the add 1 ... makes sinc(0.0) = 1.0 (not 1)
         (t 
           (let ((z (ftake '%sin x)))
             (cond ((eq t (mnqp x 0)) (div z x)) ;x =/=0, so sinc(x) = sin(x)/
                   ((great (neg x) x) (ftake '%sinc (neg x))) ; even reflection rule: sinc(-x) = sinc(x)
                   (t (give-up))))))) ; sinc nounform return

Example:

(%i1)   load("sinc.lisp")$

(%i3)   declare(n,integer)$

(%i4)   xxx : integrate(cos(n*x),x,-%pi,%pi);
(xxx)   2*%pi*sinc(%pi*n)

(%i5)   makelist(''xxx,n,0,8);
(%o5)   [2*%pi,0,0,0,0,0,0,0,0]