Services & Resources / Wolfram Forums
-----
 /
MathGroup Archive
1996
*January
*February
*March
*April
*May
*June
*July
*August
*September
*October
*November
*December
*Archive Index
*Ask about this page
*Print this page
*Give us feedback
*Sign up for the Wolfram Insider

MathGroup Archive 1996

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

Search the Archive

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] ====


  • Prev by Date: Replacement
  • Next by Date: help: Write data to file into many columns
  • Previous by thread: Compression Algorithm wanted
  • Next by thread: help: Write data to file into many columns