Main Page

From piecewisefunc

(Difference between revisions)
Jump to: navigation, search
m
m
 
(15 intermediate revisions not shown)
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 most recent earlier versions. Pw.mac is a free program and is open source.   
+
Did you know that you can integrate the abs(x) function in one step?  Did you know you can also differentiate it too, also in one step?  Well, you can.  This package, pw.mac, an add-on to Maxima can do it.  No need to break the problem up into steps, pw.mac does it all for you.  Pw.mac is the first software package for [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 one call to 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. One function call does it in most cases with pw.mac.  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 very 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] for convenience, more current help is included in the compressed download file.  The version of pw.mac is now 6.4 and works with Maxima 5.25.1 and most recent earlier versions of Maxima. 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.
Line 10: Line 10:
<pre>
<pre>
-
(%i2) piecewise([minf,-x^2,0,x^2-x,inf],x) /* pw() is an alias for piecewise */;
+
(%i1) piecewise([minf,-x^2,0,x^2-x,inf],x) /* pw() is an alias for piecewise */;
                     2                        2
                     2                        2
                   (x  - x) (signum(x) + 1)  x  (1 - signum(x))
                   (x  - x) (signum(x) + 1)  x  (1 - signum(x))
Line 16: Line 16:
                             2                      2
                             2                      2
 +
(%i2) pwsimp(%,x,'array);  // the most readable form of representation of the piecewise function.
 +
 +
 +
                                                [                                          2  ]
 +
                                                [ If  x  in  (  minf  ,  0  )  then  - x  ]
 +
                                                [                                            ]
 +
                                                [ If  x  in  [  0    ,  0  ]  then    0    ]
 +
                                                [                                            ]
 +
                                                [                                      2    ]
 +
                                                [ If  x  in  (  0    ,  inf  )  then  x  - x ]
 +
 +
                                                                    Done
(%i3) piecewise([minf,-x^2,0,x^2-x,inf],x,'iif);  /* you can simplify this a little with pulliniif() */
(%i3) piecewise([minf,-x^2,0,x^2-x,inf],x,'iif);  /* you can simplify this a little with pulliniif() */
                                                             2                                            2
                                                             2                                            2
Line 31: Line 43:
/* by default piecewise functions are expressed in terms of the signum() and unit_spike() functions but
/* 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.   
   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. */
+
   The "if then else" form is not very useful for further manipulation however it can be turned back into a sum
-
 
+
  with iif2sum(ifthen2iif(%)); */
-
(%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() */
+
(%i5) piecewise([minf, - x^2, 0, x^2  - x, inf], x, 'ifthen);
-
                                                  3    2                   3    2
+
                                                        2                                                            2
-
                                                  x   x                   x   x
+
                          2                            x                            2                              x  - x
-
                                  3            (-- - --) signum(x)   3   -- - --
+
        (if x < 0 then - x  else (if equal(x, 0) then - -- else 0)) + (if x > 0 then x  - x else (if equal(x, 0) then ------ else 0))
-
                                  signum(x)   3   2              x   3    2
+
                                                        2                                                              2
-
                                  ------------ + ------------------- - -- + -------
+
(%i6) ifthen2iif(%);
-
                                        6                  2            6      2
+
                                                                  2                                            2
 +
                                            2                    x                    2                      x  - x
 +
                              iif(x < 0, - x , iif(equal(x, 0), - --, 0)) + iif(x > 0, x - x, iif(equal(x, 0), ------, 0))
 +
                                                                  2                                              2
 +
(%i7) pulliniif(%);
 +
                                                                2        2
 +
                                        2              2      x  - x  x
 +
                          iif(x < 0, - x , iif(x > 0, x  - x, ------ - --))
 +
                                                                2     2
 +
(%i8) iif2ifthen(%);
 +
                                                                        2        2
 +
                                    2                      2          x - x   x
 +
                    if x < 0 then - x  else (if x > 0 then x  - x else ------ - --)
 +
                                                                        2      2
 +
(%i9) ifthen2iif(%);
 +
                                                                2        2
 +
                                        2              2      x  - x  x
 +
                          iif(x < 0, - x , iif(x > 0, x  - x, ------ - --))
 +
                                                                2      2
 +
(%i10) pwsimp(iif2sum(%), x);
 +
                                        2                        2
 +
                                      (x  - x) (signum(x) + 1)   x  (1 - signum(x))
 +
                                      ------------------------ - ------------------
 +
                                                  2                      2
 +
(%i11) signum2abs(pwint(%, x));
 +
                                                2                    2
 +
                                            (4 - 3 x) abs(x) - 3 x
 +
                                              --------------------------
 +
                                                          12
 +
(%i12) ratsimp(signum2abs(pwint(piecewise([minf,-x^2,0,x^2-x,inf],x,'abs) ,x)));
 +
                                                2                    2
 +
                                            (4 x  - 3 x) abs(x) - 3 x
 +
                                              --------------------------
 +
                                                          12
-
(%i7) periodic(sin(a*x),x,0,2*%pi/a); /* you can write periodic piecewise continuous functions using pw.mac */
+
(%i13) periodic(sin(a*x),x,0,2*%pi/a); /* you can write periodic piecewise continuous functions using pw.mac */
                             a x          a x
                             a x          a x
Line 53: Line 93:
                             2 %pi        2 %pi
                             2 %pi        2 %pi
-
(%i8) assume(a>0);
+
(%i14) assume(a>0);
                                           [a > 0]
                                           [a > 0]
-
(%i9) intperiodic(sin(a*x),x,0,2*%pi/a);  /* and integrate them at least once */
+
(%i15) intperiodic(sin(a*x),x,0,2*%pi/a);  /* and integrate them at least once */
                               a x          a x
                               a x          a x
                   cos(2 %pi (----- - floor(-----)))
                   cos(2 %pi (----- - floor(-----)))
Line 63: Line 103:
                                 a
                                 a
-
(%i10) pwint((x+signum(x-a))^5,x);   
+
(%i16) pwint((x+signum(x-a))^5,x);   
                                         3      3                                          6      4      2
                                         3      3                                          6      4      2
           5    5                  10 x    10 a                                          x    5 x    5 x
           5    5                  10 x    10 a                                          x    5 x    5 x
Line 70: Line 110:
-
(%i11) pwsimp((signum(x-4)*signum(x-5))+(signum(x-4)-signum(x-5)-1),x);  /*  you can show two forms are the same  
+
(%i17) 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() */
                                                                             with pwsimp() */
Line 76: Line 116:
-
(%i12) pwint(x^2*pwdelta(x-a),x,b,c);  
+
(%i18) pwint(x^2*pwdelta(x-a),x,b,c);  
       /* pwdelta is a curious animal, similar to the diracdelta() function.  Use with caution.
       /* pwdelta is a curious animal, similar to the diracdelta() function.  Use with caution.
           See the help for more information */
           See the help for more information */
Line 83: Line 123:
                           -----------------------------------
                           -----------------------------------
                                             2
                                             2
-
(%i13) pwint(f(x)*pwdelta(x-a),x,minf,inf);  /* pwdelta picks out f(a) */
+
(%i19) pwint(f(x)*pwdelta(x-a),x,minf,inf);  /* pwdelta picks out f(a) */
                                
                                
                                           f(a)
                                           f(a)
    
    
-
(%i14) pwint(f(x)*pwdelta(x-a),x,a,a);  /* but not in this case.  You can't integrate from a to a and  
+
(%i20) 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 */
                                               get a nonzero result with pw.mac */
                                             0
                                             0
-
(%i15) 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  
+
(%i21) 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 */
then you would get -a^2 */
Line 99: Line 139:
                               a  signum(e)
                               a  signum(e)
-
(%i16) is(equal(pwint(pwdelta(x-a),x),(signum(x-a)+1)/2));  /* true by definition of pwdelta */
+
(%i22) is(equal(pwint(pwdelta(x-a),x),(signum(x-a)+1)/2));  /* true by definition of pwdelta */
                                   true
                                   true
-
(%i17) 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
+
(%i23) 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  */;
                                                                                           involving signum() functions  */;
                                   0
                                   0
-
(%i18) pushout((s-signum(x-s))^5,'signum);  /* factor out the signum's */
+
(%i24) pushout((s-signum(x-s))^5,'signum);  /* factor out the signum's */
                 5                    4              2      3              3      2            4                  5
                 5                    4              2      3              3      2            4                  5
         - signum (x - s) + 5 s signum (x - s) - 10 s  signum (x - s) + 10 s  signum (x - s) - 5 s  signum(x - s) + s
         - signum (x - s) + 5 s signum (x - s) - 10 s  signum (x - s) + 10 s  signum (x - s) - 5 s  signum(x - s) + s
-
(%i19) linearize(%); /* signum(s)^2=1 unless s = 0, so assume s # 0, This is okay in case the function is continuous. */  
+
(%i25) 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
                         4                    2                                  5      3
                   - 5 s  signum(x - s) - 10 s  signum(x - s) - signum(x - s) + s  + 10 s  + 5 s
                   - 5 s  signum(x - s) - 10 s  signum(x - s) - signum(x - s) + s  + 10 s  + 5 s
-
(%i20) pushout(%,'signum);  /* factor out the signum */
+
(%i26) pushout(%,'signum);  /* factor out the signum */
                                 4      2                      5      3
                                 4      2                      5      3
                           (- 5 s  - 10 s  - 1) signum(x - s) + s  + 10 s  + 5 s
                           (- 5 s  - 10 s  - 1) signum(x - s) + s  + 10 s  + 5 s
-
(%i21) sin(abs(x-a))*sin(abs(x-b));  /* you can convert */
+
(%i27) sin(abs(x-a))*sin(abs(x-b));  /* you can convert */
                   sin(abs(x - a)) sin(abs(x - b))
                   sin(abs(x - a)) sin(abs(x - b))
-
(%i23) abs2iif(%); /* to iif() form */
+
(%i28) abs2iif(%); /* to iif() form */
                 sin(iif(x - a > 0, x - a, a - x)) sin(iif(x - b > 0, x - b, b - x))
                 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() */
+
(%i29) pulliniif(%);  /* do this to pull the abs()'s inside the iif() instead of inside the sin() */
       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),  
       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)))
                                                     sin(x - a) sin(x - b)))
-
(%i25) signum2abs(simpspikes(iif2sum(iif(x>a,x-a,-x+a)))); /* here is how to convert from an iif() form to an abs() form */
+
(%i30) signum2abs(simpspikes(iif2sum(iif(x>a,x-a,-x+a)))); /* here is how to convert from an iif() form to an abs() form */
                                                             abs(x-a)
                                                             abs(x-a)
-
(%i26) pulliniif(abs2iif(max2abs(max(1,x))));  /* max() and min() are not supported forms for expressing piecewise functions
+
(%i31) 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 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 */
           but true. Max2abs() and min2abs() are Barton Willis' idea and I use them with his permission in pw.mac */
Line 144: Line 184:
                                                         iif(x - 1 > 0, x, 1)
                                                         iif(x - 1 > 0, x, 1)
-
(%i27) pulliniif(abs2iif(min2abs(min(x^2,x))));  /* min() to iif() */
+
(%i32) pulliniif(abs2iif(min2abs(min(x^2,x))));  /* min() to iif() */
                                                               2              2
                                                               2              2
                                                         iif(x  - x > 0, x, x )
                                                         iif(x  - x > 0, x, x )
-
(%i28) pwsimp(%o2,x,'ifthen);      /*  you can use pwsimp to convert from one form to another equivalent form in most cases */
+
(%i33) 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
   - signum(x) - unit_spike(x) + 1            2                                  signum(x) - unit_spike(x) + 1          2
Line 155: Line 195:
                   2                                                                            2
                   2                                                                            2
-
(%i29) pwsimp(%o2,x,'abs);  /* abs() form */
+
(%i34) pwsimp(%o2,x,'abs);  /* abs() form */
                                       (2 x - 1) abs(x) - x
                                       (2 x - 1) abs(x) - x
Line 161: Line 201:
                                               2
                                               2
-
(%i30) pwsimp(%o2,x,'iif);  /* iif() form */
+
(%i35) pwsimp(%o2,x,'iif);  /* iif() form */
                                                     2                  2
                                                     2                  2
                                       iif(x < 0, - x , 0) + iif(x > 0, x  - x, 0)
                                       iif(x < 0, - x , 0) + iif(x > 0, x  - x, 0)
-
(%i31) plot2d([pw([minf,x^2,0,-x^2,inf],x)], [x,-4,4])$ /* this is one inefficient way to plot a piecewise function */
+
(%i36) plot2d([pw([minf,x^2,0,-x^2,inf],x)], [x,-4,4])$ /* this is one inefficient way to plot a piecewise function */
-
(%i32) h(x):=x^2*between(x,0,inf)-x^2*between(x,minf,0);  /* this is not very efficient for plotting either */
+
(%i37) h(x):=x^2*between(x,0,inf)-x^2*between(x,minf,0);  /* this is not very efficient for plotting either */
                                               2                      2
                                               2                      2
                                       h(x) := x  between(x, 0, inf) - x  between(x, minf, 0)
                                       h(x) := x  between(x, 0, inf) - x  between(x, minf, 0)
-
(%i33) plot2d([h(x)],[x,-4,4]);
+
(%i38) plot2d([h(x)],[x,-4,4]);
-
(%i34) h(x):=''(x^2*between(x,0,inf)-x^2*between(x,minf,0));  /* this is much better */
+
(%i39) h(x):=''(x^2*between(x,0,inf)-x^2*between(x,minf,0));  /* this is much better */
                                               2                    2
                                               2                    2
                                               x  (signum(x) + 1)  x  (1 - signum(x))
                                               x  (signum(x) + 1)  x  (1 - signum(x))
                                       h(x) := ------------------ - ------------------
                                       h(x) := ------------------ - ------------------
                                                       2                    2
                                                       2                    2
-
(%i35) plot2d([h(x)],[x,-4,4]);
+
(%i40) plot2d([h(x)],[x,-4,4]);
</pre>
</pre>
Output of plot2d:
Output of plot2d:
Line 185: Line 225:
<pre>
<pre>
-
(%i36) 1/(1+sqrt(1-2*x+x^2)); /* this function is piecewise continuous if you assume sqrt(x^2) = abs(x) */
+
(%i41) 1/(1+sqrt(1-2*x+x^2)); /* this function is piecewise continuous if you assume sqrt(x^2) = abs(x) */
-
(%i37) pwint(%,x);  /* pwint recognizes this fact */
+
(%i42) pwint(%,x);  /* pwint recognizes this fact */
                                         signum(x - 1) log(x)  log(x)  log(2 - x) signum(x - 1)  log(2 - x)
                                         signum(x - 1) log(x)  log(x)  log(2 - x) signum(x - 1)  log(2 - x)
Line 193: Line 233:
                                                   2              2                2                  2
                                                   2              2                2                  2
-
(%38) pwint(sum(x*abs(a*x+b[i]), i, 1, n),x);  /* pwint automatically distributes over sums */
+
(%i43) pwint(sum(x*abs(a*x+b[i]), i, 1, n),x);  /* pwint automatically distributes over sums */
                                           n          3        2    3
                                           n          3        2    3
Line 204: Line 244:
                   ^
                   ^
-
(%i39) 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),  
+
(%i44) 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 */
                                         unit_pulse can be integrated directly */
         3            x              x      3              x      2              x      2
         3            x              x      3              x      2              x      2
Line 217: Line 257:
                                                                       12                  12                    5        6
                                                                       12                  12                    5        6
-
(%i40)  pwint(prod(unit_step(x-a[i]),i, 1, 10),x);  /* I tried this in Maple 14, be prepared to wait a long time */  
+
(%i45)  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  ))
             (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
                       1  2  3  4  5  6  7  8  9  10                  1  2  3  4  5  6  7  8  9  10    x
-
(%o40)      ----------------------------------------------------------------------------------------------------------- + -
+
(%o46)      ----------------------------------------------------------------------------------------------------------- + -
                                                                   2                                                        2
                                                                   2                                                        2
</pre>
</pre>
Line 340: Line 380:
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 [http://mysite.verizon.net/res11w2yb/Gaussian_or_Binomial.html a simple piecewise function] for n = 50 (plus 2 for the first define()) and the code to create it [http://mysite.verizon.net/res11w2yb/convolution_animation.mac.txt is here].   
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 [http://mysite.verizon.net/res11w2yb/Gaussian_or_Binomial.html a simple piecewise function] for n = 50 (plus 2 for the first define()) and the code to create it [http://mysite.verizon.net/res11w2yb/convolution_animation.mac.txt is here].   
-
Anyway, try pw.mac, it is free. I hope you find it useful.     
+
So what is the indefinite integral of abs(x)?  The answer is x*abs(x)/2 when expressed in terms the abs(x).  Pwint(abs(x),x) returns x^2*signum(x)/2, which is also correct but expressed in terms of signum(x).  Anyway, try pw.mac, it is free and I hope you find it useful.     
--[[User:Richhennessy610|Richhennessy610]] 12:49, 4 March 2011 (UTC)
--[[User:Richhennessy610|Richhennessy610]] 12:49, 4 March 2011 (UTC)

Current revision as of 16:35, 6 July 2013

Pw.mac

Did you know that you can integrate the abs(x) function in one step? Did you know you can also differentiate it too, also in one step? Well, you can. This package, pw.mac, an add-on to Maxima can do it. No need to break the problem up into steps, pw.mac does it all for you. Pw.mac is the first software package for 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 one call to 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. One function call does it in most cases with pw.mac. 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 very 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 for convenience, more current help is included in the compressed download file. The version of pw.mac is now 6.4 and works with Maxima 5.25.1 and most recent earlier versions of Maxima. 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.

(%i1) 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

(%i2) pwsimp(%,x,'array);  // the most readable form of representation of the piecewise function.


                                                [                                          2  ]
                                                [ If  x  in  (  minf  ,   0   )  then   - x   ]
                                                [                                             ]
                                                [ If  x  in  [   0    ,   0   ]  then    0    ]
                                                [                                             ]
                                                [                                       2     ]
                                                [ If  x  in  (   0    ,  inf  )  then  x  - x ]

                                                                     Done
(%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 it can be turned back into a sum 
   with iif2sum(ifthen2iif(%)); */

(%i5) piecewise([minf, - x^2, 0, x^2  - x, inf], x, 'ifthen);
                                                         2                                                             2
                          2                             x                             2                               x  - x
        (if x < 0 then - x  else (if equal(x, 0) then - -- else 0)) + (if x > 0 then x  - x else (if equal(x, 0) then ------ else 0))
                                                        2                                                               2
(%i6) ifthen2iif(%);
                                                                   2                                             2
                                            2                     x                     2                       x  - x
                              iif(x < 0, - x , iif(equal(x, 0), - --, 0)) + iif(x > 0, x  - x, iif(equal(x, 0), ------, 0))
                                                                  2                                               2
(%i7) pulliniif(%);
                                                                2        2
                                         2              2      x  - x   x
                           iif(x < 0, - x , iif(x > 0, x  - x, ------ - --))
                                                                 2      2
(%i8) iif2ifthen(%);
                                                                        2        2
                                     2                      2          x  - x   x
                    if x < 0 then - x  else (if x > 0 then x  - x else ------ - --)
                                                                         2      2
(%i9) ifthen2iif(%);
                                                                2        2
                                         2              2      x  - x   x
                           iif(x < 0, - x , iif(x > 0, x  - x, ------ - --))
                                                                 2      2
(%i10) pwsimp(iif2sum(%), x);
                                         2                         2
                                       (x  - x) (signum(x) + 1)   x  (1 - signum(x))
                                       ------------------------ - ------------------
                                                  2                       2
(%i11) signum2abs(pwint(%, x));
                                                 2                    2
                                             (4 x  - 3 x) abs(x) - 3 x
                                              --------------------------
                                                           12
(%i12) ratsimp(signum2abs(pwint(piecewise([minf,-x^2,0,x^2-x,inf],x,'abs) ,x)));
                                                 2                    2
                                             (4 x  - 3 x) abs(x) - 3 x
                                              --------------------------
                                                           12

(%i13) 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

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

(%i15) 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

(%i16) 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


(%i17) 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


(%i18) 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
(%i19) pwint(f(x)*pwdelta(x-a),x,minf,inf);  /* pwdelta picks out f(a) */

                               
                                           f(a)
  
(%i20) 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

(%i21) 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)

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

                                   true

(%i23) 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  */;

                                   0

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

(%i25) 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
                   - 5 s  signum(x - s) - 10 s  signum(x - s) - signum(x - s) + s  + 10 s  + 5 s

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

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

(%i27) sin(abs(x-a))*sin(abs(x-b));  /* you can convert */

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

(%i28) abs2iif(%); /* to iif() form */

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

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

       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)))

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

                                                            abs(x-a)

(%i31) 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 */

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

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

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

(%i33) 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

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

                                      (2 x - 1) abs(x) - x
                                      --------------------
                                               2

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

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

(%i37) 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)
(%i38) plot2d([h(x)],[x,-4,4]);

(%i39) 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
(%i40) plot2d([h(x)],[x,-4,4]);

Output of plot2d:

x Plot of h(x)


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

(%i42) 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

(%i43) 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
                                           >    (---------------- - ----) signum(a x + b )
                                          /             6              2                i
                                          ====                      6 a
                                          i = 1
                  ^

(%i44) 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
       ---------------- - -------------------- + -------------------- + -------------------- - 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

(%i45)  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
(%o46)       ----------------------------------------------------------------------------------------------------------- + -
                                                                  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.

So what is the indefinite integral of abs(x)? The answer is x*abs(x)/2 when expressed in terms the abs(x). Pwint(abs(x),x) returns x^2*signum(x)/2, which is also correct but expressed in terms of signum(x). Anyway, try pw.mac, it is free and I hope you find it useful.

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

Personal tools