Re: Gradient of a List
- To: mathgroup at smc.vnet.net
- Subject: [mg82614] Re: Gradient of a List
- From: Scott Hemphill <hemphill at hemphills.net>
- Date: Fri, 26 Oct 2007 05:21:23 -0400 (EDT)
- References: <10340059.1193232789565.JavaMail.root@m35> <ffpps9$l54$1@smc.vnet.net>
- Reply-to: hemphill at alumni.caltech.edu
DrMajorBob <drmajorbob at bigfoot.com> writes:
> data = Table[{x + RandomReal[], Sin@x + 0.1 RandomReal[]}, {x, 0, Pi,
> 0.1}];
> f = Interpolation[data, InterpolationOrder -> 3];
> {min, max} = data[[Ordering[data][[{1, -1}]], 1]];
> Quiet@Plot[f'[x], {x, min, max}, PlotRange -> All]
>
> I use Quiet because Plot sometimes samples outside the data range and
> throws the InterpolatingFunction::dmval message.
>
> Notice, however, the result isn't even close to Cos[x], and it changes
> quite a bit if you change the InterpolationOrder.
Of course, these problems are because of the noise in both the x and y
data values. Since Interpolation insists on passing exactly through
the points given, the interpolating function has to wiggle around a
lot to fit all the noise. The OP may not have any noise in his
independent variables (x,y) and may have little or none in his
function values.
Still, yours is an interesting problem. One way of handling it would
be to interpolate via weighted averages. For example, you could
assign a Gaussian weight to all the function values based on how close
the x value is to the x coordinates of the data:
(* Gaussian centered at x0, with standard deviation sig *)
pdf[x_,x0_,sig_] := 1/(Sqrt[2Pi]sig) Exp[-(x-x0)^2/(2sig^2)];
(* Gaussian weighted average of data, using sig = 0.5 *)
(* try using other values for sig *)
g[x_] = Block[{x0,y0,w},
x0 = data[[All,1]]; (* x-coordinates of data *)
y0 = data[[All,2]]; (* y-coordinates of data *)
w = pdf[x,#,0.5]& /@ x0; (* weight the x-coordinates *)
w /= Plus @@ w; (* Normalize the weights *)
w . y0 (* Return interpolated function value *)
];
Now you have a continuous function g[x], and you can plot it as well
as g'[x]. (Of course it is inefficient, since it recalculates all the
weights every time you call it. You could enter "foo[x_]=g[x];" and
then "foo[x]" wouldn't have that problem.)
One nice feature of Gaussian weights happens if you have a set of
equally spaced data points. If they extend infinitely in both
positive and negative directions, or *equivalently* the function
values are zero beyond the region of interest, then you can omit the
normalization step. (One example is in image processing, where the
region beyond the boundaries of the image may be assumed to be black.)
Then the Gaussian weighting is equivalent to convolution(1) with a
Gaussian kernel. This convolution has some nice properties. For
example, it is infinitely differentiable, because the Gaussian is.
Also, you can express its derivatives in the same form as the
convolution itself, i.e. the convolution of a Gaussian with a set of
data points.
The OP might be able to use the two-dimensional version of the the
Gaussian weighted interpolation above, but Mathematica's built-in
polynomial interpolation might work perfectly well.
(1) when I speak of convolution with a "data point" (x0,y0), I really
mean convolving with the function y0*DiracDelta[x0]. The result
is a y0 times a Gaussian centered at x0. Convolution with a
collection of data points gives the sum of all the Gaussians.
> On Wed, 24 Oct 2007 03:34:28 -0500, olalla <operez009 at ikasle.ehu.es> wrote:
>
> > Hi everybody,
> >
> > Does anybody know how can I get the "gradient" of a list of points?
> >
> > My real problem is:
> >
> > I have a scalar field previously obtained numerically that for a
> > given point (xi,yi) takes a value f(xi,yi). What I want to do is an
> > estimation of the gradient of this scalar field BUT I haven't got any
> > analytical function that expresses my field so I can't use the Grad
> > function.
> >
> > How can I solve this using Mathematica?
> >
> > Thanks in advance
> >
> > Olalla, Bilbao UPV/EHU
> >
> >
>
>
>
> --
>
> DrMajorBob at bigfoot.com
>
--
Scott Hemphill hemphill at alumni.caltech.edu
"This isn't flying. This is falling, with style." -- Buzz Lightyear