windowed functions, integration errors, and numerical stability

On Fri, Apr 27, 2007 at 06:10:38PM +0900, Valery Pipin wrote:
> I use
> phi(x):= .5*(1.+erf(a*(x-xb)));
> with a and xb as parameters to control the stepnees
> also (tanh((x - xb)*a) + 1.) is good

The sigmoid, or tanh based function above is extremely convenient
analytically, but produces problems numerically.

For example modelling the US tax code (married filing jointly)

turnon(var,a) := (tanh(20*(var-a))+1)/2;

taxrate(x) := turnon(x,0)*.1+turnon(x,15650)*.05+turnon(x,63700)*.1+turnon(x,128500)*.03+turnon(x,195850)*.05+turnon(x,349700)*.02;

if we try to integrate taxrate to get the amount of taxes paid when
your income is x... 

first off the maxima integrator can't do it... 

If you force it to do each term in the sum seperately it will handle
it fine (ie. map(lambda([a],integrate(a,x)),taxrate(x))) ... but it
will produce terms which are numerically no good, specifically it
produces log(cosh(*bignumbers*)) which triggers cosh(*bignumber*) to
overflow, even though log(cosh(*bignumber*)) is perfectly fine.

I can think of several methods to deal with this, but none of them are
ideal. One example is to use tellsimp and friends to replace the
numerically unstable expression with yet another approximation...

Another is to give up on analytical solutions and use the step
functions or piecewise splines with numerical integration...

Anyone else tried doing something like this? Any suggestions?

-- 
Daniel Lakeland
dlakelan at street-artists.org
http://www.street-artists.org/~dlakelan