Re: Re: Gradient of a List

*To*: mathgroup at smc.vnet.net*Subject*: [mg82647] Re: [mg82614] Re: Gradient of a List*From*: DrMajorBob <drmajorbob at bigfoot.com>*Date*: Sat, 27 Oct 2007 06:00:27 -0400 (EDT)*References*: <10340059.1193232789565.JavaMail.root@m35> <ffpps9$l54$1@smc.vnet.net> <10446382.1193425327958.JavaMail.root@m35>*Reply-to*: drmajorbob at bigfoot.com

Scott, I hadn't noticed I was perturbing the x-values and using UNperturbed x to calculate y. No wonder the result looked so awful! Your Gaussian smoothing works amazingly well in this situation. I think I'd write it this way, however: ClearAll[pdf, g, x, data] pdf[xIdeal_, sig_] = Exp[-((# - xIdeal)/sig)^2/2]/(Sqrt[2 Pi] sig) &; g[data_, sigma_] := g[data, sigma] = Function[{x}, Evaluate@Block[{x0, y0}, {x0, y0} = Transpose@data; y0.Normalize[pdf[x, sigma] /@ x0, Total]]] data = Table[{x + RandomReal[], Sin@x + 0.1 RandomReal[]}, {x, 0, Pi, 0.1}]; Plot[g[data, 0.5]'[x], {x, min, max}, PlotRange -> All] Notice the "just-in-time" definition of g. Bobby On Fri, 26 Oct 2007 04:21:23 -0500, Scott Hemphill <hemphill at hemphills.net> wrote: > 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 >> > -- DrMajorBob at bigfoot.com