MathGroup Archive 1993

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

Search the Archive

Integrate[] and InterpolatingFunction's

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.

      Integrate[f_InterpolatingFunction[x_], x_] := int[f][x];

    (** 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};

    (* 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";

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
  • Next by thread: accelerated degradation test data analysis