Skip to main content

splines - How to make the result of Piecewise be a closed interval?


Description


This question comes from two questions. Namely Q1 and Q2


The defintion of B-Spline basis function as shown below:


Let →U={u0,u1,…,um} a nondecreasing sequence of real numbers,i.e, ui≤ui+1i=0,1,2…m−1


$$N_{i,0}(u)= \begin{cases} 1 & u_i\leq u

  NBSpline[i_Integer, 0, knots_?(VectorQ[#, NumericQ] && OrderedQ[#] &),u_] /;
   i <= Length[knots] - 2 :=

   Piecewise[
    {{1, knots[[i + 1]] <= u < knots[[i + 2]]}},0]

  coeff[u_, i_, j_, knots_] /; knots[[i]] == knots[[j]] := 0;
  coeff[u_, i_, j_, knots_] := (u - knots[[i]])/(knots[[j]] - knots[[i]])

  NBSpline[i_Integer, p_Integer, knots_?(VectorQ[#, NumericQ] && OrderedQ[#] &), 
    u_] /;p > 0 && i + p <= Length[knots] - 2 :=
   Module[{init, res},
    init = Table[NBSpline[j, 0, knots, u], {j, i, i + p}];

    res = First@
    Nest[
     Dot @@@
      (Thread@
        {Partition[#, 2, 1],
         With[{m = p + 2 - Length@#},
          Table[
           {coeff[u, k + 1, k + m + 1, knots],
            coeff[u, k + m + 2, k + 2, knots]}, {k, i, i + Length@# - 2}]]}) &,
    init, p]

  ]

Compare to built-in function BSplineBasis


 sortResult[x_ /; x == 0] := 0;
 sortResult[res_] := MapAt[SortBy[#, Last] &, res, 1]

 pts = {{0, 0}, {0, 2}, {2, 3}, {4, 0}, {6, 3}, {8, 2}, {8, 0}};
 knots = {0, 0, 0, 1/5, 2/5, 3/5, 4/5, 1, 1, 1};

Comparsion 1



 NBSpline[#, 0, knots, u] & /@ Range[0, 8] // Simplify


enter image description here



 PiecewiseExpand@BSplineBasis[{0, knots}, #, u] & /@ Range[0, 8]


enter image description here




Comparsion 2


sortResult /@ (NBSpline[#, 1, knots, u] & /@ Range[0, 7] // Simplify)


enter image description here



sortResult /@ (PiecewiseExpand@BSplineBasis[{1, knots}, #, u] & /@ 
Range[0, 7] // Simplify)



enter image description here




By the comparsion, I found that the result of BSplineBasis is always a closed interval, whereas the result my function NBSpline is a half closed- half open interval.


My trail:


postProcess[res_ /; res == 0] := 0;

postProcess[res_] :=
 Module[{interval, pos, expr},
  interval =

   First@Simplify@res;
  pos =
   Position[
   SortBy[interval, Last], Less];
  expr =
   ReplacePart[
    SortBy[interval, Last], pos[[-1]] -> LessEqual];
  Piecewise[expr, 0]
]


Test


 postProcess /@ (NBSpline[#, 1, knots, u] & /@ Range[0, 7])


enter image description here



Question:


Is there any good/simple method(strategy) to deal with this problem? I think my method is complex.



enter image description here





Answer



a fix, specially treating the end interval:


 NBSpline[i_Integer, 0, knots_?(VectorQ[#, NumericQ] && OrderedQ[#] &),
      u_] /; i <= Length[knots] - 2 := (
   If[ knots[[i + 2]] != knots[[-1]],
       Piecewise[{{1, knots[[i + 1]] <= u < knots[[i + 2]]}}, 0],
       Piecewise[{{1, knots[[i + 1]] <= u <= knots[[i + 2]]}}, 0] ]);
 coeff[u_, i_, j_, knots_] /; knots[[i]] == knots[[j]] := 0;
 coeff[u_, i_, j_, knots_] := (u - knots[[i]])/(knots[[j]] - knots[[i]])

 NBSpline[i_Integer, p_Integer, 
      knots_?(VectorQ[#, NumericQ] && OrderedQ[#] &), u_] /; 
           p > 0 && i + p <= Length[knots] - 2 := 
 Module[{init, res}, 
      init = Table[NBSpline[j, 0, knots, u], {j, i, i + p}];
      res = First@
      Nest[Dot @@@ (Thread@{Partition[#, 2, 1], 
          With[{m = p + 2 - Length@#}, 
              Table[{coeff[u, k + 1, k + m + 1, knots], 
              coeff[u, k + m + 2, k + 2, knots]}, {k, i, 

                  i + Length@# - 2}]]}) &, init, p]]
  pts = {{1/10, 2/10}, {0, 2}, {2, 3}, {4, 0}, {6, 3}, {8, 2}, {8, 0}};
  knots = {0, 0, 0, 1/5, 2/5, 3/5, 4/5, 1, 1, 1};
  BSplinePlot1[pts : {{_, _} ..}, knots_, opts : OptionsPattern[Plot]] :=
  Module[{p = Length@First@Split[knots] - 1, a, b},
        {a, b} = {First[knots], Last[knots]};
    f[u_] = 
        Evaluate@
           Simplify@
               Total@MapIndexed[NBSpline[First@#2 - 1, p, knots, u] #1 &, 

                      pts] /. x_Piecewise :> Piecewise[x[[1]], Null];
    ParametricPlot[f[u], {u, a, b}, opts, PlotPoints -> 1000]]
    BSplinePlot1[pts, knots, ImageSize -> 600]

this bit: /. x_Piecewise :> Piecewise[x[[1]], Null] fixes the catchall behavior of Piecewise so you will get an error if you try to evaluate out of bounds. (not necessary for this example )


here is f[u] , note that we have closed just the end interval..


enter image description here


note as a practical matter @belisarius approach works as well, so long as you apply at the last stage:


 f[u_] = Evaluate@
          Simplify@

             Total@MapIndexed[NBSpline[First@#2 - 1, p, knots, u] #1 &, 
                pts] /. Less :> LessEqual /. x_Piecewise :> Piecewise[x[[1]], Null]

enter image description here


since the function is continuous it doesn't really matter if you close all the internal intervals.


Comments

Post a Comment

Popular posts from this blog

functions - Get leading series expansion term?

Given a function f[x] , I would like to have a function leadingSeries that returns just the leading term in the series around x=0 . For example: leadingSeries[(1/x + 2)/(4 + 1/x^2 + x)] x and leadingSeries[(1/x + 2 + (1 - 1/x^3)/4)/(4 + x)] -(1/(16 x^3)) Is there such a function in Mathematica? Or maybe one can implement it efficiently? EDIT I finally went with the following implementation, based on Carl Woll 's answer: lds[ex_,x_]:=( (ex/.x->(x+O[x]^2))/.SeriesData[U_,Z_,L_List,Mi_,Ma_,De_]:>SeriesData[U,Z,{L[[1]]},Mi,Mi+1,De]//Quiet//Normal) The advantage is, that this one also properly works with functions whose leading term is a constant: lds[Exp[x],x] 1 Answer Update 1 Updated to eliminate SeriesData and to not return additional terms Perhaps you could use: leadingSeries[expr_, x_] := Normal[expr /. x->(x+O[x]^2) /. a_List :> Take[a, 1]] Then for your examples: leadingSeries[(1/x + 2)/(4 + 1/x^2 + x), x] leadingSeries[Exp[x], x] leadingSeries[(1/x + 2 + (1 - 1/x...
Post a Comment
Read more

mathematical optimization - Minimizing using indices, error: Part::pkspec1: The expression cannot be used as a part specification

I want to use Minimize where the variables to minimize are indices pointing into an array. Here a MWE that hopefully shows what my problem is. vars = u@# & /@ Range[3]; cons = Flatten@ { Table[(u[j] != #) & /@ vars[[j + 1 ;; -1]], {j, 1, 3 - 1}], 1 vec1 = {1, 2, 3}; vec2 = {1, 2, 3}; Minimize[{Total@((vec1[[#]] - vec2[[u[#]]])^2 & /@ Range[1, 3]), cons}, vars, Integers] The error I get: Part::pkspec1: The expression u[1] cannot be used as a part specification. >> Answer Ok, it seems that one can get around Mathematica trying to evaluate vec2[[u[1]]] too early by using the function Indexed[vec2,u[1]] . The working MWE would then look like the following: vars = u@# & /@ Range[3]; cons = Flatten@{ Table[(u[j] != #) & /@ vars[[j + 1 ;; -1]], {j, 1, 3 - 1}], 1 vec1 = {1, 2, 3}; vec2 = {1, 2, 3}; NMinimize[ {Total@((vec1[[#]] - Indexed[vec2, u[#]])^2 & /@ R...
Post a Comment
Read more

What is and isn't a valid variable specification for Manipulate?

I have an expression whose terms have arguments (representing subscripts), like this: myExpr = A[0] + V[1,T] I would like to put it inside a Manipulate to see its value as I move around the parameters. (The goal is eventually to plot it wrt one of the variables inside.) However, Mathematica complains when I set V[1,T] as a manipulated variable: Manipulate[Evaluate[myExpr], {A[0], 0, 1}, {V[1, T], 0, 1}] (*Manipulate::vsform: Manipulate argument {V[1,T],0,1} does not have the correct form for a variable specification. >> *) As a workaround, if I get rid of the symbol T inside the argument, it works fine: Manipulate[ Evaluate[myExpr /. T -> 15], {A[0], 0, 1}, {V[1, 15], 0, 1}] Why this behavior? Can anyone point me to the documentation that says what counts as a valid variable? And is there a way to get Manpiulate to accept an expression with a symbolic argument as a variable? Investigations I've done so far: I tried using variableQ from this answer , but it says V[1...
Post a Comment
Read more