Maxima - differentiating a piecewise function

Question

Suppose you have a function defined by intervals, such as

f(x):=block(if x<0 then x^2 else x^3);

When we differentiate it with

diff(f(x),x);

we get

d/dx (if x<0 then x^2 else x^3)

whereas I'd like to get

(if x<0 then 2*x else 3*x^2)

Is there a way to obtain such result?

Answer 1

This may help in a simple case:

(%i1) f(x):= charfun(x<0)*x^2 + charfun(x>=0)*x^3$

(%i2) gradef(charfun(y), 0)$

(%i3) diff(f(x),x);
                                           2
(%o3)              2 x charfun(x < 0) + 3 x  charfun(x >= 0)

charfun, gradef

You can try also Pw.mac package from Richard Hennessy.

Answer 2

Here's a different approach using a simplification rule for "if" expressions. The unsolved part here is to detect discontinuities and generate delta functions for those locations. If you want to ignore those, you can define FOO to return 0. Note that I didn't attempt to implement the function discontinuities; that part is unsolved here. I can give it a try if there is interest.

(%i1) display2d : false $
(%i2) matchdeclare ([aa, bb, cc], all, xx, symbolp) $
(%i3) 'diff (if aa then bb else cc, xx) $
(%i4) tellsimpafter (''%, apply ("if", [aa, diff (bb, xx), true, diff (cc, xx)]) + FOO (aa, bb, cc, xx)) $
(%i5) FOO (a, b, c, x) := 'lsum ((ev (c, x = d) - ev (b, x = d)) * delta (d, x), d, discontinuities (a, x)) $
(%i6) diff (if x > 0 then x^2 else x^3, x);
(%o6) (if x > 0 then 2*x else 3*x^2)+'lsum((d^3-d^2)*delta(d,x),d,
                                           discontinuities(x > 0,x))

Answer 3

Building on slitinov's answer I wrote this quite naive implementation for functions with more than two "pieces":

gradef(charfun(dummy),0)$

/* piecewise function definition */
itv: [[x<0],[x>=0,x<1], [x>=1]]; /* intervals */
fi:  [ 1,    x^2+1,      2*x  ]; /* local functions */

/* creation of global function f and its derivative df */
f:0;
for i: 1 thru 3 do f:f+charfun(apply("and",itv[i]))*fi[i];
df:diff(f,x);

/* display of local functions and derivatives */
for i: 1 thru 3 do (
  apply(assume,itv[i]),
  newline(),
  print(itv[i]),
  print("f = ",ev(f)),
  print("df = ",ev(df)),
  apply(forget,itv[i])
  );

plot2d([f,df],[x,-2,3],[y,-1,5],[style,[lines,4,3],[lines,2,2]]);