Re: [Maxima-discuss] fourie.mac package fails for simple square wave

[Barton W. <wi...@un...>]

Here is a modest reworking of the function fourcoeff. It integrates from -p to 0 and from o to p separately. Possibly doing so eliminates the need for detecting even or odd integrands.

It doesn't use the function adefint. I'm not confident that I understand adefint—the way it disposes of the absolute value function seems sketchy. ​

The option variables  sinnpiflag and cosnpiflag activate the simplifications sin(n %pi) = 0 and cos(n %pi) = (-1)^n.  Guess: the purpose of these flags is a clumsy workaround to the bug

(%i3) foursimple(integrate(sin(x)*sin(n*x),x,0,%pi)), sinnpiflag : true;
(%o3) 0

When sinnpiflag is false, the denominator tells an alert reader that n=1 is a special case; for example

(%i3) foursimple(integrate(sin(x)*sin(n*x),x,0,%pi)), sinnpiflag : true;
(%o3) 0

(%i4) foursimple(integrate(sin(x)*sin(n*x),x,0,%pi)), sinnpiflag : false;
(%o4) -sin(%pi*n)/((n-1)*(n+1))

My modest rewrite makes the assumption that n is an integer. That effectively removes the purpose of  the flags sinnpiflag and cosnpiflag.

A good fix to this bug is to (a) have integrate return a conditional expression (b) have fourcoeff special case trigonometric polynomials (c) ?

load("unwind_protect");

fourcoeff(f%,x,p) := block([a0,an,bn,n,w,sgn,cntx : supcontext()],
      unwind_protect((
          assume(n >= 1),
      declare(n,integer),
      w : n*%pi*x/p,
      sgn : asksign(p),
      if sgn = 'pos then assume(p > 0) elseif sgn = 'neg then assume(p < 0),
      a0 : (defint(f%,x,-p,0) + defint(f%,x,0,p))/(2*p),
      an : (defint(f%*cos(w),x,-p,0) + defint(f%*cos(w),x,0,p))/p,
      bn : (defint(f%*sin(w),x,-p,0) + defint(f%*sin(w),x,0,p))/p,
      append(ldisp(a[0] = a0),ldisp(a[n] = an),ldisp(b[n] = bn))),
   killcontext(cntx)))$

--Barton
________________________________
From: Barton Willis <wi...@un...>
Sent: Friday, March 15, 2024 21:18
To: Majzoub, Eric <eri...@um...>; rob...@gm... <rob...@gm...>; Barton Willis <wi...@un...>
Cc: max...@li... <max...@li...>
Subject: Re: [Maxima-discuss] fourie.mac package fails for simple square wave


(%i1) fourier_coeffs(e,x,p) := block([n : new_variable('integer), cntx : supcontext(), z : gensym(), w,a0,an,bn],
          unwind_protect((
                 assume(n > 0),
                 w : n*%pi* x / p,
                 a0  : (defint(e,x,-p,0) + defint(e,x,0,p))/(2*p),
                 an :  (defint(e*cos(w),x,-p,0) + defint(e*cos(w),x,0,p))/p,
                 bn :  (defint(e*sin(w),x,-p,0) + defint(e*sin(w),x,0,p))/p,
                  [a0,an,bn]),
              killcontext(cntx)))$

(%i2) (load(unwind_protect),load(to_poly))$

(%i3) fourier_coeffs(1-2*unit_step(x),x,1);
(%o3) [0,0,(2*(-1)^%z1)/(%pi*%z1)-2/(%pi*%z1)]

This is a reported bug, I think:

(%i7) fourier_coeffs(cos(2*%pi*x),x,1);
(%o7) [0,0,0]

--Barton
________________________________
From: Barton Willis via Maxima-discuss <max...@li...>
Sent: Friday, March 15, 2024 17:17
To: Majzoub, Eric <eri...@um...>; rob...@gm... <rob...@gm...>
Cc: max...@li... <max...@li...>
Subject: Re: [Maxima-discuss] fourie.mac package fails for simple square wave

Caution: Non-NU Email

Here is a way to the simplification unit_step(X) + unit_step(-X) --> 1 without using matchdeclare. This simplification is OK for X real and nonzero.  Maybe the freeof(x,ax) check isn't wanted.

(%i1) load("opsubst")$

(%i2) unit_step_reduce(e,x) := block([a : gatherargs(e,'unit_step)],
           for ax in a do (
              ax : first(ax),
              if equal(imagpart(ax),0) and not freeof(x,ax) then (
                  e : ratsubst(1, unit_step(ax)+unit_step(-ax),e))),
          e)$

OK:

(%i13)      unit_step_reduce(x+unit_step(x) + unit_step(-x),x);
(%o13)      x+1

OK:

(%i14)      unit_step_reduce(x + unit_step(x) + unit_step(-x),z);
(%o14)      unit_step(x)+x+unit_step(-x)


Not wrong, but maybe not what you all want:

(%i15)      unit_step_reduce(x + 5*unit_step(x) - unit_step(-x),x);
(%o15)      6*unit_step(x)+x-1

--Barton
________________________________
From: Majzoub, Eric <eri...@um...>
Sent: Friday, March 15, 2024 15:57
To: rob...@gm... <rob...@gm...>
Cc: max...@li... <max...@li...>
Subject: Re: [Maxima-discuss] fourie.mac package fails for simple square wave

Caution: Non-NU Email


Hi Robert: Thanks for the reply.

I tried using tellsimpafter by itself and your suggestion with
matchdeclare as well. But the code in fourie.mac still returns the
wrong answer -- the coefficients are returned as functions of x, which
should integrated out.

In any case, the pdefourier package works and there is code
infrastructure for piecewise continuous functions there.

-Eric


_______________________________________________
Maxima-discuss mailing list
Max...@li...
https://lists.sourceforge.net/lists/listinfo/maxima-discuss