MathGroup Archive 2003

[Date Index] [Thread Index] [Author Index]

Search the Archive

Re: Integrate with piecewise function

Daniel Lichtblau <danl at> wrote in message news:<bn2iko$o0d$1 at>...
> Other approaches, already mentioned in this thread, include rewriting
> the function in terms of Abs or UnitStep.
Among 'other approaches', here is my own approach :

To rewrite piecewise constant functions or
piecewise linear functions in terms of UnitSteps
I wrote the two modules included below.

I hope my little DIY will be useful for those who want
to replace the Ifs and Whichs by UnitSteps.


pieceConstant::usage = "pieceConstant[y0, piecePoints, x] \n
returns a piecewise constant function expressed in terms of UnitSteps.\n
arguments : \n
   y0 = leftmost value  \n
   piecePoints = {{x1,y1},...,{xn,yn}} where yi are right values \n
   at points of discontinuity \n
   x = variable used to generate output expression \n
example : \n
   pieceConstant[0,{{0,1},{1,0}},x] returns the square 'bump' function \n
    UnitStep[x] - UnitStep[x-1] ";

Module[{lp = Length[piecePoints],f,a,b,eq,vars},
f[u_] := Sum[a[i]UnitStep[u-piecePoints[[i,1]]],{i,1,lp}]+b;
eq = Append[f[#[[1]]] == #[[2]]& /@ piecePoints, f[-Infinity] == y0]; 
vars = Append[Table[a[i],{i,1,lp}],b]; 

pieceLinear::usage="pieceLinear[leftSlope0, piecePoints, x] \n
returns a piecewise linear function expressed in terms of UnitSteps \n
(function may be discontinuous).\n
arguments : \n
   leftSlope0 : leftmost slope \n
   piecePoints format : \n
   {{x1, lefty1, righty1, rightSlope1},..., \n
    {xn, leftyn, rightyn, rightSlopen}} \n
   x : variable \n
example : \n
   pieceLinear[0,{{0,0,0,1},{1,1,0,1},{2,1,0,1},{3,1,0,0}},x] \n
   returns the 'sawtooth' function \n
   x UnitStep[x] - (x-2)UnitStep[x-3] - UnitStep[x-2] - UnitStep[x-1]";

pieceLinear[leftSlope0_, piecePoints_List, x_]:=
Module[{lp = Length[piecePoints],
(* use delta functions to deal with discontinuities : *)
deltaFunctions = (#[[3]]-#[[2]])DiracDelta[x-#[[1]]]& /@ piecePoints;
pieceConstantPoints = {#[[1]],#[[4]]}& /@ piecePoints;
(* build the derivative, then integrate : *)
pieceDerivative = pieceConstant[leftSlope0, pieceConstantPoints,x]+
pieceLinearFunction = Integrate[pieceDerivative,x]+k;
(* compute the integration constant k : *)
sol=If[lp == 1,
(* specific case if there is only one point : *)
Solve[(pieceLinearFunction/.x -> piecePoints[[1,1]]) == 
(* generic case : *)
Solve[(pieceLinearFunction /. x ->
(piecePoints[[1,1]] + piecePoints[[2,1]])/2) ==
(piecePoints[[1,3]] + piecePoints[[2,2]])/2,k]];
pieceLinearFunction /. sol[[1]]

  • Prev by Date: LogLogListPlot, LogLinearListPlot- Ticks and TickMarks
  • Next by Date: Re: recode procedural algorithm to faster functional module
  • Previous by thread: Re: Integrate with piecewise function
  • Next by thread: non linear least square minimization....