[Maxima] Fourier Series using fourie.mac

[Richard Hennessy]

My package pw.mac can do the integration necessary to solve the Fourier
series directly without using "fourie.mac" as long as you remember how to do
it.
Here is a solution using pw.mac.  All the needed functions are documented on
my web site.

This result is based on the fact that for x in the interval (-%pi/2, %pi/2)

abs(sin(x)) = pw([-%pi/2,-sin(x),0,sin(x),%pi/2],x) = signum(x) * sin(x)

(%i1) (load(pw), ratprint:false);
(%o1) false
(%i2) now:elapsed_run_time();
(%o2) 2.95
(%i3) declare(n,integer);
(%o3) done
(%i4) assume(n>0);
(%o4) [n>0]
(%i5) f:[-%pi/2, -sin(x), 0, sin(x), %pi/2];
(%o5) [-%pi/2,-sin(x),0,sin(x),%pi/2]
(%i6) a0 : pwdefint(pw(f,x), x, -%pi/2, %pi/2)/%pi;
(%o6) 2/%pi
(%i7) an : (s:1/(%pi/2)*pwdefint(pw(f,x)*cos(n*%pi*x/(%pi/2)), x, -%pi/2,
%pi/2),radcan(s));
(%o7)
-((4*n-2)*cos((2*%pi*n+%pi)/2)+(-4*n-2)*cos((2*%pi*n-%pi)/2)+4)/(4*%pi*n^2-%pi)
(%i8) bn : (s:(1/(%pi/2))*pwdefint(pw(f,x)*sin(n*%pi*x/(%pi/2)), x, -%pi/2,
%pi/2), radcan(s));
(%o8) 0
(%i9) (define(a(n), an),define(b(n), bn))$
(%i10) g(x, terms)
:=a0+sum(b(m)*sin(m*%pi*x/(%pi/2))+a(m)*cos(m*%pi*x/(%pi/2)), m, 1, terms)$
(%i11) define(h(x), periodic(pw(f,x),x,-%pi/2,%pi/2))$
(%i12) n:5;
(%o12) 5
(%i13) load(draw)$
(%i14) draw2d(
    color=brown,
    yrange=[-0,1],
    label_alignment=left,
    explicit(h(x),x,-%pi,%pi),
    color=blue,
    label(["f(x)=abs(sin(x))",-2,.3]),
    label([concat("with fourier series for n = ",n), -2,.25]),
    explicit(g(x,5),x,-%pi,%pi))$
(%i15) elapsed_run_time()-now;
(%o15) 2.44


To prove the identity

(%i1) display2d:false;

(o1) false
(%i2) abs(sin(x));
(o2) abs(sin(x))
(%i3) abs2signum(%);
(o3) sin(x)*signum(sin(x))
(%i4) periodic(sin(x),x,-%pi/2,%pi/2);
(o4) -cos(%pi*((x+%pi/2)/%pi-floor((x+%pi/2)/%pi)))
(%i5) integrate(%,x);
(o5) -'integrate(cos(%pi*((x+%pi/2)/%pi-floor((x+%pi/2)/%pi))),x)
(%i6) changevar(%,y-(x+%pi/2)/%pi,y,x);
(%o6) sin(%pi*floor(y)-%pi*y)
(%i7) subst((x+%pi/2)/%pi,y,%);
(o7) sin(%pi*floor((x+%pi/2)/%pi)-x-%pi/2)
(%i8) abs(sin(x))=pw([-%pi/2,-sin(x),0,sin(x),%pi/2],x);
(o8) abs(sin(x)) =
sin(x)*(signum(x)-signum(x-%pi/2))/2-sin(x)*(signum(x+%pi/2)-signum(x))/2
(%i9) simp_assuming(%,x<%pi/2,x>-%pi/2);
(o9) abs(sin(x)) = (signum(x)+1)*sin(x)/2-(1-signum(x))*sin(x)/2
(%i10) radcan(%);
(o10) abs(sin(x)) = signum(x)*sin(x)
(%i11)

You can get pw.mac from my site.
http://mysite.verizon.net/res11w2yb/id2.html

Hope this helps.

Richard Hennessy



2009/7/28 Leo Butler <l.butler at ed.ac.uk>

>
>
> On Tue, 28 Jul 2009, Rafa Topolnicki wrote:
>
> < Leo Butler pisze:
> < >
> < > On Mon, 27 Jul 2009, Rafa? Topolnicki wrote:
> < >
> < > < Hi,
> < > <
> < > < During testing Fourier Series Expansion in maxima I found a problem I
> < > < can't cope with.
> < > <
> < > < (%i3)                            load(fourie)
> < > < (%o3)         /usr/share/maxima/5.18.1/share/calculus/fourie.mac
> < > < (%i4)                         f(x) := abs(sin(x))
> < > < (%o4)                         f(x) := abs(sin(x))
> < > < (%i5)                    flist : fourier(f(x), x, %pi)
> < > <                                           2
> < > < (%t5)                              a  = ---
> < > <                                      0   %pi
> < > <
> < > <                       cos(%pi n)   cos(%pi n)      1         1
> < > <                    2 (---------- - ---------- + ------- - -------)
> < > <                        2 n + 2      2 n - 2     2 n + 2   2 n - 2
> < > < (%t6)        a  = -----------------------------------------------
> < > <                n                         %pi
> < > <
> < >
> < > The easiest thing to do would be to add a line
> < >
> < > %t6 : a[n] = if n=1 then 'limit( subst(n=_n,rhs(%t6)),_n,n ) else
> < > %rhs(%t6);
>
> You need to overcome the evaluation rules in sum. Here is one way to do
> this using a lambda-function:
>
> (%i2) load("fourie");
> (%o2) "/home/work/maxima/maxima-5.18.1-clisp/share/calculus/fourie.mac"
> (%i3) f(x):=abs(sin(x));
> (%o3) f(x):=abs(sin(x))
> (%i4) flist : fourier(f(x),x,%pi);
> (%t4) a[0] = 2/%pi
>
> (%t5) a[n] =
> 2*(cos(%pi*n)/(2*n+2)-cos(%pi*n)/(2*n-2)+1/(2*n+2)-1/(2*n-2))/%pi
>
> (%t6) b[n] = 0
>
> (%o6) [%t4,%t5,%t6]
> (%i7) %t5 : a[n] = 'apply(lambda([s], if s=1 then 0 else
> rhs(''%t5)),[n]);
> (%o7) a[n] = 'apply(lambda([s],
>                           if s = 1 then 0
>                               else rhs(
>                               a[n] = 2*(cos(%pi*n)/(2*n+2)
>                                        -cos(%pi*n)/(2*n-2)+1/(2*n+2)
>                                        -1/(2*n-2))
>                                    /%pi)),[n])
> (%i8) fourexpand(flist,y,%pi,3)$
> (%i9) ev(%,nouns);
> (%o9) 2/%pi-4*cos(2*y)/(3*%pi)
>
>
> Leo
>
> --
> The University of Edinburgh is a charitable body, registered in
> Scotland, with registration number SC005336.
>
> _______________________________________________
> Maxima mailing list
> Maxima at math.utexas.edu
> http://www.math.utexas.edu/mailman/listinfo/maxima
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://www.math.utexas.edu/pipermail/maxima/attachments/20090728/be2b444e/attachment.htm