BeginPackage["Piecewise`"]; (*Piecewise Package: "Defines operations for manipulating piecewise defined functions of the form f=PW[var, func_List, knot_List]; where Length[knot]=Length[func]+1 and func is with respect to the independent variable var. The header PW[] declares f of type piecewise function. Example: f=PW[x,{1+x,1-x},{-1,0,1}]; This is a preliminary package. First version written by Doug Hardin and David Roach, 3/15/98" Changes to plot: admits usual Plot options and plots several functions together. Also added trans and dil. Changes by Doug Hardin 2/16/99 Changed format of PW function (pw is now capital and there is no silly extra layer of brackets:Old format was pw[{x,{1+x,1-x},{-1,0,1}}]. Changes by Doug Hardin 3/18/99 To do: Add f[pwf_PW] without messing up plot. Add pwf_PW[g] for monotone g. Add convolution and basic spline functions. Make affine transformation Listable. Changes by Doug Hardin 10/1/99 Added getInt, eval, and overloaded Power function. To Do: rethink compositions. 1) pwf[mf] where mf is a monotone function. This would include rewriting trans, aff, dil in this form. plot will be a problem as well. Also will probably rename plot PlotPW. 2) f[pwf] where f is an arbitrary function. 3) Longer term, include symbolic knots. *) PW::usage= "Head for piecewise defined functions: Example:PW[x,{1-x,1+x},{-1,0,1}]. Built-in addition, multiplication, translation and dilation. Addition: pwf1+pwf2. Multiplicaton: pwf1*pwf2. Translation: trans[k][pwf] translates k units. Dilation: dilate[j][pwf] dilates by j. Trans & Dilation: aff[j,k][pwf] or pwf[j,k]. Note: The piecewise defined functions are assumed zero outside of the knot range. -Infinity and Infinity are acceptable knots." int::usage="Integrates piecewise defined functions over the reals. int[pwf];" ip::usage="Calculates the L2 innner product of two piecewise defined functions. ip[pwf1,pwf2];" plot::usage= "Plots piecewise defined functions over their support. Note: Infinitely supported functions are truncated. plot[pwf]; plot[pwf,opts]-produces graphics with no output." antiD::usage="antiD[pwf] gives the antiderivative: Integrate[f[t],{t,a,x}] where the a is smallest knot." Refine::usage="Refines a function given a superset of knots: Refine[pwf,{k1,k2,...,kn}]." compose::usage="compose[f_,pwf_PW]=f[pwf]" knotvals::usage="Calculates function values at the knots. knotvals[pwf];" trans::usage="trans[b][pwf] is composition of a PW function pwf with x->x-b." dilate::usage="dilate[a][pwf] is composition of a PW function pwf with x->ax." aff::usage="aff[a,b][pwf]is compostion of a PW function pwf with x->ax-b. Alternate form: pwf[a,b]." getInt::usage="getInt[a,{b_1,..., b_n] finds the value of k so that b_k-1 < a <= b_k. Note b_0=-Infinity and b_n+1=Infinity." eval::usage="eval[pwf][a] evaluates a PW function pwf at a. Alternate form: pwf[a]. This is also listable." Begin["`Private`"] pwf_PW*pwg_PW^:=mult[pwf,pwg]; c_ * pwf_PW^:=MapAt[c*#&,pwf,2]; pwf_PW+pwg_PW^:=add[pwf,pwg]; Power[f_PW,n_]^:=MapAt[Power[#,n]&,f,2] f_PW[a_]:=eval[f][a]; trans[b_][pwf_PW]:=PW[pwf[[1]], pwf[[2]]/.pwf[[1]]->pwf[[1]]-b, (pwf[[3]]+b)]//Simplify; dilate[a_][pwf_PW]:=PW[pwf[[1]], pwf[[2]]/.pwf[[1]]->a*pwf[[1]], pwf[[3]]/a]//Simplify; aff[a_,b_][pwf_PW]:=dilate[a][trans[b][pwf]]; trans[b_][pwfL_List]:=Map[trans[b][#]&,pwfL]; dilate[a_][pwfL_List]:=Map[dilate[a][#]&,pwfL]; aff[a_,b_][pwfL_List]:=Map[aff[a,b][#]&,pwfL]; pwf_PW[a_,b_]:=aff[a,b][pwf]; int[pwf_PW]:= Sum[Integrate[pwf[[2,i]],{pwf[[1]],pwf[[3,i]],pwf[[3,i+1]]}], {i,1,Length[pwf[[2]]]}]; plot[pwf_PW,opts___]:=Block[{tbl,var,vec,pts,x}, var=pwf[[1]]; vec=pwf[[2]]; pts=pwf[[3]]; If[pts[[1]]==-Infinity&&pts[[2]]==Infinity,pts={-5,5}]; If[pts[[-1]]==Infinity,pts[[-1]]=pts[[-2]]+1]; If[pts[[1]]==-Infinity,pts[[1]]=pts[[2]]-1]; vec=(vec/.var->x ); tbl=Table[ Plot[vec[[i]],{x,pts[[i]],pts[[i+1]]}, DisplayFunction->Identity, PlotRange->All], {i,1,Length[vec]}]; Show[tbl,opts,DisplayFunction->$DisplayFunction]]; plot[pwfs_List,opts___]:= Show[ Table[plot[ pwfs[[i]],DisplayFunction->Identity ],{i,1,Length[pwfs]}], opts, DisplayFunction->$DisplayFunction] D[pwf_PW]^:=PW[ pwf[[1]], D[pwf[[2]],pwf[[1]]], pwf[[3]] ]; antiD[pwf_PW]:= Block[{prev,temp,vec,pts,var}, var=pwf[[1]]; vec=pwf[[2]]; pts=pwf[[3]]; prev=0; PW[var,Table[temp=Integrate[vec[[i]],{var,pts[[i]],var}]; temp+=prev; prev=(temp/.var->pts[[i+1]]); temp,{i,1,Length[vec]}],pts] ]; ip[pwf_PW,pwg_PW]:=int[pwf*pwg]; mult[pwfunc1_PW,pwfunc2_PW]:=Block[{pwf1,pwf2}, {pwf1,pwf2}=JointKnots[pwfunc1,pwfunc2]; PW[pwf1[[1]],pwf1[[2]]*(pwf2[[2]]/.pwf2[[1]]->pwf1[[1]]),pwf1[[3]]]]; add[pwfunc1_PW,pwfunc2_PW]:=Block[{pwf1,pwf2}, {pwf1,pwf2}=JointKnots[pwfunc1,pwfunc2]; PW[pwf1[[1]],pwf1[[2]]+(pwf2[[2]]/.pwf2[[1]]->pwf1[[1]]),pwf1[[3]]]]; divide[pwfunc1_PW,pwfunc2_PW]:=Block[{pwf1,pwf2}, {pwf1,pwf2}=JointKnots[pwfunc1,pwfunc2]; PW[pwf1[[1]],pwf1[[2]]/(pwf2[[2]]/.pwf2[[1]]->pwf1[[1]]),pwf1[[3]]]]; JointKnots[pwfunc1_,pwfunc2_]:=Block[{knots1,knots2,newknots,posSym1,posSym2}, knots1=pwfunc1[[3]]; knots2=pwfunc2[[3]]; newknots=Union[Select[knots1,NumericQ],Select[knots2,NumericQ]]; posSym1=SymbolPosition[knots1]; posSym2=SymbolPosition[knots2]; Do[ j=posSym1[[i]]; k=pos[{knots1[[j-1]]},newknots][[1]]; newknots=Insert[newknots,knots1[[j]],k+1], {i,1,Length[posSym1]}]; Do[ j=posSym2[[i]]; If[Not[MemberQ[knots1,knots2[[j]]]], k=pos[{knots2[[j-1]]},newknots][[1]]; newknots=Insert[newknots,knots2[[j]],k+1]], {i,1,Length[posSym2]}]; {Refine[pwfunc1,newknots],Refine[pwfunc2,newknots]} ]; Refine[pwfunc_,newknots_]:=Block[{f,fnew,y1,y2,var,i,p}, y1=pwfunc[[3]]; f=pwfunc[[2]]; var=pwfunc[[1]]; p=pos[y1,newknots]; fnew=Join[Table[0,{p[[1]]-1}], Flatten[Table[Table[f[[i]],{p[[i+1]]-p[[i]]}],{i,1,Length[p]-1}]], Table[0,{Length[newknots]-Last[p]}]]; PW[var,fnew,newknots]]; subsetQ[x_,y_]:=Apply[And,Map[MemberQ[y,#]&,x]]; pos[x_,y_]:=Flatten[Map[Position[y,#][[1]]&,x]]; SymbolPosition[x_]:=pos[Select[x,Not[NumericQ[#]]&],x] ourN[x_]:=(N[x]/.(0->0.0)); knotvals[pwf_PW]:=Prepend[Table[pwf[[2,i]]/.pwf[[1]]->pwf[[3, i+1]], {i,1,Length[pwf[[2]]]}],pwf[[2,1]]/.pwf[[1]]->pwf[[3,1]]]; compose[f_,pwf_PW]:=PW[pwf[[1]],Map[f,pwf[[2]]],pwf[[3]]]; getInt[a_,b_List]:= Block[{num}, num=1; While[ num<=Length[b]&&a>=b[[num]],++num]; num-1]; eval[f_PW][a_]:=Block[{intNum}, intNum=getInt[a,f[[3]]]; If[0a), 0]] End[]; EndPackage[];