Re: Merging InterpolatingFunctions
- To: mathgroup at smc.vnet.net
- Subject: [mg9809] Re: Merging InterpolatingFunctions
- From: Robert Knapp <rknapp>
- Date: Fri, 28 Nov 1997 05:35:21 -0500
- Organization: Wolfram Research, Inc.
- Sender: owner-wri-mathgroup at wolfram.com
Christian Zemlin wrote:
>
> Hi,
>
> I have a set of InterpolatingFunctions, such that the endpoint of each
> is the starting point of the next. Is there any way to merge them into
> one InterpolationFunction? I tried myself and produced the following
> program:
>
> Merge[f1_,f2_] := Module[{Start,Stop}, Start=Extract[f1,{1,1,1}];
> Stop = Extract[f2,{1,1,2}];
> Part1 = Extract[f1,2];
> Part21=Extract[f1,{3,1}];
> Part22=Delete[Extract[f2,{3,1}],1]; Part31= Extract[f1,4];
> Part32=Delete[Extract[f2,4],1];
> InterpolatingFunction[{{Start,Stop}},Part1,{Join[Part21,Part22]},
> Join[Part31, Part32]]
> ];
>
> I guessed the meaning of the entries in an InterpolatingFunction from
> looking at examples, and the program does work in many cases. But
> sometimes it produces an expression that looks to me like an
> InterpolationFunction, but is not treated as one by Mathematica.
> Unfortunately, I did not find a full documentation of
> InterpolatingFunction Objects.
> Can anybody help me with this?
>
The internals of InterpolatingFunction objects are not documented for
the simple reason that most users should treat them as just that, data
objects which have limited access to the data. There are several
contructors for this data object: Interpolation and ListInterpolation
will suit your purpose fine. As to getting data out of the object,
evaluation at a point is the best way.
The internal form of an InterpolatingFunction is subject to change from
version to version. This is not done without good reason. but it does
happen. For example, they changed substantially between V2.2 and V3.0
so that multidimensional InterpolatingFunctions could be differentiated
and integrated. In the next versionto be released there are some
slight changes to take better advantage of mroe compact data
structures.
However, if you want to reconstruct without doing any approximation, you
need to know two things. First, at what points is the interpolation
exact? Second, how many derivatives were specified at each point? It
is not hard to get these two pieces of information from the function.
Below, I will show you an object oriented approach which can easily be
modified if the internal data structure does change:
First, lets make 3 contiguous objects to works with
In[1]:=
if1 = x /. First[NDSolve[{x'[t] == x[t],x[0] == 1},x,{t,0,1}]]
Out[1]=
InterpolatingFunction[{{0.,1.}},"<>"]
In[2]:=
if2 = x /. First[NDSolve[{x'[t] == -if1[t-1]x[t],x[1] ==
if1[1]},x,{t,1,2}]]
Out[2]=
InterpolatingFunction[{{1.,2.}},"<>"]
In[3]:=
if3 = x /. First[NDSolve[{x'[t] == -if2[t-1]x[t],x[2] ==
if2[2]},x,{t,2,3}]]
Out[3]=
InterpolatingFunction[{{2.,3.}},"<>"]
Define a function "Knots" which will give the interpolation points. I
have restricted the definition to 1-dimensionsal objects
In[4]:=
Knots[ifun_InterpolatingFunction] := ifun[[3,1]] /; Length[ifun[[1]]] ==
1
In[5]:=
Knots[if1]
Out[5]=
{0.,0.000894427,0.00178885,0.0223764,0.042964,0.0635516,0.122802,0.182053,
0.241303,0.300553,0.399063,0.497572,0.596081,0.69459,0.793099,0.945494,1.}
Here is another function which gives the number of derivatives which are
given at each knot. 0 means just the function value.
In[6]:=
InterpolatedDerivatives[
ifun_InterpolatingFunction] := ((Length /@ ifun[[4]])- 1) /;
Length[ifun[[1]]] == 1
In[7]:=
InterpolatedDerivatives[if1]
Out[7]=
{0,1,2,2,2,3,3,3,3,4,4,4,4,4,5,5,5}
You can prepare data suitable for Interpolation using just these
functions, and evaluation of the InterpolatingFunction at the knot
points by using MapThread:
In[8]:=
MapThread[{#1,Table[Derivative[j][if1][#1],{j,0,#2}]}&,{Knots[if1],
InterpolatedDerivatives[if1]}]
Out[8]=
{{0.,{1.}},{0.000894427,{1.0009,1.0009}},{
0.00178885,{1.00179,1.00179,1.00179}},{
0.0223764,{1.02263,1.02263,1.01221}},{0.042964,{1.0439,1.0439,1.03327}},{
0.0635516,{1.06562,1.06562,1.05476,1.04406}},{
0.122802,{1.13066,1.13066,1.13474,1.24789}},{
0.182053,{1.19968,1.19968,1.19839,1.13218}},{
0.241303,{1.27291,1.27291,1.27148,1.19976}},{
0.300553,{1.35061,1.35061,1.34911,1.2734,1.2429}},{
0.399063,{1.49043,1.49043,1.49071,1.508,1.81217}},{
0.497572,{1.64473,1.64473,1.64512,1.65398,1.64704}},{
0.596081,{1.815,1.815,1.81461,1.80037,1.56657}},{
0.69459,{2.00289,2.00289,2.00247,1.98685,1.72978}},{
0.793099,{2.21024,2.21024,2.20978,2.19265,1.90946,1.82406}},{
0.945494,{2.57409,2.57409,2.57351,2.56299,2.51458,2.68272}},{
1.,{2.71829,2.71829,2.71896,2.75402,3.44689,8.45155}}}
InterpolationOrder->1 is used here because thats what the value is in
NDSolve output (The reasons are complicated, but ultimately, it is to
keep approximations as local as possible)
In[9]:=
Interpolation[%, InterpolationOrder->1]
Out[9]=
InterpolatingFunction[{{0.,1.}},"<>"]
This is the same as the original object.
In[10]:=
% == if1
Out[10]=
True
Now you can write a function to merge two coniguous
InterpolatingFunction objects
In[11]:=
MergeInterpolatingFunction[if1_InterpolatingFunction,
if2_InterpolatingFunction] := Module[{k1, k2, d1, d2},
k1 = Knots[if1];
k2 = Knots[if2];
If[Last[k1] != First[k2], Return[$Failed]];
d1 = InterpolatedDerivatives[if1];
d2 = InterpolatedDerivatives[if2];
If[Last[d1] > First[d2],
k2 = Drop[k2,1];
d2 = Drop[d2,1],
(* else *)
k1 = Drop[k1,-1];
d2 = Drop[d1, -1]];
kd1 = MapThread[{#1,Table[Derivative[j][if1][#1],{j,0,#2}]}&,{k1,d1}];
kd2 = MapThread[{#1,Table[Derivative[j][if2][#1],{j,0,#2}]}&,{k2,d2}];
Interpolation[Join[kd1,kd2]]]
It can easily be generalized to several
In[12]:=
MergeInterpolatingFunction[ifs__InterpolatingFunction] :=
With[{ifl = {ifs}},
MergeInterpolatingFunction[MergeInterpolatingFunction[ifl[[1]],ifl[[2]]],
Sequence @@ Drop[ifl,2]]]
Combining the three solutions makes an interesting plot.
In[13]:=
MergeInterpolatingFunction[if1,if2,if3]
Out[13]=
InterpolatingFunction[{{0.,3.}},"<>"]
In[14]:=
Plot[%[t],{t,0,3}]