Fourier series of a 2step function?



It could be better, but it's late.

Rich


----- Original Message ----- 
From: "Richard Hennessy" <rvh2007 at comcast.net>
To: "Edwin Woollett" <woollett at charter.net>; "maxima mailing list" <maxima at math.utexas.edu>
Cc: <dastew at sympatico.ca>
Sent: Friday, September 19, 2008 4:13 PM
Subject: Re: [Maxima] Fourier series of a 2step function?


I was working on declaring functions peicewise and I came up with this (probably not very good solution).  I have
defined some functions that allow you to declare a function in sections (pw) like you want and integrate (pwdefint and
pwintegrate) and differentiate (pwdiff).

(%i1) set_display('ascii)$
(%i2) pwintegrate(expr,x):=
block
    ([__i, __k,__l],
        if matrixp(expr) then
            (
                __l:[],
                if length(expr[1]) = 11 then
                (
                    for __i: 1 thru length(expr) do __l:endcons([expr[__i,5],expr[__i,7], integrate(expr[__i,11],
x)],__l),
                    pw(__l,x)
                )
            )
        else
            integrate(expr,x)
    )$
(%i3) pwdefint(expr,x,a,b):=
block
    ([__i, __s],
        if matrixp(expr) then
            (
                if length(expr[1]) = 11 then
                (
                    __s:0,
                    for __i: 1 thru length(expr) do __s:__s+integrate(expr[__i,11],
x,max(a,expr[__i,5]),min(b,expr[__i,7])),
                    __s
                )
            )
        else
            integrate(expr,x,a,b)
    )$
(%i4) pwdiff(expr,x,n):=
block
    ([__i, __k,__l],
        if matrixp(expr) then
            (
                __l:[],
                if length(expr[1]) = 11 then
                (
                    for __i: 1 thru length(expr) do __l:endcons([expr[__i,5],expr[__i,7], diff(expr[__i,11], x,n)],__l),
                    pw(__l,x)
                )
            )
        else
            diff(expr,x,n)
    )$

(%i5)
pw(L,x):=
block
([__M,__k],
    if listp(L[1]) then
        (
            __M:matrix(["If", x,"in","[",L[1][1], ",", L[1][2],"]",'pw(x),"=",L[1][3]]),
            for __k: 2 thru length(L) do (if listp(L[__k]) then if length(L[__k])=3 then __M:addrow(__M,["If",
x,"in","[",L[__k][1], "," ,L[__k][2],"]",'pw(x),"=", L[__k][3]]))
        )
    else
    (
        __M:matrix(["If",x,"in","[",L[1], ",", L[3],"]",'pw(x),"=",L[2]]),
        for __k:2 thru length(L)-3 step 2 do
            (
            __M:addrow(__M,["If",x,"in","[",L[1+__k],",",L[3+__k],"]",'pw(x),"=",L[2+__k]])
            )
    ),
    __M
)$

(%i6) f:[-inf, 0, -15, x, -10, sin(x^2)*x*20+cos(3*x)*10, -5, x^2+x+6, -2, 10*x, 0, x, 10, 0, inf]$

(%i7) pw(f,x);
                                 [ If  x  in  [  - inf  ,  - 15  ]  pw(x)  =              0              ]
                                 [                                                                       ]
                                 [ If  x  in  [  - 15   ,  - 10  ]  pw(x)  =              x              ]
                                 [                                                                       ]
                                 [                                                      2                ]
                                 [ If  x  in  [  - 10   ,  - 5   ]  pw(x)  =  20 x sin(x ) + 10 cos(3 x) ]
                                 [                                                                       ]
(%o7)                            [                                                     2                 ]
                                 [ If  x  in  [   - 5   ,  - 2   ]  pw(x)  =          x  + x + 6         ]
                                 [                                                                       ]
                                 [ If  x  in  [   - 2   ,   0    ]  pw(x)  =             10 x            ]
                                 [                                                                       ]
                                 [ If  x  in  [    0    ,   10   ]  pw(x)  =              x              ]
                                 [                                                                       ]
                                 [ If  x  in  [   10    ,  inf   ]  pw(x)  =              0              ]

(%i8) pwintegrate(pw(f,x),x);
                                  [ If  x  in  [  - inf  ,  - 15  ]  pw(x)  =             0             ]
                                  [                                                                     ]
                                  [                                                        2            ]
                                  [                                                       x             ]
                                  [ If  x  in  [  - 15   ,  - 10  ]  pw(x)  =             --            ]
                                  [                                                       2             ]
                                  [                                                                     ]
                                  [                                            10 sin(3 x)           2  ]
                                  [ If  x  in  [  - 10   ,  - 5   ]  pw(x)  =  ----------- - 10 cos(x ) ]
                                  [                                                 3                   ]
                                  [                                                                     ]
                                  [                                                  3    2             ]
(%o8)                             [                                                 x    x              ]
                                  [ If  x  in  [   - 5   ,  - 2   ]  pw(x)  =       -- + -- + 6 x       ]
                                  [                                                 3    2              ]
                                  [                                                                     ]
                                  [                                                         2           ]
                                  [ If  x  in  [   - 2   ,   0    ]  pw(x)  =            5 x            ]
                                  [                                                                     ]
                                  [                                                        2            ]
                                  [                                                       x             ]
                                  [ If  x  in  [    0    ,   10   ]  pw(x)  =             --            ]
                                  [                                                       2             ]
                                  [                                                                     ]
                                  [ If  x  in  [   10    ,  inf   ]  pw(x)  =             0             ]

(%i9) pwdefint(pw(f,x),x,-inf,inf), numer,ratprint=false;
(%o9)                                                        7.250095748755271

Rich



----- Original Message ----- 
From: "Edwin Woollett" <woollett at charter.net>
To: "maxima mailing list" <maxima at math.utexas.edu>
Cc: <dastew at sympatico.ca>
Sent: Thursday, September 18, 2008 12:01 PM
Subject: Re: [Maxima] Fourier series of a 2step function?


On Sept. 17, Edwin Woollett wrote:

>(%i7) bn : ( (1/10)*(integrate( -5*sin(n*%pi*x/10), x, -5, 0 ) +
>            integrate( 5*sin(n*%pi*x/10), x, 0, 5 )  ),
>            ratsimp(%%) );
>(%o7) -( 10*cos(%pi*n/2) - 10 )/(%pi*n)
>/* get a Maxima function using expressions bn  */
>(%i8) define( b(n), bn );
>(%o8) b(n) := -( 10*cos(%pi*n/2) - 10 )/(%pi*n)
>
>/* what are the first seven coefficients?  */
>
>(%i9) (bnlist : [], for i thru 7 do bnlist : cons( b(i), bnlist ),
>                                reverse( bnlist )  );
>(%o9) [ 10/%pi,  10/%pi,  10/(3*%pi),  0,
>                       2/%pi,  10/(3*%pi),  10/(7*%pi)  ]
>
getting that list of b(n) values was too complicated!!

Just use map and makelist:

(%i9) map( 'b, makelist( i, i, 1, 7 ) );
(%o9) [10/%pi, 10/%pi, 10/(3*%pi),
               0, 2/%pi, 10/(3*%pi), 10/(7*%pi) ]

Likewise with fs(nmax):

(%i10) fs(nmax) :=  sum( b(m)*sin(m*%pi*x/10), m, 1, nmax )$
(%i11) map( 'fs, [ 1, 2, 3 ] );

(%o11) [ 10*sin(%pi*x/10)/%pi,
              10*sin(%pi*x/5)/%pi + 10*sin(%pi*x/10)/%pi,
        10*sin(3*%pi*x/10)/(3*%pi) + 10*sin(%pi*x/5)/%pi +
                                           10*sin(%pi*x/10)/%pi ]

Ted Woollett


_______________________________________________
Maxima mailing list
Maxima at math.utexas.edu
http://www.math.utexas.edu/mailman/listinfo/maxima

_______________________________________________
Maxima mailing list
Maxima at math.utexas.edu
http://www.math.utexas.edu/mailman/listinfo/maxima