12
$\begingroup$

I need to find the "Maxima" and "Minima" on a B-Spline or more correct the points where the 2nd components of the derivate equal zero.

For example: 

g = BSplineFunction[{{1, 2}, {2, 4}, {3, -1}, {4, 2}}] ; 
dg=g';
Solve[{dg[t][[2]] == 0, 0 <= t <= 1}, t]

The problem is that "Solve" wont work for this kind of application, and "Minimize" or similar functions stop at the first finding.

Any ideas?
CX

|improve this question
$\endgroup$
20
$\begingroup$

You can use MeshFunctions to do the trick:

g = BSplineFunction[{RandomReal[1, 20], RandomReal[1, 20]}\[Transpose]];

dg = g';

ParametricPlot[
    g[t], {t, 0, 1},
    MeshFunctions -> Function[{x, y, t}, dg[t].{0, 1}],
    Mesh -> {{0}},
    MeshStyle -> Directive[AbsolutePointSize[5], Red]
    ]

extrema meshed spline

Here the MeshFunctions specifies the value of dg[t].{0, 1}, i.e. the $y$ component of the tangential vector of $g(t)$ at $t$, is used to generate mesh levels. Then Mesh -> {{0}} specifies that we only draw meshes where dg[t].{0, 1} == 0, which is exactly the extrema of $g(t)$.

To tell maxima from minima, use the RegionFunction:

maximapart = ParametricPlot[
                 g[t], {t, 0, 1},
                 PlotStyle -> Lighter[Blue, .6],
                 RegionFunction -> Function[{x, y, t}, g''[t].{0, 1} < 0],
                 MeshFunctions -> Function[{x, y, t}, dg[t].{0, 1}],
                 Mesh -> {{0}},
                 MeshStyle -> Directive[AbsolutePointSize[5], Red]
                 ]

minimapart = ParametricPlot[
                 g[t], {t, 0, 1},
                 PlotStyle -> Lighter[Brown, .6],
                 RegionFunction -> Function[{x, y, t}, g''[t].{0, 1} > 0],
                 MeshFunctions -> Function[{x, y, t}, dg[t].{0, 1}],
                 Mesh -> {{0}},
                 MeshStyle -> Directive[AbsolutePointSize[5], Blue]
                 ]

Show[{maximapart, minimapart}]

styled maxima and minima

Extracting the points is straightforward:

maximaptSet = Cases[maximapart, GraphicsComplex[pt_, __] :> pt, ∞][[1]];
maximaIdx = Cases[maximapart, Point[pt_] :> pt, ∞][[1]];
maximaptSet[[maximaIdx]]

minimaptSet = Cases[minimapart, GraphicsComplex[pt_, __] :> pt, ∞][[1]];
minimaIdx = Cases[minimapart, Point[pt_] :> pt, ∞][[1]];
minimaptSet[[minimaIdx]]
|improve this answer
$\endgroup$
  • $\begingroup$ Very nice use of 1D MeshFunction $\endgroup$ – chris Apr 21 '14 at 19:25
  • $\begingroup$ @chris Thanks. I realized this trick only a couple of days ago! :D $\endgroup$ – Silvia Apr 21 '14 at 19:27
  • $\begingroup$ could it be used to my problem? mathematica.stackexchange.com/q/9928/1089 just a thought? $\endgroup$ – chris Apr 21 '14 at 19:29
  • $\begingroup$ @chris Ah That's an interesting question. I upvoted at that time but didn't come up with any useful thought. Please allow me try tomorrow. Now I'm going to have some sleep :) $\endgroup$ – Silvia Apr 21 '14 at 19:32
  • $\begingroup$ Something like that ought to work? dat = GaussianRandomField[nn = 32] // Chop; dat /= Max[dat]; dat *= nn/2; dat2 = Table[{i, j, dat[[i, j]]}, {i, nn}, {j, nn}]; bs = BSplineFunction[dat2]; dbs = Function[{u, v}, Sqrt[Derivative[1, 0][bs][u, v][[3]]^2 + Derivative[0, 1][bs][u, v][[3]]^2] // Evaluate]; ParametricPlot3D[bs[u, v], {u, 0, 1}, {v, 0, 1}, MeshFunctions -> Function[{x, y, z, u, v}, dbs[u, v]], MeshStyle -> Directive[AbsolutePointSize[Small], Red], Mesh -> {{0}}] $\endgroup$ – chris Apr 21 '14 at 20:13
8
$\begingroup$

A more interesting example with multiple extrema..

 g = BSplineFunction[{{1, 2}, {2, 4}, {3, -1}, {4, 2}, {5, 0}, {6, 1}}];
 gp = g';
 gpy[t_?NumericQ] := gp[t][[2]];

This is utilising Plot to generate the curve and look for zero crossings, which we then pass as starting points to FindRoot

 loc = Flatten[
     t /. # & /@ 
        FindRoot[gpy[t] , Evaluate[ {t, Sequence @@ #[[;; , 1]]}]] & /@ 
          Select[ Partition[ 
           Cases[Plot[ gpy[t], {t, 0, 1}], Line[pts_] :> List[pts], 
               Infinity][[1, 1]] , 2, 1] , #[[1, 2]] #[[2, 2]] <= 0 & ] ]
 Show[{ ParametricPlot[ g[t], {t, 0, 1}], 
      Graphics@{PointSize[.02], Point[g[#] & /@ loc]}} ]

enter image description here

now distinguish min/max by looking at the second derivative:

 gpp = g'';
 min = Select[ loc, gpp[#][[2]] > 0 &]
 max = Select[ loc, gpp[#][[2]] < 0 &]
 Show[{ ParametricPlot[ g[t], {t, 0, 1}], 
     Graphics@{Red, PointSize[.02], Point[g[#] & /@ min], Blue, 
     PointSize[.02], Point[g[#] & /@ max]}} ]

enter image description here

|improve this answer
$\endgroup$
  • 1
    $\begingroup$ I'd never have hit on that. The problem is solved, thank you. But it's a quiet unusual way to solve such problems, isn't it? Doesn't Mathematica offer a simpler solutions? $\endgroup$ – cxkoda Apr 21 '14 at 16:59
  • $\begingroup$ I agree its odd there isn't a built in way to do this. For this case by the way you could as well do Table[gpy[t], {t, 0, 1,.01}] in place of the Cases[Plot..][[1,1]] because the function is nice enough we dont really need the adaptive sampling that comes from using Plot $\endgroup$ – george2079 Apr 21 '14 at 17:45
1
$\begingroup$

How about

 g = BSplineFunction[{{1, 2}, {2, 4}, {3, -1}, {4, 2}}];
 gN[t_?NumericQ] := g[t][[2]]

So that gN[0.1] returns a number.

Then 

 NMinimize[{gN[t], t > 0, t < 1}, t] 
 NMaximize[{gN[t], t > 0, t < 1}, t]

(* {1.01494,{t->0.75726}} {2.48728,{t->0.176073}} *)

works.

|improve this answer
$\endgroup$
  • $\begingroup$ Your method finds the point with the smallest gradient, but i need to find lokal extrema like in NMinimize[{Abs@dgN[t], t > 0, t < 1}, t] but there exist two of them. And thats my problem. Minimize stops at the first finding $\endgroup$ – cxkoda Apr 21 '14 at 13:50
  • $\begingroup$ Just using g instead of g' should do what you want (at least in this case with only one local min and max ) $\endgroup$ – george2079 Apr 21 '14 at 13:59
  • $\begingroup$ the example has two, the first at t~0.17 and the second at t~0.75. I need this bspline method to evaluate some data, therefore it's important for me to know every extrema. $\endgroup$ – cxkoda Apr 21 '14 at 14:06

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.