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

```

• Prev by Date: InterpolatingFunction's and NIntegrate
• Next by Date: accelerated degradation test data analysis
• Previous by thread: InterpolatingFunction's and NIntegrate