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

Main Page

From piecewisefunc

Revision as of 16:35, 6 July 2013 by Richhennessy610 (Talk | contribs)
(diff) ← Older revision | Current revision (diff) | Newer revision → (diff)
Jump to: navigation, search

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