DInterpolatingFunction.m
- To: mathgroup at smc.vnet.net
- Subject: [mg4440] DInterpolatingFunction.m
- From: Andrei Constantinescu <constant at athena.polytechnique.fr>
- Date: Mon, 22 Jul 1996 03:10:40 -0400
- Sender: owner-wri-mathgroup at wolfram.com
>From news at styx.uwa.edu.au Mon Jul 1 06:12 WET 1996 (1.38.193.4/16.2) id AA20836; Mon, 1 Jul 1996 06:12:09 +0100 >From: Paul Abbott <paul at physics.uwa.edu.au> >To: mathgroup at smc.vnet.net >Subject: Re: Derivative of InterpolatingFunction >Organization: University of Western Australia 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 (tam at 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 paul at physics.uwa.edu.au AUSTRALIA http://www.pd.uwa.edu.au/Paul _________________________________________________________________ >From: Paul Abbott <paul at physics.uwa.edu.au> >To: mathgroup at smc.vnet.net >Organization: University of Western Australia >Subject: Re: Derivative of InterpolatingFunction 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 (tam at 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 paul at physics.uwa.edu.au AUSTRALIA http://www.pd.uwa.edu.au/Paul _________________________________________________________________ >From mrj at cs.usyd.edu.au Mon Jul 1 06:52 WET 1996 (1.38.193.4/16.2) id AA20885; Mon, 1 Jul 1996 06:52:54 +0100 >From: Mark James <mrj at cs.usyd.edu.au> >To: mathgroup at smc.vnet.net >Organization: The University of Sydney >Subject: Re: Derivative of InterpolatingFunction 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. -- Mark James | EMAIL : mrj at cs.usyd.edu.au | Basser Department of Computer Science, F09 | PHONE : +61-2-351-3423 | The University of Sydney NSW 2006 AUSTRALIA | FAX : +61-2-351-3838 | ==== [MESSAGE SEPARATOR] ====