1. Summary
  2. Files
  3. Support
  4. Report Spam
  5. Create account
  6. Log in

Main Page

From piecewisefunc

(Difference between revisions)
Jump to: navigation, search
m
(boo boo)
Line 1: Line 1:
<big>'''Pw.mac'''</big>
<big>'''Pw.mac'''</big>
-
Did you know that you can integrate the abs(x) function?  Did you know you can also differentiate it too?  Well, you can.  This package, pw.mac, can do it.  Pw.mac is the only math software add-on to [http://maxima.sourceforge.net Maxima] for piecewise defined functions. Pw.mac allows Maxima to do algebra, calculus and plotting of '''any piecewise continuous function''' including products, quotients and sums of two or more piecewise functions.  You can integrate piecewise continuous functions with pwint(), simplify them with pwsimp(), simpunitstep(), simpsignum(), create them with piecewise(), between(), iif(), differentiate them with diff(), plot them with plot2d(), plot3d() or draw(), find limits of them with pwlimit(), change them from one form to another with a large collection of conversion functions. Algebra, plotting and calculus of piecewise continuous functions are all easier with pw.mac. You can compute a Fourier series of a piecewise continuous function easily.  If you have a desire to explore piecewise type functions pw.mac is available from sourceforge.net. Somewhat out of date help for pw.mac is available from [http://mysite.verizon.net/res11w2yb/Pw.html here].  The test suite is rtest_pw.mac.  The version of pw.mac is now 6.4 and works with Maxima 5.25.1 and some earlier versions. Pw.mac is a free program and is open source.   
+
Did you know that you can integrate the abs(x) function?  Did you know you can also differentiate it too?  Well, you can.  This package, pw.mac, can do it.  Pw.mac is the only math software add-on to [http://maxima.sourceforge.net Maxima] for piecewise defined functions. Pw.mac allows Maxima to do algebra, calculus and plotting of '''any piecewise continuous function''' including products, quotients and sums of two or more piecewise functions.  You can integrate piecewise continuous functions with pwint(), simplify them with pwsimp(), simpunitstep(), simpsignum(), create them with piecewise(), between(), iif(), differentiate them with diff(), plot them with plot2d(), plot3d() or draw(), find limits of them with pwlimit(), change them from one form to another with a large collection of conversion functions. Algebra, plotting and calculus of piecewise continuous functions are all easier with pw.mac. You can compute a Fourier series of a piecewise continuous function easily.  If you have a desire to explore piecewise type functions pw.mac is available from sourceforge.net. Somewhat out of date help for pw.mac is available from [http://mysite.verizon.net/res11w2yb/Pw.html here].  The test suite is rtest_pw.mac.  The version of pw.mac is now 6.3 and works with Maxima 5.25.1 and some earlier versions. Pw.mac is a free program and is open source.   
To use this package first download it from my site, then unzip the compressed files preferably to the share/contrib directory and then type.
To use this package first download it from my site, then unzip the compressed files preferably to the share/contrib directory and then type.

Revision as of 03:43, 13 November 2011

Pw.mac

Did you know that you can integrate the abs(x) function? Did you know you can also differentiate it too? Well, you can. This package, pw.mac, can do it. Pw.mac is the only math software add-on to Maxima for piecewise defined functions. Pw.mac allows Maxima to do algebra, calculus and plotting of any piecewise continuous function including products, quotients and sums of two or more piecewise functions. You can integrate piecewise continuous functions with pwint(), simplify them with pwsimp(), simpunitstep(), simpsignum(), create them with piecewise(), between(), iif(), differentiate them with diff(), plot them with plot2d(), plot3d() or draw(), find limits of them with pwlimit(), change them from one form to another with a large collection of conversion functions. Algebra, plotting and calculus of piecewise continuous functions are all easier with pw.mac. You can compute a Fourier series of a piecewise continuous function easily. If you have a desire to explore piecewise type functions pw.mac is available from sourceforge.net. Somewhat out of date help for pw.mac is available from here. The test suite is rtest_pw.mac. The version of pw.mac is now 6.3 and works with Maxima 5.25.1 and some earlier versions. Pw.mac is a free program and is open source.

To use this package first download it from my site, then unzip the compressed files preferably to the share/contrib directory and then type.

(%i1) load(pw)$

Some functions that are available are demonstrated here.

(%i2) piecewise([minf,-x^2,0,x^2-x,inf],x) /* pw() is an alias for piecewise */;
                    2                         2
                  (x  - x) (signum(x) + 1)   x  (1 - signum(x))
                  ------------------------ - ------------------
                             2                       2

(%i3) piecewise([minf,-x^2,0,x^2-x,inf],x,'iif);  /* you can simplify this a little with pulliniif() */
                                                             2                                             2
                                      2                     x                     2                       x  - x
(%o3)                   iif(x < 0, - x , iif(equal(x, 0), - --, 0)) + iif(x > 0, x  - x, iif(equal(x, 0), ------, 0))
                                                             2                                               2

(%i4) piecewise([minf,-x^2,0,x^2-x,inf],x,'abs); /* You can express all piecewise continuous functions in terms 
of the abs() function instead of the signum() function if you prefer. */

                                  (2 x - 1) abs(x) - x
(%o4)                             --------------------
                                           2

/* by default piecewise functions are expressed in terms of the signum() and unit_spike() functions but
   between(), iif(), %if(), abs() and if then else can be used too.  
   The "if then else" form is not very useful for further manipulation however. */

(%i5) signum2abs(pwint(%,x));
                                   2
                                  x  (4 abs(x) - 3) - 3 x abs(x)
                                  ------------------------------
                                                12

(%i6) pwint(%o3,x);  /* integrate %o3 and leave it in the default form, you can simplify this with ratsimp() */
                                                   3    2                    3    2
                                                  x    x                    x    x
                                   3             (-- - --) signum(x)    3   -- - --
                                  x  signum(x)    3    2               x    3    2
                                  ------------ + ------------------- - -- + -------
                                        6                  2            6       2

(%i7) periodic(sin(a*x),x,0,2*%pi/a); /* you can write periodic piecewise continuous functions using pw.mac */

                             a x           a x
                 sin(2 %pi (----- - floor(-----)))
                            2 %pi         2 %pi

(%i8) assume(a>0);
                                           [a > 0]

(%i9) intperiodic(sin(a*x),x,0,2*%pi/a);  /* and integrate them at least once */
                              a x           a x
                  cos(2 %pi (----- - floor(-----)))
                              2 %pi         2 %pi
                - ---------------------------------
                                 a

(%i10) pwint((x+signum(x-a))^5,x);  
                                        3       3                                           6      4      2
           5    5                   10 x    10 a                                           x    5 x    5 x
         (x  - a ) signum(x - a) + (----- - -----) signum(x - a) + (x - a) signum(x - a) + -- + ---- + ----
                                      3       3                                            6     2      2


(%i11) pwsimp((x+signum(x-a))^5,x);      /*  you can use pwsimp to show the function in a different but 
                                             equivalent form  */
                                   5                              5
                            (x + 1)  (signum(x - a) + 1)   (x - 1)  (1 - signum(x - a))
                            ---------------------------- + ----------------------------
                                         2                              2
(%i12) pwsimp((signum(x-4)*signum(x-5))+(signum(x-4)-signum(x-5)-1),x);  /*  you can show two forms are the same 
                                                                             with pwsimp() */

                                               0


(%i13) pwint(x^2*pwdelta(x-a),x,b,c); 
       /* pwdelta is a curious animal, similar to the diracdelta() function.  Use with caution.
          See the help for more information */
                            2                  2
                           a  signum(c - a) - a  signum(b - a)
                           -----------------------------------
                                            2
(%i14) pwint(f(x)*pwdelta(x-a),x,minf,inf);  /* pwdelta picks out f(a) */

                               
                                           f(a)
  
(%i15) pwint(f(x)*pwdelta(x-a),x,a,a);  /* but not in this case.  You can't integrate from a to a and 
                                              get a nonzero result with pw.mac */

                                             0

(%i16) ratsimp(pwint(x^2*pwdelta(x-a),x,a-e,a+e)); /* if e is positive then this will pick out a^2 if e < 0 
then you would get -a^2 */

                                2
                               a  signum(e)

(%i17) is(equal(pwint(pwdelta(x-a),x),(signum(x-a)+1)/2));   /* true by definition of pwdelta */

(%o17)                             true

(%i18) simp_given(pwlimit(x^2*(signum(a*x^5+b)-signum(a*x-b)),x,inf),notequal(a,0)) /* pwlimit tries to find limits of expressions
                                                                                          involving signum() functions  */;

(%o18)                             0

(%i19) pushout((s-signum(x-s))^5,'signum);  /* factor out the signum's */
                5                    4              2       3              3       2             4                  5
(%o19)  - signum (x - s) + 5 s signum (x - s) - 10 s  signum (x - s) + 10 s  signum (x - s) - 5 s  signum(x - s) + s

(%i20) linearize(%); /* signum(s)^2=1 unless s = 0, so assume s # 0, This is okay in case the function is continuous. */ 
                        4                     2                                  5       3
(%o20)             - 5 s  signum(x - s) - 10 s  signum(x - s) - signum(x - s) + s  + 10 s  + 5 s

(%i21) pushout(%,'signum);  /* factor out the signum */

                                4       2                       5       3
(%o21)                    (- 5 s  - 10 s  - 1) signum(x - s) + s  + 10 s  + 5 s

(%i22) sin(abs(x-a))*sin(abs(x-b));  /* iif() and pulliniif() can sometimes extend integrate to do problems otherwise not possible */

(%o22)             sin(abs(x - a)) sin(abs(x - b))

(%i23) abs2iif(%); /* convert to iif() form */

(%o23)           sin(iif(x - a > 0, x - a, a - x)) sin(iif(x - b > 0, x - b, b - x))

(%i24) pulliniif(%);  /* do this to pull the abs()'s inside the iif() instead of inside the sin() */

(%o24) iif(x - b > 0, iif(x - a > 0, sin(x - a) sin(x - b), - sin(x - a) sin(x - b)), iif(x - a > 0, - sin(x - a) sin(x - b), 
                                                     sin(x - a) sin(x - b)))

(%i25) pwint(%,x);  /* now integrate and you get the right answer */

                                      sin(b - a) - 2 b cos(b - a)   sin(2 x - b - a) - 2 cos(b - a) x
(%o25)  signum(x - a) (signum(x - b) (--------------------------- - ---------------------------------)
                                                   4                                4

                                                           sin(b - a) - 2 b cos(b - a)   - sin(b - a) - 2 a cos(b - a)
                                          + signum(b - a) (--------------------------- - -----------------------------))
                                                                        4                              4

(%i26) signum2abs(simpspikes(iif2sum(iif(x>a,x-a,-x+a)))); /* here is how to convert from an iif() form to an abs() form */

(%o26)                                                      abs(x-a)

(%i27) pulliniif(abs2iif(max2abs(max(1,x))));  /* max() and min() are not supported forms for expressing piecewise functions
          but they can be converted to other forms.  In this case I convert to an iif() form which is simple and maybe not so obvious 
          but true. Max2abs() and min2abs() are Barton Willis' idea and I use them with his permission in pw.mac */

(%o28)                                                   iif(x - 1 > 0, x, 1)

(%i29) pulliniif(abs2iif(min2abs(min(x^2,x))));  /* min() to iif() */

                                                              2              2
                                                         iif(x  - x > 0, x, x )

(%i30) pwsimp(%o2,x,'ifthen);      /*  you can use pwsimp to convert from one form to another equivalent form in most cases */

   - signum(x) - unit_spike(x) + 1             2                                  signum(x) - unit_spike(x) + 1           2
if ------------------------------- > 0 then - x  elseif equal(x, 0) then 0 elseif ----------------------------- > 0 then x  - x else 0
                  2                                                                             2

(%i31) pwsimp(%o2,x,'abs);  /* abs() form */

                                      (2 x - 1) abs(x) - x
(%o31)                                --------------------
                                               2

(%i32) pwsimp(%o2,x,'iif);  /* iif() form */
                                                    2                   2
(%o32)                                iif(x < 0, - x , 0) + iif(x > 0, x  - x, 0)

(%i33) plot2d([pw([minf,x^2,0,-x^2,inf],x)], [x,-4,4])$ /* this is one inefficient way to plot a piecewise function */

(%i34) h(x):=x^2*between(x,0,inf)-x^2*between(x,minf,0);  /* this is not very efficient for plotting either */
                                               2                       2
                                      h(x) := x  between(x, 0, inf) - x  between(x, minf, 0)
(%i35) plot2d([h(x)],[x,-4,4]);

(%i36) h(x):=''(x^2*between(x,0,inf)-x^2*between(x,minf,0));  /* this is much better */
                                               2                    2
                                              x  (signum(x) + 1)   x  (1 - signum(x))
                                      h(x) := ------------------ - ------------------
                                                      2                    2
(%i37) plot2d([h(x)],[x,-4,4]);

Output of plot2d:

x Plot of h(x)


(%i38) 1/(1+sqrt(1-2*x+x^2)); /* this function is piecewise continuous if you assume sqrt(x^2) = abs(x) */

(%39) pwint(%,x);  /* pwint recognizes this fact */

                                         signum(x - 1) log(x)   log(x)   log(2 - x) signum(x - 1)   log(2 - x)
                                         -------------------- + ------ + ------------------------ - ----------
                                                  2               2                 2                   2

(%40) pwint(sum(x*abs(a*x+b[i]), i, 1, n),x);  /* pwint automatically distributes over sums */

                                           n          3         2     3
                                          ====   2 a x  + 3 b  x     b
                                          \                  i        i
(%o40)                                     >    (---------------- - ----) signum(a x + b )
                                          /             6              2                i
                                          ====                      6 a
                                          i = 1
                  ^

(%i1) pushout(pwint(pw([-10,x,-5,-x^2,0,x^2,5,x,10],x,'pulse), x),'signum,'unit_pulse);  /* preview of new feature (not available yet), 
                                        unit_pulse can be integrated directly */
        3            x               x       3              x       2              x       2
       x  unit_pulse(-)   unit_pulse(- + 1) x    unit_pulse(- + 2) x    unit_pulse(- - 1) x
                     5               5                      5                      5                     x
(%o1)  ---------------- - -------------------- + -------------------- + -------------------- - 25 signum(- + 2)
              3                    3                      2                      2                       5
                                                                          x                   x
                                                               175 signum(- + 1)   175 signum(- - 1)
                                                                          5                   5                  x        625
                                                             - ----------------- + ----------------- + 25 signum(- - 2) + ---
                                                                      12                  12                     5         6

(%i2)  pwint(prod(unit_step(x-a[i]),i, 1, 10),x);  /* I tried this in Maple 14, be prepared to wait a long time */ 
             (x - max(a , a , a , a , a , a , a , a , a , a  )) signum(x - max(a , a , a , a , a , a , a , a , a , a  ))
                       1   2   3   4   5   6   7   8   9   10                   1   2   3   4   5   6   7   8   9   10     x
(%o2)        ----------------------------------------------------------------------------------------------------------- + -
                                                                  2                                                        2

Pw.mac can do algebra, differentiation, integration, plotting, Convolutions, Fourier series and Fourier tranforms of piecewise continuous functions expressed in most of the forms mentioned above. No special plotting function is needed for plotting, just use plot2d(), plot3d() or draw(). For integration of functions pw.mac converts the function to signum() or between() form and then integrates these functions. Pw.mac cannot differentiate iif() yet. Try iif(x>0, x,-x), this function can be converted to other forms and then differentiated and converted back to iif() form. There are a lot of conversion functions for changing the form of the output from one type to another supplied. I have already shown how to convert iif(x>0, x,-x) to abs(x) and max() to iif().

In the next section we compute a Fourier series of a piecewise continuous function. It is so easy you will hopefully be impressed.

Example Session Follows:

(%i1) (kill(all),now:elapsed_run_time(),showtime:true,load(pw),use_between:true,display2d:false);
Evaluation took 0.5700 seconds (0.5700 elapsed)
(%o0) false

(%i1) (load(draw),wxplot_size:[900,400]);
Evaluation took 0.1100 seconds (0.1100 elapsed)
(%o1) [900,400]

(%i2) declare(n,integer);assume(n>0);
Evaluation took 0.0000 seconds (0.0000 elapsed)
(%o2) done
(%i3) Evaluation took 0.0000 seconds (0.0000 elapsed)
(%o3) [n > 0]

(%i4) f:[-10, -(x^2-25)/15, -5,  -5*sin(3*x), 0, x^2/6, 5, -(8*x-65)/6, 10];  /* piecewise() or its alias pw()
 uses this function list to express a piecewise function */
Evaluation took 0.0000 seconds (0.0000 elapsed)
(%o4) [-10,(25-x^2)/15,-5,-5*sin(3*x),0,x^2/6,5,(65-8*x)/6,10]

(%i5) pw(f,x); /* this creates the piecewise function from the list f */
Evaluation took 0.0000 seconds (0.0000 elapsed)
(%o5) -5*between(x,-5,0)*sin(3*x)+(65-8*x)*between(x,5,10)/6+x^2*between(x,0,5)/6+(25-x^2)*between(x,-10,-5)/15

(%i6) a0 : radcan(pwint(pw(f,x), x,'minf,'inf)/20);
Evaluation took 0.1300 seconds (0.1300 elapsed)
(%o6) -(cos(15)-1)/12

(%i7) bn : (s:(1/10)*pwint(pw(f,x)*sin(n*%pi*x/10), x, 'minf, 'inf),radcan(%%));
Evaluation took 0.0500 seconds (0.0500 elapsed)
(%o7) ((15*%pi^4*n^4-450*%pi^3*n^3)*sin((%pi*n+30)/2)+
(-15*%pi^4*n^4-450*%pi^3*n^3)*sin((%pi*n-30)/2)+(140*%pi^3*n^3-126000*%pi*n)*sin(%pi*n/2)
 +(120*%pi^2*n^2-108000)*cos(%pi*n/2)+(-15*%pi^4*n^4+13580*%pi^2*n^2-72000)*(-1)^n-200*%pi^2*n^2+180000)
 /(6*%pi^5*n^5-5400*%pi^3*n^3)

(%i8) an : (s:(1/10)*pwint(pw(f,x)*cos(n*%pi*x/10), x, 'minf, 'inf), radcan(s));
Evaluation took 0.0500 seconds (0.0500 elapsed)
(%o8) -((15*%pi^4*n^4-450*%pi^3*n^3)*cos((%pi*n+30)/2)+
(-15*%pi^4*n^4-450*%pi^3*n^3)*cos((%pi*n-30)/2)+(280*%pi^2*n^2-252000)*sin(%pi*n/2)
 +(198000*%pi*n-220*%pi^3*n^3)*cos(%pi*n/2)+(160*%pi^3*n^3-144000*%pi*n)*(-1)^n+900*%pi^3*n^3)
  /(6*%pi^5*n^5-5400*%pi^3*n^3)

(%i9) define( b(n), bn)$ define(a(n), an)$ define(h(x), pw(f,x))$
Evaluation took 0.0000 seconds (0.0000 elapsed)
(%i10) Evaluation took 0.0000 seconds (0.0000 elapsed)
(%i11) Evaluation took 0.0000 seconds (0.0000 elapsed)
(%i12) g(x, terms) :=a0+sum(b(n)*sin(n*%pi*x/10)+a(n)*cos(n*%pi*x/10), n, 1, terms)$
Evaluation took 0.0000 seconds (0.0000 elapsed)
(%i13) terms:40;
Evaluation took 0.0000 seconds (0.0000 elapsed)
(%o13) 40
(%i14) draw2d(
        color = black,
        explicit(h(x),x,-10,10),
        yrange=[-8,6],
        color=blue,
        explicit(g(x,terms),x,-10,10),
        label_alignment=right,
        label(["f(x) = pw([-10, -(x^2-25)/15, -5,  -5*sin(3*x), 0, x^2/6, 5, -(8*x-65)/6, 10],x)",9.5, -7]),
        label([concat("Fourier series for n = ", terms), 9.5, -7.5]))$

Output from draw2d()

x Piecewise continuous function and its Fourier Series Approximation

Lastly I have added example code which shows that pwint and diff are self consistent for the piecewise function above.


(%i15) s:radcan(pw(f,x));  /* the original function */
Evaluation took 0.0100 seconds (0.0100 elapsed)
(%o15) -(150*between(x,-5,0)*sin(3*x)+(40*x-325)*between(x,5,10)-5*x^2*between(x,0,5)+(2*x^2-50)*between(x,-10,-5))/30

(%i16) t:radcan(diff(pwint(pw(f,x),x),x));  /* int and then diff and you get the function back */
Evaluation took 0.0300 seconds (0.0300 elapsed)
(%o16) -(150*between(x,-5,0)*sin(3*x)+(40*x-325)*between(x,5,10)-5*x^2*between(x,0,5)+(2*x^2-50)*between(x,-10,-5))/30

(%i17) pwsimp(s-t,x);  /* proof */
(%o17) 0

(%i18) use_pwdelta(true);  /* this is important for the next step */
Evaluation took 0.0000 seconds (0.0000 elapsed)
(%o18) "pwdeltas() will be used"

(%i19) u:radcan(pwint(diff(pw(f,x),x),x));  /* diff and then integrate. This must be the same as pw(f,x) except for a constant.  */
Evaluation took 0.0500 seconds (0.0500 elapsed)
(%o19) (50*signum(x+10)-50*signum(x+5)-300*between(x,-5,0)*sin(3*x)-80*x*between(x,5,10)+
    10*x^2*between(x,0,5)-4*x^2*between(x,-10,-5)+325*signum(x-5)-325*signum(x-10)+300*sin(15)-800)/60

(%i20) pwsimp(s-u,x);  /*  Proof */
Evaluation took 0.0200 seconds (0.0200 elapsed)
(%o20) (300*sin(15)-800)/60

(%i21) elapsed_run_time()-now;  /* this does it */
Evaluation took 0.0000 seconds (0.0000 elapsed)
(%o21) 23.36

You can also do convolution integrals of piecewise continuous functions with pw.mac.

Type:

define(p(x),pw([-a,1/(2*a),a],x)); define(pdf(x), pwint(p(p)*p(x-p),p,minf,inf)); for i : 1 thru n do define(pdf(x), pwint(pdf(p)*p(x-p),p,minf,inf));

This can be quite time consuming for large n or complex piecewise functions and can easily exceed Maxima's memory limits even for low n. Of course the functions integral must be finite from minf to inf. For the case where pwint(p(x),x,minf,inf) = 1 there are some interesting uses in statistical engineering. I have a graphic of a convolution integral of a simple piecewise function for n = 50 (plus 2 for the first define()) and the code to create it is here.

Anyway, try pw.mac, it is free. I hope you find it useful.

--Richhennessy610 12:49, 4 March 2011 (UTC)

Personal tools