Google Groups no longer supports new Usenet posts or subscriptions. Historical content remains viewable.

Dismiss

Derivative of InterpolatingFunction

244 views

Skip to first unread message

Paul Abbott

unread,

Jul 8, 1996, 3:00:00 AM7/8/96

Andrei Constantinescu wrote:

> I have the following problem:
>
> fct = InterpolatingFunction[{{0., 3.14159}, {0., 3.14159}}, <>]
>
> ... so its an InterpolatingFunction of 2 variables .
>
> What I want now to do is the derivative of this function.
> But I fail !
>
> I get for example:
>
> In[69]:= Derivative[1,0][fct]
>
> (1,0)
> Out[69]= (InterpolatingFunction[{{25., 210.}, {0., 17.}}, <>])
> In[71]:= %69[100., 4.]

In The Mathematica Journal 4(2):31 the following appears:

Partial Derivatives

Presently, Mathematica cannot handle partial derivatives of
InterpolatingFunctions. The package DInterpolatingFunction.m, provided
by Hon Wah Tam (t...@wri.com) and included in the electronic
supplement, computes partial derivatives of two-dimensional
InterpolatingFunctions. Enhancements for higher dimensions will
eventually be incorporated into Mathematica.

After placing DInterpolatingFunction.m in your home directory, it can
be loaded using the command:

<< DInterpolatingFunction`

Define an array of data points:

g[x_, y_] := Sin[3 x] Cos[4 y]

tb = Table[{x, y, g[x, y]}, {x, 0, 1, .1}, {y, 0, 1, .05}];

and then interpolate it:

ifunc = Interpolation[Flatten[tb, 1]];

You can evaluate the interpolated function directly:

ifunc[0.23, 0.41]

-0.0440066

and plot its graph:

Plot3D[ifunc[x, y], {x, 0, 1}, {y, 0, 1}];

With the DInterpolatingFunction.m package, you can also compute partial
derivatives. Here is the derivative with respect to the first variable:

dx = Derivative[1, 0][ifunc];

dx[0.23, 0.41]

-0.160227

which can be compared with the expected value:

Derivative[1, 0][g][0.23, 0.41]

-0.159991

Similarly:

dy = Derivative[0, 1][ifunc];

dy[0.23, 0.41]

-2.53594

Derivative[0, 1][g][0.23, 0.41]

-2.54005

Here is the package:

(* DInterpolatingFunction.m *)

Unprotect[ InterpolatingFunction ];

BeginPackage[ "DInterpolatingFunction`", "Global`" ]

Begin[ "`Private`" ]

DiffElem[ element_, zerolist_ ] :=
Module[ { coeffs, dvdfs, order, var,
polynomial, j, newelement },

dvdfs = element[[1]];
coeffs = element[[2]];
order = Length[ coeffs ];

polynomial = dvdfs[[ order + 1 ]];
For[ j = order, j >= 1, j--,
polynomial = ( var + coeffs[[j]] ) * polynomial + dvdfs[[j]]
];

polynomial = D[ Expand[ polynomial ], var ];
dvdfs = CoefficientList[ polynomial, var ];

If[ Length[ dvdfs ] < order,
dvdfs = Join[ dvdfs, zerolist[[ order - Length[ dvdfs ] ]] ]
];
newelement = { dvdfs, zerolist[[order-1]] };

Return[ newelement ];

];

Interpolate[ x_, dvdfs_, coeffs_, tn_, order_ ] :=
Module[ { answer, xmtn, j },

answer = dvdfs[[order+1]];
xmtn = x - tn;
For[ j = order, j >= 1, j--,
answer = ( xmtn + coeffs[[j]] ) * answer + dvdfs[[j]]
];

Return[ answer ];

]

DiffLine[ line_, zerolist_, lastdimpts_ ] :=
Module[ { i, newline, dvdfs, coeffs, order, tn, x, df, element },

newline = Table[ DiffElem[ line[[i]], zerolist ],
{ i, 2, Length[ line ] } ];

{ dvdfs, coeffs } = newline[[1]];
tn = lastdimpts[[2]];
x = lastdimpts[[1]];
order = Length[ coeffs ];

df = Interpolate[ x, dvdfs, coeffs, tn, order ];
element = { { df, 0 }, { 0 } };
PrependTo[ newline, element ];

Return[ newline ];

]

InterpolatingFunction /:
Derivative[0,1][
InterpolatingFunction[ range_, table_ ] ] :=
Module[ { gridpoints, tbl, i, j, order,
zerolist, lastdimpts, newtable,answer },

gridpoints = table[[1]];
lastdimpts = Last[ gridpoints ];
tbl = table[[2]];

order = Max[ Table[ Length[ tbl[[1,i,2]] ],
{ i, Length[ tbl[[1]] ] } ] ];
zerolist = Table[ Table[ 0, { j, i } ], { i, order } ];

newtable = Table[ DiffLine[ tbl[[i]], zerolist, lastdimpts ],
{ i, Length[ tbl ] } ];

answer = InterpolatingFunction[ range, { gridpoints, newtable }
];
Return[ answer ];

];

InterpolatingFunction /:
Derivative[1,0][InterpolatingFunction[ range_, table_ ] ] :=
Module[ { gridpoints, ifunc, data, i, j, x, y },

gridpoints = table[[1]];

ifunc = InterpolatingFunction[ range, table ];
data = Table[ { y = gridpoints[[2,j]], x = gridpoints[[1,i]],
ifunc[x,y] },
{ j, Length[ gridpoints[[2]] ] },
{ i, Length[ gridpoints[[1]] ] } ];

data = Flatten[ data, 1 ];
ifunc = Interpolation[ data ];

ifunc = Derivative[0,1][ifunc];

data = Table[ { x = gridpoints[[1,i]], y = gridpoints[[2,j]],
ifunc[y,x] },
{ i, Length[ gridpoints[[1]] ] },
{ j, Length[ gridpoints[[2]] ] } ];
data = Flatten[ data, 1 ];
ifunc = Interpolation[ data ];

Return[ ifunc ];

];

End[]

EndPackage[]

Cheers,
Paul

_________________________________________________________________
Paul Abbott
Department of Physics Phone: +61-9-380-2734
The University of Western Australia Fax: +61-9-380-1014
Nedlands WA 6907 pa...@physics.uwa.edu.au
AUSTRALIA http://www.pd.uwa.edu.au/Paul
_________________________________________________________________

Mark James

unread,

Jul 8, 1996, 3:00:00 AM7/8/96

Andrei Constantinescu wrote:
>
> I have the following problem:
>
> fct = InterpolatingFunction[{{0., 3.14159}, {0., 3.14159}}, <>]
>
> ... so its an InterpolatingFunction of 2 variables .
>
> What I want now to do is the derivative of this function.

For the time being you will have to use the ND function as part
of NumericalMath`NLimit` in the standard packages.

You can use ND every time you need a derivative, or for greater
speed you may like to pre-compute the partial derivatives as
interpolating functions themselves using ND.

0 new messages