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