Integrate[] and InterpolatingFunction's

*To*: mathgroup at yoda.physics.unc.edu*Subject*: Integrate[] and InterpolatingFunction's*From*: TDR at vaxc.cc.monash.edu.au*Date*: 01 Jul 1993 22:56:02 +1100

I've just re-read the posting by Paul Abbott which contains a description of the internal structure of 1D InterpolatingFunction's. >From it I have worked out how (see below) to integrate them. One useful thing that would require integrating InterpolatingFunction's is computing the time-average of a solution generated by NDSolve. I'm sure there are other useful applications too. Anyway here's a sample session that uses the extended Integrate[] to integrate an InterpolatingFunction. See line In[5] which does the integration; and also line In[7] which checks that it differentiates back correctly. In[3]:= NDSolve[{y''[t] == -Sin[y[t]], y[0]==Pi/4, y'[0]==0}, y, {t,0,9}][[1]] Out[3]= {y -> InterpolatingFunction[{0., 9.}, <>]} In[4]:= solution = y[t] /. % Out[4]= InterpolatingFunction[{0., 9.}, <>][t] In[5]:= integral = Integrate[solution,t] Out[5]= InterpolatingFunction[{0., 9.}, <>][t] In[6]:= derivative = D[solution,t] Out[6]= InterpolatingFunction[{0., 9.}, <>][t] In[7]:= D[integral,t] == solution Out[7]= True In[8]:= Plot[{solution,integral,derivative}, {t,0,9}] Out[8]= -Graphics- The code I used to extend Integrate[] is given below. It's just my first attempt to prove that it can be done - so the code is not very elegant. Begin["`private`"]; Unprotect[Integrate]; Integrate[f_InterpolatingFunction[x_], x_] := int[f][x]; Protect[Integrate]; (** CONST is passed globally. This is bad style I know. **) int[f_] := Block[{intf=f, CONST=0}, intf[[2 ]] = Map[doit, intf[[2]]]; intf[[2,1,3]] = {0,0}; intf[[2,1,4]] = {0}; intf ]; (* the NDSolve[] case : all coefs are zero *) doit[{ta_, tb_, d_, c:{(0)..}}] := Block[{d1=d/Range[Length[d]], dt=ta-tb}, CONST += dt*Fold[(#2-#1*dt)&, Last[d1], Rest[Reverse[d1]]]; {ta, tb, Prepend[d1,CONST], Prepend[c,0]} ]; (* the Interpolation[] case : coefs are non-zero *) doit[{ta_, tb_, d_, c_}] := "to be implemented"; End[]; I have formulas for the non-implemented case above. The formulas are resonably simple, but I can't think of any easy way to code them in Mathematica. I'll send the formulas if anyone wants to know them. Terry Robb