Re: Re: DiffMaps package

*To*: mathgroup at smc.vnet.net*Subject*: [mg91420] Re: Re: DiffMaps package*From*: "Jean-Marc Gulliet" <jeanmarc.gulliet at gmail.com>*Date*: Thu, 21 Aug 2008 04:14:47 -0400 (EDT)*References*: <20080820143021.4EB9F70@mach.vub.ac.be>

On Wed, Aug 20, 2008 at 4:30 PM, eabad wrote: > a couple of months ago, you were kind enough to reply to my email on how to compute extrema of 2D functions (see below). Imagine now that I have a list of data for the values of a two variable function, from which I can get an Interpolating function f[x,y] in mathematica following the standard procedure. When trying to apply your code below to f[x,y], the Solve or Nsolve command fails to provide the solutions of the equation D[f[x,y],{{x,y}},1]==0 for the zeros of the gradient in a specified x-y range. I wonder if you know a getaround for this problem. I am also somewhat worried about the discontinuity of the derivatives of f[x,y], I wonder how to obtain a smoother interpolation than the one provided by the default procedure. I'd be very grateful if you could briefly answer to these questions. [... Message crossposted to MathGroup ...] (Do not forget to post your query to the list/newsgroup MathGroup (comp.soft-sys.math.newsgroup) to get a broader on how to solve your problem.) With an interpolating function, you can use *FindRoot* to find zeros. Pay particular attention to the various ways of entering the initial value(s) for each variable for they influence the choice of a search method and allow to restrict the search to a specific interval. You can change the interpolating order with the option *InterpolatingOrder->n* and also tell Mathematica whether the fuctiion is periodic with the option *PeriodicInterpolation -> True*. The following example illustrate some technics you could use and some pitfall you may encounter. f[x_, y_] := Sin[x ] Cos[y] Plot3D[f[x, y], {x, 0, 2 Pi}, {y, 0, 2 Pi}] Plot3D[f[x, y] + 1/10 RandomReal[], {x, 0, 2 Pi}, {y, 0, 2 Pi}] data = Flatten[ Table[{{x, y}, f[x, y] + 1/10 RandomReal[]}, {x, 0, 2 Pi}, {y, 0, 2 Pi}], 1] g = Interpolation[data, InterpolationOrder -> 4] Plot3D[g[x, y], {x, 0, 2 Pi}, {y, 0, 2 Pi}] ListPlot3D[Flatten /@ data] grad = D[g[x, y], {{x, y}, 1}] FindRoot[grad, {{x, 1, 0, 2 Pi}, {y, 1, 0, 2 Pi}}] grad /. % // Chop Table[FindRoot[grad, {{x, x0, 0, 2 Pi}, {y, y0, 0, 2 Pi}}], {x0, 0, 2 Pi}, {y0, 0, 2 Pi}] grad /. % // Chop Regards, - Jean-Marc >>Modeler wrote: >> >>> Hi everybody, >>> >>> I need to determine the critical points of a 3D surface parameterized by a function f[x,y]. I was told there is a package called diffMaps (apparently not installed by default) which allows one to easily compute local minima and maxima as well as saddle points. Does anyone know how to get it or else can give any ideas about computation of critical points? Thanks a lot. >> >> >>Say that the function we want to study is the following: >> >> f[x_, y_] := 3*x*y - x^3 - y^3 >> >>We compute the first derivatives of f w.r.t. to x and y: >> >> D[f[x, y], {{x, y}, 1}] >> >> {-3 x^2 + 3 y, 3 x - 3 y^2} >> >>We equate each derivative to zero and solve for x and y: >> >> Solve[% == 0] >> >> {{x -> 0, y -> 0}, {x -> 1, y -> 1}, {x -> -(-1)^(1/3), >> y -> (-1)^(2/3)}, {x -> (-1)^(2/3), y -> -(-1)^(1/3)}} >> >>We want only the real points, so we discard the complex solutions: >> >> pts = >> Cases[%, {x -> vx_, y -> vy_} /; Im[vx] == 0 && Im[vy] == 0] >> >> {{x -> 0, y -> 0}, {x -> 1, y -> 1}} >> >>So we have two critical points: (0, 0) and (1, 1). >> >>Now, we compute the Hessian matrix, >> >> D[f[x, y], {{x, y}, 2}] >> >> {{-6 x, 3}, {3, -6 y}} >> >>and take its determinant: >> >> % // Det >> >> -9 + 36 x y >> >>Then we compute the values of the Hessian determinant at the critical >>points: >> >> % /. pts >> >> {-9, 27} >> >>Since the Hessian is negative at the point (0, 0), we conclude that this >>is a saddle point. >> >>Since the Hessian is positive at the point (1, 1), we conclude that this >>is a local maximum. >> >>We can "check" (or anticipate) these results by plotting the function: >> >> Plot3D[f[x, y], {x, -2, 2}, {y, -2, 2}, PlotRange -> All] >> >> ContourPlot[f[x, y], {x, -2, 2}, {y, -2, 2}, >> ColorFunction -> "SunsetColors"]