Quelle bernstein_polynomials.pvs Sprache: PVS

bernstein_polynomials: THEORY

%--------------------------------------------------------------------------
% Bernstein Polynomials
%
% Version 1      April 2010
%
% Anthony Narkawicz       NASA
%--------------------------------------------------------------------------

BEGIN

  IMPORTING polynomials,
            sigma_swap[nat],
            binomial_identities

  x,y,t,Ma   : VAR real
  i,j,m,p,bi,n : VAR nat
  k,r      : VAR posnat

  Bernstein_Polynomial: TYPE = [# index: nat, bern_seq: [upto(index) -> real] #]

  b,b1,b2 : VAR Bernstein_Polynomial
  a,a1,a2 : VAR sequence[real] %---- An ordinary polynomial ----%

  Bern(i:nat,n:above(i-1))(x) : real = C(n,i)*x^i*(1-x)^(n-i)

  Bern_top: LEMMA Bern(p,p)(x) = x^p

  Bern_bottom: LEMMA Bern(0,k)(x) = (1-x)^k

  Bern_middle_zero: LEMMA FORALL (i:nat,nn:above(i-1)): i/=0 IMPLIES Bern(i,nn)(0) = 0

  Bern_middle_one: LEMMA FORALL (i:nat,nn:above(i-1)): i/=nn IMPLIES Bern(i,nn)(1) = 0

  Bern_Polynomial: LEMMA
        FORALL (i:nat, nn:above(i-1)):
        Bern(i,nn) =
        polynomial((LAMBDA (j:nat): IF j<i OR j>nn THEN 0 ELSE (-1)^(j-i)*C(nn,j)*C(j,i) ENDIF),nn)

  Bernstein_Recursion: LEMMA k>r IMPLIES Bern(r,k)(x) = (1-x)*Bern(r,k-1)(x) + x*Bern(r-1,k-1)(x)

  Bern_nonnegative: LEMMA
        0<=x AND x<=1 IMPLIES FORALL (i:nat, nn:above(i-1)): Bern(i,nn)(x) >= 0

  Bern_positive: LEMMA
       0<x AND x<1 IMPLIES FORALL (i:nat, nn:above(i-1)): Bern(i,nn)(x) > 0

  % Degree raising for Bernstein polynomials

  Bernstein_degree_raise: LEMMA
    k>p IMPLIES Bern(p,k-1)(x) = ((k-p)/k)*Bern(p,k)(x) + ((p+1)/k)*Bern(p+1,k)(x)

  degree_raise(b): Bernstein_Polynomial =
    LET
      m = b`index,
      bseq = b`bern_seq
    IN
      (# index := m+1,
      bern_seq := (LAMBDA (i:upto(m+1)):
        IF i = 0 THEN bseq(0)
ELSIF i = m+1 THEN bseq(m)
ELSE
   (bseq(i)*(m+1-i) + bseq(i-1)*i)/(m+1)
ENDIF) #)

  degree_raise_power(b,j): RECURSIVE Bernstein_Polynomial =
    IF j=0 THEN b
    ELSE degree_raise(degree_raise_power(b,j-1))
    ENDIF
    MEASURE j

  Bernstein_index_degree_raise: LEMMA k>=p IMPLIES x*Bern(p,k)(x) = ((p+1)/(k+1))*Bern(p+1,k+1)(x)

  % The Bernstein polynomials form a partition of unity

  Bernstein_partition_of_unity: LEMMA
    sigma(0,p,(LAMBDA (i:nat): IF i>p THEN 0 ELSE Bern(i,p)(x) ENDIF)) = 1

  Bern_le_1: LEMMA 0<=x AND x<=1 IMPLIES FORALL (i:nat, nn:above(i-1)): Bern(i,nn)(x) <= 1

  Bernstein_conversion: LEMMA
    k>=p IMPLIES
    x^p = sigma(p,k,(LAMBDA (i:nat): IF i<p OR i>k THEN 0 ELSE (C(i,p)/C(k,p))*Bern(i,k)(x) ENDIF))

  Bernstein_conversion_2: LEMMA
    k>=p IMPLIES
    x^p = sigma(0,k,(LAMBDA (i:nat): IF i<p OR i>k THEN 0 ELSE (C(i,p)/C(k,p))*Bern(i,k)(x) ENDIF))

  % Bernstein polynomials - linear combinations of the basic Bernstein polynomials B_ik

  Bern_poly(b): [real -> real] = LAMBDA (x:real):
    sigma(0,b`index,(LAMBDA (i:nat): IF i > b`index THEN 0 ELSE b`bern_seq(i)*Bern(i,b`index)(x) ENDIF))

  degree_raise_id: LEMMA Bern_poly(b) = Bern_poly(degree_raise(b))

  degree_raise_power_id: LEMMA Bern_poly(b) = Bern_poly(degree_raise_power(b,j))

  Bern_poly_reverse(b): Bernstein_Polynomial =
    (# index:=b`index, bern_seq := (LAMBDA (i:upto(b`index)): b`bern_seq(b`index-i)) #)

  Bern_poly_reverse_sym_half: LEMMA Bern_poly(b)(1/2 + x) = Bern_poly(Bern_poly_reverse(b))(1/2-x)

  Bern_poly_reverse_sym: LEMMA Bern_poly(b)(x) = Bern_poly(Bern_poly_reverse(b))(1-x)

  Bern_poly_inverse(b): sequence[real] = LAMBDA (j:nat):
        IF j>b`index THEN
   0
ELSE
   sigma(0,j,(LAMBDA (i:nat): IF i > j THEN 0 ELSE b`bern_seq(i)*(-1)^(j-i)*(C(b`index,j)*C(j,i)) ENDIF))
ENDIF

  % Structural Property Decomposition For Bernstein Polynomials

  P: VAR [Bernstein_Polynomial -> bool]

  Bern_structural_decomp: LEMMA
    ((FORALL (b1,b2: Bernstein_Polynomial): b1`index = b2`index AND P(b1) AND P(b2)
      IMPLIES P((# index:= b1`index,bern_seq := (LAMBDA (i:upto(b1`index)): b1`bern_seq(i) + b2`bern_seq(i)) #)))
    AND
    (FORALL (i,m:nat,x:real): i<=m IMPLIES
     P((# index := m, bern_seq := (LAMBDA (j:upto(m)): IF j = i THEN x ELSE 0 ENDIF) #))))
    IMPLIES
    (FORALL (b:Bernstein_Polynomial): P(b))

  % The function that takes a sequence corresponing to a polynomial in the monomial
  % basis and returns the sequence corresponing to the polynomial in the Bernstein basis. The number p
  % is the "degree" of the polynomial in the monomial basis, and m is the Bernstein degree.

  Bernstein_polynomial(a,p,m): Bernstein_Polynomial =
   (# index := m,
   bern_seq := LAMBDA (s: upto(m)): sigma(0,s,LAMBDA(i:nat):IF i>s OR i>p THEN 0 ELSE a(i)*(C(s,i)/C(m,i)) ENDIF) #)

  Bernstein_equivalence: LEMMA
    m>=p IMPLIES polynomial(a,p) = Bern_poly(Bernstein_polynomial(a,p,m))

  Bern_poly_inverse_def: LEMMA
    polynomial(Bern_poly_inverse(b),b`index) = Bern_poly(b)

  % Splitting the Interval

  Bern_subdiv_left(b): Bernstein_Polynomial =
    (# index := b`index,
       bern_seq := LAMBDA (i:upto(b`index)): (1/2^i)*sigma(0,i,
                   LAMBDA (j:nat): IF j > i THEN 0 ELSE C(i,j)*b`bern_seq(j) ENDIF) #)

  Bern_subdiv_right(b): Bernstein_Polynomial =
    (# index := b`index,
    bern_seq := LAMBDA (i:upto(b`index)): (1/2^(b`index-i))*sigma(0,b`index-i,
                LAMBDA (j:nat): IF j > b`index-i THEN 0 ELSE C(b`index-i,j)*b`bern_seq(b`index-j) ENDIF) #)

  Bern_subdiv_left_id: LEMMA Bern_poly(Bern_subdiv_left(b)) = Bern_poly(b) o (LAMBDA (x:real): (1/2)*x)

  Bern_subdiv_left_reverse: LEMMA Bern_poly(Bern_poly_reverse(Bern_subdiv_left(Bern_poly_reverse(b)))) =
    Bern_poly(b) o (LAMBDA (x:real): (1/2)*(x+1))

  Bern_subdiv_right_id: LEMMA Bern_poly(Bern_subdiv_right(b)) = Bern_poly(b) o (LAMBDA (x:real): (1/2)*(x+1))

  Bern_subdiv_ge: LEMMA
      (FORALL (y:real): 0<=y AND y<=1 IMPLIES Bern_poly(Bern_subdiv_right(b))(y) >= Ma)
    AND
    (FORALL (y:real): 0<=y AND y<=1 IMPLIES Bern_poly(Bern_subdiv_left(b))(y) >= Ma)
    IMPLIES
    (0<=x AND x<=1 IMPLIES Bern_poly(b)(x) >= Ma)

  %-------------------------------------------------------%
  %-------------------------------------------------------%
  % Splitting the interval another way
  %-------------------------------------------------------%
  %-------------------------------------------------------%

  Bern_sweep(b)(p:upto(b`index)): RECURSIVE {bp: Bernstein_Polynomial | bp`index = b`index}
    = (# index := b`index,
      bern_seq :=
      IF p = 0 THEN b`bern_seq
      ELSE
        (LAMBDA (i:upto(b`index)):
IF i<p THEN Bern_sweep(b)(p-1)`bern_seq(i)
ELSE (1/2)*(Bern_sweep(b)(p-1)`bern_seq(i-1) + Bern_sweep(b)(p-1)`bern_seq(i))
ENDIF)
      ENDIF #) MEASURE p

  Bern_sweep_expand: LEMMA
    p<=b`index IMPLIES
    Bern_sweep(b)(p) =
      (# index := b`index,
      bern_seq :=
      (LAMBDA (i:upto(b`index)):
IF i<=p THEN
   (1/(2^i))*sigma(0,i,(LAMBDA (j:nat): IF j > i THEN 0 ELSE b`bern_seq(j)*C(i,j) ENDIF))
ELSE
   (1/(2^p))*sigma(0,p,(LAMBDA (j:nat): IF j > p THEN 0 ELSE b`bern_seq(j+(i-p))*C(p,j) ENDIF))
ENDIF) #)

  Bern_subdiv_l(b): Bernstein_Polynomial = Bern_sweep(b)(b`index)

  Bern_subdiv_r(b): Bernstein_Polynomial =
        (# index := b`index,
      bern_seq := (LAMBDA (i:upto(b`index)): Bern_sweep(b)(b`index - i)`bern_seq(b`index)) #)

  Bern_subdiv_l_id: LEMMA Bern_poly(Bern_subdiv_l(b)) = Bern_poly(b) o (LAMBDA (x:real): (1/2)*x)

  Bern_subdiv_r_id: LEMMA Bern_poly(Bern_subdiv_r(b)) = Bern_poly(b) o (LAMBDA (x:real): (1/2)*(x+1))

  %-------------------------------------------------------%
  %-------------------------------------------------------%
  % End of splitting the interval another way
  %-------------------------------------------------------%
  %-------------------------------------------------------%

  % Bernstein Coefficients

  Bern_coeff(a,p,m)(s: upto(m)): real =
    LET msp = min(s,p) IN sigma(0,msp,LAMBDA(i:nat):IF i>msp THEN 0 ELSE a(i)*(C(s,i)/C(m,i)) ENDIF)

  Bernstein_le: LEMMA
    (FORALL (s:upto(m)): Bern_coeff(a,p,m)(s) <= Ma) AND 0<=x AND x<=1 AND m>=p IMPLIES polynomial(a,p)(x) <= Ma

  Bernstein_lt: LEMMA
    (FORALL (s:upto(m)): Bern_coeff(a,p,m)(s) < Ma) AND 0<=x AND x<=1 AND m>=p IMPLIES polynomial(a,p)(x) < Ma

  Bernstein_ge: LEMMA
    (FORALL (s:upto(m)): Bern_coeff(a,p,m)(s) >= Ma) AND 0<=x AND x<=1 AND m>=p IMPLIES polynomial(a,p)(x) >= Ma

  Bernstein_gt: LEMMA
    (FORALL (s:upto(m)): Bern_coeff(a,p,m)(s) > Ma) AND 0<=x AND x<=1 AND m>=p IMPLIES polynomial(a,p)(x) > Ma

  % For Bernstein Polynomials

  Bern_poly_le: LEMMA
    (FORALL (s:upto(b`index)): b`bern_seq(s) <= Ma) AND 0<=x AND x<=1 IMPLIES Bern_poly(b)(x) <= Ma

  Bern_poly_lt: LEMMA
    (FORALL (s:upto(b`index)): b`bern_seq(s) < Ma) AND 0<=x AND x<=1 IMPLIES Bern_poly(b)(x) < Ma

  Bern_poly_ge: LEMMA
    (FORALL (s:upto(b`index)): b`bern_seq(s) >= Ma) AND 0<=x AND x<=1 IMPLIES Bern_poly(b)(x) >= Ma

  Bern_poly_gt: LEMMA
    (FORALL (s:upto(b`index)): b`bern_seq(s) > Ma) AND 0<=x AND x<=1 IMPLIES Bern_poly(b)(x) > Ma

  % Quotients

  Bern_poly_quotient_le: LEMMA
    b1`index = b2`index AND
    (FORALL (s:upto(b2`index)): b2`bern_seq(s) > 0) AND
    (FORALL (s:upto(b1`index)): b1`bern_seq(s)/b2`bern_seq(s) <= Ma)
    AND 0<=x AND x<=1 IMPLIES
    (Bern_poly(b2)(x) > 0 AND
    Bern_poly(b1)(x)/Bern_poly(b2)(x) <= Ma)

  % Splitting

  Bern_poly_split_le: LEMMA
     LET bsl = Bern_subdiv_l(b),bsr = Bern_subdiv_r(b) IN
        (FORALL (x:real): 0<=x AND x<=1 IMPLIES Bern_poly(bsr)(x)<=Ma) AND
        (FORALL (x:real): 0<=x AND x<=1 IMPLIES Bern_poly(bsl)(x)<=Ma)
        IMPLIES (FORALL (x:real): 0<=x AND x<=1 IMPLIES Bern_poly(b)(x)<=Ma)

  Bern_poly_split_lt: LEMMA
     LET bsl = Bern_subdiv_l(b),bsr = Bern_subdiv_r(b) IN
        (FORALL (x:real): 0<=x AND x<=1 IMPLIES Bern_poly(bsr)(x)<Ma) AND
        (FORALL (x:real): 0<=x AND x<=1 IMPLIES Bern_poly(bsl)(x)<Ma)
        IMPLIES (FORALL (x:real): 0<=x AND x<=1 IMPLIES Bern_poly(b)(x)<Ma)

  Bern_poly_split_ge: LEMMA
     LET bsl = Bern_subdiv_l(b),bsr = Bern_subdiv_r(b) IN
        (FORALL (x:real): 0<=x AND x<=1 IMPLIES Bern_poly(bsr)(x)>=Ma) AND
        (FORALL (x:real): 0<=x AND x<=1 IMPLIES Bern_poly(bsl)(x)>=Ma)
        IMPLIES (FORALL (x:real): 0<=x AND x<=1 IMPLIES Bern_poly(b)(x)>=Ma)

  Bern_poly_split_gt: LEMMA
     LET bsl = Bern_subdiv_l(b),bsr = Bern_subdiv_r(b) IN
        (FORALL (x:real): 0<=x AND x<=1 IMPLIES Bern_poly(bsr)(x)>Ma) AND
        (FORALL (x:real): 0<=x AND x<=1 IMPLIES Bern_poly(bsl)(x)>Ma)
        IMPLIES (FORALL (x:real): 0<=x AND x<=1 IMPLIES Bern_poly(b)(x)>Ma)

  % COLLECTIONS OF BERNSTEIN POLYNOMIALS

  % Next are recursive definitions of min and max for bernstein polynomials

  Bern_min_rec(b)(q:upto(b`index)): RECURSIVE
    {r:real | (EXISTS (iq:upto(q)): r = b`bern_seq(iq)) AND (FORALL (iq:upto(q)): r <= b`bern_seq(iq))} =
    LET y = b`bern_seq(q),
        x = IF q = 0 THEN y ELSE Bern_min_rec(b)(q-1) ENDIF
    IN min(y,x)
    MEASURE q

  Bern_max_rec(b)(q:upto(b`index)): RECURSIVE
   {r:real | (EXISTS (iq:upto(q)): r = b`bern_seq(iq)) AND (FORALL (iq:upto(q)): r >= b`bern_seq(iq))} =
   LET y = b`bern_seq(q),
       x = IF q = 0 THEN y ELSE Bern_max_rec(b)(q-1) ENDIF
   IN max(y,x)
   MEASURE q

  Bern_min(b) :
   {r : real | (EXISTS (ib:upto(b`index)): r = b`bern_seq(ib)) AND
               (FORALL (ib:upto(b`index)): r<=b`bern_seq(ib))} =
   Bern_min_rec(b)(b`index)

  Bern_max(b) :
   {r : real | (EXISTS (ib:upto(b`index)): r = b`bern_seq(ib)) AND
               (FORALL (ib:upto(b`index)): r>=b`bern_seq(ib))} =
    Bern_max_rec(b)(b`index)

  Bern_min_bound : LEMMA 0<=x AND x<=1 IMPLIES Bern_poly(b)(x)>=Bern_min(b)

  Bern_max_bound : LEMMA 0<=x AND x<=1 IMPLIES Bern_poly(b)(x)<=Bern_max(b)

  % A function to compute the current min and max estimates follows. One should note that if b_0,...,b_n are
  % Bernstein coefficients for a polynomial, then the polynomial attains the values of b_0 and b_n at t = 0
  % and t = 1.

  Bern_min_est(b): real = min(b`bern_seq(0),b`bern_seq(b`index))

  Bern_max_est(b): real = max(b`bern_seq(0),b`bern_seq(b`index))

  % The next construction, Bernstein_Collection is a type. An element of this type is a collection of
  % Bernstein Polynomials, and associated to each one is a pair consisting of its minimum and maximum entries

  Bernstein_MinMax: TYPE = [# bernpoly: Bernstein_Polynomial, bernmin:real, bernmax:real #]

  Bern_minmax(b) : Bernstein_MinMax = (# bernpoly := b, bernmin := Bern_min(b), bernmax := Bern_max(b) #)

  Bernstein_MinMax_Collection: TYPE = [# number_polys : posnat,
                                       bern_poly_number : [below(number_polys) -> Bernstein_MinMax] #]

  BC, BC1, BC2: VAR Bernstein_MinMax_Collection

  Bern_concat(BC1,BC2): Bernstein_MinMax_Collection =
    LET NN = BC1`number_polys + BC2`number_polys IN
    (# number_polys := NN,
       bern_poly_number := LAMBDA (i:below(NN)): IF i<=BC1`number_polys-1 THEN BC1`bern_poly_number(i)
                                                 ELSE BC2`bern_poly_number(i-BC1`number_polys) ENDIF #)

  Bern_collection_min_est_rec(BC)(np:below(BC`number_polys)) : RECURSIVE real =
    IF np = 0 THEN Bern_min_est(BC`bern_poly_number(0)`bernpoly)
    ELSE min(Bern_min_est(BC`bern_poly_number(0)`bernpoly),Bern_collection_min_est_rec(BC)(np-1)) ENDIF
  MEASURE np

  Bern_collection_min_est(BC): real =
    IF BC`number_polys /= 0 THEN Bern_collection_min_est_rec(BC)(BC`number_polys-1) ELSE 0 ENDIF

  Bern_collection_reduce_min_rec(BC)(np:below(BC`number_polys)): RECURSIVE Bernstein_MinMax_Collection =
    LET minest = Bern_collection_min_est(BC),
        thispoly:Bernstein_MinMax_Collection = (# number_polys := 1,
                                                  bern_poly_number := LAMBDA (i:below(1)): BC`bern_poly_number(np) #)
    IN
    IF np = 0 THEN thispoly
    ELSIF BC`bern_poly_number(np)`bernmin > minest THEN Bern_collection_reduce_min_rec(BC)(np-1)
    ELSE  Bern_concat(Bern_collection_reduce_min_rec(BC)(np-1),thispoly) ENDIF
  MEASURE np

  Bern_collection_reduce_min(BC) : Bernstein_MinMax_Collection =
    Bern_collection_reduce_min_rec(BC)(BC`number_polys-1)

  % The minimum coefficient

  Bern_min_coeff_rec(BC)(np:below(BC`number_polys)) : RECURSIVE real =
    LET current_poly_min = BC`bern_poly_number(np)`bernmin
    IN
    IF np = 0 THEN current_poly_min
    ELSE min(current_poly_min,Bern_min_coeff_rec(BC)(np-1)) ENDIF
  MEASURE np

  Bern_min_coeff(BC): real = Bern_min_coeff_rec(BC)(BC`number_polys-1)

  % Next, we subdivide each bernstein polynomial in the collection.

  Bern_coll_subdiv_total_left(BC) : Bernstein_MinMax_Collection =
    (# number_polys := BC`number_polys,
       bern_poly_number := LAMBDA (np: below(BC`number_polys)):
                           Bern_minmax(Bern_subdiv_left(BC`bern_poly_number(np)`bernpoly)) #)

  Bern_coll_subdiv_total_right(BC) : Bernstein_MinMax_Collection =
    (# number_polys := BC`number_polys,
       bern_poly_number := LAMBDA (np: below(BC`number_polys)):
                           Bern_minmax(Bern_subdiv_right(BC`bern_poly_number(np)`bernpoly)) #)

  Bern_coll_subdiv_total(BC) : Bernstein_MinMax_Collection =
   Bern_concat(Bern_coll_subdiv_total_left(BC),Bern_coll_subdiv_total_right(BC))

  Bern_coll_subdiv_reduced(BC) : Bernstein_MinMax_Collection =
    Bern_collection_reduce_min(Bern_coll_subdiv_total(BC))

  Bern_split(b)(n :nat) : recursive Bernstein_MinMax_Collection =
    IF n=0 THEN (# number_polys := 1, bern_poly_number := (LAMBDA (i:below(1)): Bern_minmax(b)) #)
    ELSE
    Bern_coll_subdiv_reduced(Bern_split(b)(n-1))
    ENDIF
    MEASURE n

  Bern_split_min_coeff(a,p,m)(n:nat) : real =
    LET b = Bernstein_polynomial(a,p,m)
    IN Bern_min_coeff(Bern_split(b)(n))

END bernstein_polynomials

quality60%

¤ Dauer der Verarbeitung: 0.14 Sekunden (vorverarbeitet) ¤

Wurzel

Suchen

Beweissystem der NASA

Beweissystem Isabelle

NIST Cobol Testsuite

Cephes Mathematical Library

Wiener Entwicklungsmethode

Haftungshinweis

Die Informationen auf dieser Webseite wurden nach bestem Wissen sorgfältig zusammengestellt. Es wird jedoch weder Vollständigkeit, noch Richtigkeit, noch Qualität der bereit gestellten Informationen zugesichert.

Bemerkung:

Die farbliche Syntaxdarstellung ist noch experimentell.