with(combinat); #==================================================================== print(Posted: Nov. 18, 2010 ): print(Last updated: Nov. 29, 2010): print( This is ProveGAS, a Maple package ): print(to prove Global Asymptotic Stability of certain rational difference equations); print(using the method of proving that some polynomial is positive.); print(): print(Please report bugs to eahogan at math dot rutgers dot edu): print(): print(The most current version of this package): print( is available from): print(http://www.math.rutgers.edu/~eahogan/maple/ProveGAS.txt .): print(For more information about the accompanying paper and Web Books); print(see the website:); print(http://www.math.rutgers.edu/~eahogan/GAS.html); print(For a list of the main procedures type Help();): #==================================================================== #IsPosNoPrint(H,Z) #Inputs: # H - a polynomial # Z - variables in H #Outputs: # true if H(Z)>0 when all variables are >=0 (except H(Z)>=0 when all vars=0) # false if there is an assignment, Z', of the variables in the positive # orthant s.t. H(Z')<=0 (or H(0)<0) # FAIL if the sufficient conditions checked are inconclusive #Method: # Checks two sufficient conditions: # 1. All coefficients are positive, and the constant term is positive or 0. # If the constant term is 0 then check that there is no subset of the # variables s.t. setting that subset 0 yields the 0 polynomial. # 2. If the only negative coefficient is on a term from the set # {Z[i]*Z[j], Z[i], Z[i]^2} then checks that the sub-polynomial # consisting of the terms from this set (for all i,j) is still positive # using the quadratic formula. #Also: # "NoPrint" means that nothing is printed in this implementation of IsPos IsPosNoPrint := proc(H,Z) local Os,C1,neg1,const,n,sets,flag,II,h,LargestNegDeg,tDeg,look,ii,loc, checkPos,flag2,a,b,c,disc,Cd,j,sets2,flag3,JJ,disc0,CheckPos0,d; n := nops(Z); sets := powerset(n); Os := [op(sort(collect(H, Z, 'distributed'),order=tdeg(op(Z))))]; #list of terms in H C1 := [seq([coeffs(Os[i],Z),Os[i]/coeffs(Os[i],Z),i],i=1..nops(Os))]; #list of [coeff,vars,term#] neg1 := select(i->evalb(evalf(i[1])<=0),C1); #list of terms w/ negative coeff const := coeftayl(H, Z=[0$n], [0$n]); #constant term in H if nops(neg1)=0 and const>0 then #all coeffs positive and constant term positive RETURN(true); elif nops(neg1)=0 and const=0 then #all coeffs positive and constant term = 0 for II in sets minus {{seq(i,i=1..n)}} do h[II] := subs({seq(Z[i]=0, i in II)},H); if h[II]=0 then RETURN(false); fi; od; RETURN(true); elif evalf(const)<0 then #constant term is negative RETURN(false); else #there are negative coefficients and constant term is >= 0 if max({seq(degree(neg1[1][2],Z[i]),i=1..n)})>2 then RETURN(FAIL); else LargestNegDeg := add(degree(neg1[1][2],Z[i]),i=1..n); look := neg1[1][3]; #the index in Os of the first term w/ neg. coeff loc := look; #if look=1 then this is needed because we won't enter into the for loop below tDeg := LargestNegDeg; #will provide ending condition for the for loop for ii from 1 to look-1 while tDeg = LargestNegDeg do tDeg := add(degree(C1[look-ii][2],Z[j]), j=1..n); loc := C1[look-ii][3]; #will be the position in Os from which to start checkPos od; if LargestNegDeg=tDeg then #we got all the way to look-1 with the same total degree checkPos := add(Os[j], j=loc..nops(Os)); else #we stopped before ii=look-1 so we must add 1 to loc checkPos := add(Os[j], j=loc+1..nops(Os)); fi; flag2 := true; #will serve as "flag" in the pseudocode algorithm for j from 1 to n while flag2 do if degree(checkPos,Z[j])=2 then a := coeff(checkPos,Z[j],2); b := coeff(checkPos,Z[j],1); c := coeff(checkPos,Z[j],0); disc := b^2-4*a*c; Cd := [coeffs(collect(disc,Z,'distributed'),Z)]; #list of coeffs in disc if {seq(sign(evalf(Cd[i])), i=1..nops(Cd))}={-1} then #all coeffs are negative sets2 := powerset({seq(i, i=1..j-1), seq(i,i=j+1..n)}) minus {{seq(i, i=1..j-1), seq(i,i=j+1..n)}}; flag3 := true; for JJ in sets2 while flag3 do d[JJ] := subs({seq(Z[i]=0, i in JJ)},disc); if d[JJ]=0 then flag3 := false; fi; od; if flag3 then #there are no subsets of {1,...,j-1,j+1,...,n} s.t. setting those vars to 0 yields disc=0 disc0 := subs({seq(Z[i]=0, i=1..j-1),seq(Z[i]=0,i=j+1..n)},disc); CheckPos0 := subs({seq(Z[i]=0, i=1..j-1),seq(Z[i]=0,i=j+1..n)},checkPos); if disc0 < 0 or (disc0 = 0 and coeff(CheckPos0,Z[j])*coeff(CheckPos0,Z[j],2)>=0) then #keep flag=true else flag2 := false; fi; else flag2 := false; fi; else #there are positive coefficients in disc flag2 := false; fi; else #degree(checkPos,Z[j]) <> 2 flag2 := false; fi; od; if flag2 then RETURN(true) fi; fi; RETURN(FAIL); fi; end: #==================================================================== #SmallBoxMethodNoPrint(h,Z,B,N) #Inputs: # h - a polynomial # Z - variables in h # B - list of positive real numbers of length nops(Z) # N - an integer #Outputs: # true if h(Z)>0 for 0<=Z[i]<=B[i] for all i, except possibly # h(0,...,0)>=0 and h(B)>=0 # false if there is a n assignment, Z', of variables s.t. h(Z')<0 # FAIL if the method is inconclusive #Method: # Splits the hyperbox 0<=Z[i]<=B[i] for all i into 2^N smaller hyperboxes # (split each axis into N intervals) and uses (the correct version of) # IsPos on each of the smaller hyperboxes. #Also: # "NoPrint" means that almost nothing is printed in this implementation of # SmallBoxMethod SmallBoxMethodNoPrint := proc(h,Z,B,NN) local n,interval,i,j,L,Li,Y,H,Hstar,Hstar0,Os,C1,neg1,const,tryFlips,k,flag,flag2,smaller, disc,checkPos,tDegNeg,look,tDeg,ii,loc,a,b,c,Cd,solNumer,d,Bad,Pos,N,nt,h2,B2; N := 2; if nargs=5 then nt := args[5]; else nt := 1; fi; if nt>NN then RETURN(FAIL); fi; flag := true; n := nops(Z); for i from 1 to n do interval[i] := [seq([(j-1)*B[i]/N, j*B[i]/N], j=1..N)]; od; Li := Nlists(n,N); #All lists of length n consisting of integers between 1 and N Bad := {}; #Will be the set of all bad hyperboxes for L in Li do Y := [seq(interval[j][L[j]],j=1..n)]; if nt=1 then print(SBM,Y); fi; H := subs({seq(Z[ii]=Z[ii]+Y[ii][1], ii=1..n)} , h); Hstar := subs({seq(Z[ii]=1/Z[ii], ii=1..n)}, H)*mul(Z[jj]^degree(H,Z[jj]), jj=1..n); Hstar0 := simplify(subs({seq(Z[ii]=Z[ii]+1/(Y[ii][2]-Y[ii][1]), ii=1..n)}, Hstar)); const := coeftayl(Hstar0, Z=[0$n], [0$n]); #the constant term in Hstar0 Pos := IsPosNoPrint(Hstar0,Z); if Pos=true then if const=0 then if not [seq(Y[i][2],i=1..n)]=B then if not [seq(Y[i][1],i=1..n)]=[0$n] then Bad := Bad union {Y}; flag := false; RETURN(false); else #keep flag what it is if nt=1 then print(true); fi; fi; else #keep flag what it is if nt=1 then print(true); fi; fi; else #keep flag what it is if nt=1 then print(true); fi; fi; elif Pos=false then Bad := Bad union {Y}; flag := false; if nt=1 then print(false); fi; else #Pos=FAIL then #don't add Y to Bad since we don't know if Hstar0 is < 0 on Y flag := false; if nt=1 then print(FAIL); fi; h2 := subs({seq(Z[ii]=Z[ii]+Y[ii][1],ii=1..n)},h); B2 := [seq(B[ii]-Y[ii][1],ii=1..n)]; if nt=1 then print(Running SBM on ,Y, with smaller subdivision.); fi; smaller := SmallBoxMethodNoPrint(h2,Z,B2,N,nt+1); if nt=1 then if smaller=true then print(On ,Y, we were redeemed.); elif smaller=false then print(FAIL turned into false on ,Y); elif smaller=FAIL then print(We still FAILed.); fi; fi; fi; od; if flag then RETURN(true); else if Bad={} then #flag was set to false but only because IsPos couldn't #determine the sign of Hstar0 for some Y RETURN(FAIL); else #flag was set to false becuase Hstar0 is <0 for some point #in each of the hyperboxes in Bad #print(Bad); RETURN(false); fi; fi; end: #==================================================================== #SmallBoxMethodNoPrintNoIter(h,Z,B,N) #Inputs: # h - a polynomial # Z - variables in h # B - list of positive real numbers of length nops(Z) # N - an integer #Outputs: # true if h(Z)>0 for 0<=Z[i]<=B[i] for all i, except possibly # h(0,...,0)>=0 and h(B)>=0 # false if there is a n assignment, Z', of variables s.t. h(Z')<0 # FAIL if the method is inconclusive #Method: # Splits the hyperbox 0<=Z[i]<=B[i] for all i into 2^N smaller hyperboxes # (split each axis into N intervals) and uses (the correct version of) # IsPos on each of the smaller hyperboxes. #Also: # "NoPrint" means that almost nothing is printed in this implementation of # SmallBoxMethod SmallBoxMethodNoPrintNoIter := proc(h,Z,B,N) local n,interval,i,j,L,Li,Y,H,Hstar,Hstar0,Os,C1,neg1,const,tryFlips,k,flag,flag2, disc,checkPos,tDegNeg,look,tDeg,ii,loc,a,b,c,Cd,solNumer,d,Bad,Pos; flag := true; n := nops(Z); for i from 1 to n do interval[i] := [seq([(j-1)*B[i]/N, j*B[i]/N], j=1..N)]; od; Li := Nlists(n,N); #All lists of length n consisting of integers between 1 and N Bad := {}; #Will be the set of all bad hyperboxes for L in Li do Y := [seq(interval[j][L[j]],j=1..n)]; print(SBM,Y); H := subs({seq(Z[ii]=Z[ii]+Y[ii][1], ii=1..n)} , h); Hstar := subs({seq(Z[ii]=1/Z[ii], ii=1..n)}, H)*mul(Z[jj]^degree(H,Z[jj]), jj=1..n); Hstar0 := simplify(subs({seq(Z[ii]=Z[ii]+1/(Y[ii][2]-Y[ii][1]), ii=1..n)}, Hstar)); const := coeftayl(Hstar0, Z=[0$n], [0$n]); #the constant term in Hstar0 Pos := IsPosNoPrint(Hstar0,Z); if Pos=true then if const=0 then if not [seq(Y[i][2],i=1..n)]=B then if not [seq(Y[i][1],i=1..n)]=[0$n] then Bad := Bad union {Y}; flag := false; RETURN(false); else #keep flag what it is print(true); fi; else #keep flag what it is print(true); fi; else #keep flag what it is print(true); fi; elif Pos=false then Bad := Bad union {Y}; flag := false; print(false); else #Pos=FAIL then #don't add Y to Bad since we don't know if Hstar0 is < 0 on Y flag := false; print(FAIL); fi; od; if flag then RETURN(true); else if Bad={} then #flag was set to false but only because IsPos couldn't #determine the sign of Hstar0 for some Y RETURN(FAIL); else #flag was set to false becuase Hstar0 is <0 for some point #in each of the hyperboxes in Bad #print(Bad); RETURN(false); fi; fi; end: #==================================================================== #PolynomialPositiveNoPrint(P,Z,xbar,N) #Inputs: # P - a polynomial # Z - variables in P # xbar - a positive real number which marks the first cut of Rplus^|Z| into 2^|Z| regions # N - an integer greater than 1. Will be used when calling SmallBoxMethod if necessary #Outputs: # Set containing elements from {true, false, FAIL}. # if {true} then # P(Z)>0 when Z[i]>=0 except possibly P(0,...,0)>=0 and P(xbar,...,xbar)>=0 # elif false in output then # there is a point where P(Z')<=0 when Z'[i]>=0 (or P(0,...,0)<0 or P(xbar,...,xbar)<0 ) # else (FAIL in output and false not in output) # there is some place where using IsPos is inconclusive #Method: # Cuts the positive orthant into 2^|Z| subsets and checks on each of # these whether P(Z)>0. Each subset, I, of {1,...,|Z|} corresponds to the # subset of the positive orthant where 0<=Z[i]<=xbar if i in I, # and xbar<=Z[j] if j not in I. #Also: # "NoPrint" means that almost nothing is printed in this implementation # of PolynomialPositive PolynomialPositiveNoPrint := proc(P,Z,xbar,N) local n,sets,J,f,g,C,neg,const,Result,h,compJ,B,i,K,j,g2,g3,flag,II, Os,C1,neg1,tDegNeg,look,tDeg,loc,checkPos,ii,disc,a,b,c,d,flag2,Cd,solNumer,Pos; n := nops(Z); sets := powerset(n); if xbar=0 then C := [coeffs(collect(P, Z, 'distributed'), Z)]; neg := select(i->evalb(i<=0) ,C); const := coeftayl(P, Z=[0$n], [0$n]); if nops(neg)=0 and const>0 then RETURN({true}); elif nops(neg)=0 and const=0 then RETURN({true}); else RETURN({false}); fi; fi; for J in sets do print(PolPos,J); f[J] := subs({seq(Z[j]=1/Z[j] , j in J)}, P)*mul(Z[j]^degree(P,Z[j]), j in J); compJ := {seq(i, i=1..n)} minus J; #{1,...,n} minus J g[J] := simplify(subs({seq(Z[j]=Z[j]+1/xbar, j in J), seq(Z[cj]=Z[cj]+xbar, cj in compJ)}, f[J])); Pos := IsPosNoPrint(g[J],Z); if Pos=true then Result[J] := true; elif Pos=false then Result[J] := false; else #Pos=FAIL h[J] := simplify(subs({seq(Z[j]=1/Z[j], j in compJ)}, P) *mul(Z[j]^degree(P,Z[j]), j in compJ)); B[J] := []; for i from 1 to n do if i in J then B[J] := [op(B[J]), xbar]; else B[J] := [op(B[J]),1/xbar]; fi; od; Result[J] := SmallBoxMethodNoPrint(h[J], Z, B[J], N); fi; od; for K in sets do print(K,Result[K]); od; RETURN({seq(Result[K], K in sets)}); end: #==================================================================== #ProveKNoPrint(R,vars,K,delta,N) #Inputs: # R - a rational difference equation # vars - variables in R # K - K value to test # delta - delta value to test # N - positive integer to be used when calling PolynomialPositive #Outputs: # List of lists of length 2. First element is equilibrium point # second element is output of PolynomialPositive for that equil. #Also: # "NoPrint" means that almost nothing is printed. ProveKNoPrint := proc(R,vars,K,delta,N) local P,p,GAS,newvars; newvars := [seq(z[i], i=1..nops(vars))]; P := OrigPoly(R,vars,newvars,K,delta); for p in P do GAS[p[1]] := PolynomialPositiveNoPrint(p[2],newvars,p[1],N); od; RETURN([seq([p[1], GAS[p[1]]], p in P)]); end: #==================================================================== #ProveNoPrint(R,vars,MinK,MaxK,delta,N) #Inputs: # R - A rational difference equation # vars - variables in R # MinK - the minimum K value to try # MaxK - the maximum K value to try # delta - the delta to try # N - positive integer to be used in PolynomialPositive #Outputs: # A set of lists of length 2. For each list the first element is # the equilibrium value, the second element is the K value for that equ. # If K=-1 then that equilibrium is not LAS. If K=0 then the MaxK # value is not high enough. #Method: # For each equilibrium calls PolynomialPositive with OrigPoly for # that equilibrium. #Also: # "NoPrint" means that almost nothing is printed ProveNoPrint := proc(R,vars,MinK,MaxK,delta,N) local newvars,P,P2,LS1,p,i,flag,Ret,K,GAS,l; newvars := [seq(z[i], i=1..nops(vars))]; LS1 := LinStab(R,vars); Ret := {}; for l in LS1 do flag := false; if l[2] then for K from MinK to MaxK while not flag do print(Testing K=,K, for the equilibrium xbar=,l[1]); P := OrigPoly(R,vars,newvars,K,delta); P2 := select(i->evalb(i[1]=l[1]),P)[1][2]; #print(P2); GAS[l] := PolynomialPositiveNoPrint(P2,newvars,l[1],N); print(GAS[l]); if GAS[l]={true} then flag := true; Ret := Ret union {[l[1],K]}; fi; od; if not flag then print(The MaxK for the equilibrium xbar=,l[1], is not large enough.); Ret := Ret union {[l[1],0]}; fi; else print(The equilibrium xbar=,l[1], is not LAS); Ret := Ret union {[l[1],-1]}; fi; od; RETURN(Ret); end: #==================================================================== #WebBookNoPrint(R,vars,params,paramPoss,numParams,MinK,MaxK,delta,N) #Inputs: # R - a rational difference equation # vars - variables in R # params - parameters in R # paramPoss - possible values for parameters # numParams - number of parameter sets to test # MinK - the minimum K value to test # MaxK - the maximum K value to test # delta - the delta to test # N - a positive integer to be used when calling PolynomialPositive #Outputs: # List of lists of length 2. First element is set of parameters, second # element is output of Prove for subs(params,R). #Method: # Chooses parameters from paramPoss for which the equilibrium is # rational if possible. If "not" possible (tries 10*numParams) then # checks other parameter sets. #Also: # "NoPrint" means that almost nothing is printed. WebBookNoPrint := proc(R,vars,params,paramPoss,numParams,MinK,MaxK,delta,N) local Ps,P,i,j,Ret,R1,k,K,r,r2,x,xK,C; Ps := {}; Ret := {}; for i from 1 to 10*numParams while nops(Ps)FAIL then print(P); Ps := Ps union {P}; R1 := subs(P,R); C := {coeffs(numer(R1),vars),coeffs(denom(R1),vars)}; if not member(-1,{seq(sign(C[i]),i=1..nops(C))}) then Ret := Ret union {[P,ProveNoPrint(R1,vars,MinK,MaxK,delta,N)]}; else Ps := Ps minus {P}; fi; else fi; od; if nops(Ps)0 when all variables are >=0 (except H(Z)>=0 when all vars=0) # false if there is an assignment, Z', of the variables in the positive # orthant s.t. H(Z')<=0 (or H(0)<0) # FAIL if the sufficient conditions checked are inconclusive #Method: # Checks two sufficient conditions: # 1. All coefficients are positive, and the constant term is positive or 0. # If the constant term is 0 then check that there is no subset of the # variables s.t. setting that subset 0 yields the 0 polynomial. # 2. If the only negative coefficient is on a term from the set # {Z[i]*Z[j], Z[i], Z[i]^2} then checks that the sub-polynomial # consisting of the terms from this set (for all i,j) is still positive # using the quadratic formula. #Also: # "NoPrint" means that nothing is printed in this implementation of IsPos IsPosNoPrint := proc(H,Z) local Os,C1,neg1,const,n,sets,flag,II,h,LargestNegDeg,tDeg,look,ii,loc, checkPos,flag2,a,b,c,disc,Cd,j,sets2,flag3,JJ,disc0,CheckPos0,d; n := nops(Z); sets := powerset(n); Os := [op(sort(collect(H, Z, 'distributed'),order=tdeg(op(Z))))]; #list of terms in H C1 := [seq([coeffs(Os[i],Z),Os[i]/coeffs(Os[i],Z),i],i=1..nops(Os))]; #list of [coeff,vars,term#] neg1 := select(i->evalb(evalf(i[1])<=0),C1); #list of terms w/ negative coeff const := coeftayl(H, Z=[0$n], [0$n]); #constant term in H if nops(neg1)=0 and evalf(const)>0 then #all coeffs positive and constant term positive RETURN(true); elif nops(neg1)=0 and const=0 then #all coeffs positive and constant term = 0 for II in sets minus {{seq(i,i=1..n)}} do h[II] := subs({seq(Z[i]=0, i in II)},H); if h[II]=0 then RETURN(false); fi; od; RETURN(true); elif evalf(const)<0 then #constant term is negative RETURN(false); else #there are negative coefficients and constant term is >= 0 if max({seq(degree(neg1[1][2],Z[i]),i=1..n)})>2 then RETURN(FAIL); else LargestNegDeg := add(degree(neg1[1][2],Z[i]),i=1..n); look := neg1[1][3]; #the index in Os of the first term w/ neg. coeff loc := look; #if look=1 then this is needed because we won't enter into the for loop below tDeg := LargestNegDeg; #will provide ending condition for the for loop for ii from 1 to look-1 while tDeg = LargestNegDeg do tDeg := add(degree(C1[look-ii][2],Z[j]), j=1..n); loc := C1[look-ii][3]; #will be the position in Os from which to start checkPos od; if LargestNegDeg=tDeg then #we got all the way to look-1 with the same total degree checkPos := add(Os[j], j=loc..nops(Os)); else #we stopped before ii=look-1 so we must add 1 to loc checkPos := add(Os[j], j=loc+1..nops(Os)); fi; flag2 := true; #will serve as "flag" in the pseudocode algorithm for j from 1 to n while flag2 do if degree(checkPos,Z[j])=2 then a := coeff(checkPos,Z[j],2); b := coeff(checkPos,Z[j],1); c := coeff(checkPos,Z[j],0); disc := b^2-4*a*c; Cd := [coeffs(collect(disc,Z,'distributed'),Z)]; #list of coeffs in disc if {seq(sign(evalf(Cd[i])), i=1..nops(Cd))}={-1} then #all coeffs are negative sets2 := powerset({seq(i, i=1..j-1), seq(i,i=j+1..n)}) minus {{seq(i, i=1..j-1), seq(i,i=j+1..n)}}; flag3 := true; for JJ in sets2 while flag3 do d[JJ] := subs({seq(Z[i]=0, i in JJ)},disc); if d[JJ]=0 then flag3 := false; fi; od; if flag3 then #there are no subsets of {1,...,j-1,j+1,...,n} s.t. setting those vars to 0 yields disc=0 disc0 := subs({seq(Z[i]=0, i=1..j-1),seq(Z[i]=0,i=j+1..n)},disc); CheckPos0 := subs({seq(Z[i]=0, i=1..j-1),seq(Z[i]=0,i=j+1..n)},checkPos); if disc0 < 0 or (disc0 = 0 and coeff(CheckPos0,Z[j])*coeff(CheckPos0,Z[j],2)>=0) then #keep flag=true else flag2 := false; fi; else flag2 := false; fi; else #there are positive coefficients in disc flag2 := false; fi; else #degree(checkPos,Z[j]) <> 2 flag2 := false; fi; od; if flag2 then RETURN(true) fi; fi; RETURN(FAIL); fi; end: #==================================================================== #SmallBoxMethodAbsNoPrint(h,Z,B,N) #Inputs: # h - a polynomial # Z - variables in h # B - list of positive real numbers of length nops(Z) # N - an integer #Outputs: # true if h(Z)>0 for 0<=Z[i]<=B[i] for all i, except possibly # h(0,...,0)>=0 and h(B)>=0 # false if there is a n assignment, Z', of variables s.t. h(Z')<0 # FAIL if the method is inconclusive #Method: # Splits the hyperbox 0<=Z[i]<=B[i] for all i into 2^N smaller hyperboxes # (split each axis into N intervals) and uses (the correct version of) # IsPos on each of the smaller hyperboxes. #Also: # "AbsNoPrint" means that nothing is printed in this implementation # of SmallBoxMethod SmallBoxMethodAbsNoPrint := proc(h,Z,B,N) local n,interval,i,j,L,Li,Y,H,Hstar,Hstar0,Os,C1,neg1,const,tryFlips,k,flag,flag2, disc,checkPos,tDegNeg,look,tDeg,ii,loc,a,b,c,Cd,solNumer,d,Bad,Pos; flag := true; n := nops(Z); for i from 1 to n do interval[i] := [seq([(j-1)*B[i]/N, j*B[i]/N], j=1..N)]; od; Li := Nlists(n,N); #All lists of length n consisting of integers between 1 and N Bad := {}; #Will be the set of all bad hyperboxes for L in Li do Y := [seq(interval[j][L[j]],j=1..n)]; #print(SBM,Y); H := subs({seq(Z[ii]=Z[ii]+Y[ii][1], ii=1..n)} , h); Hstar := subs({seq(Z[ii]=1/Z[ii], ii=1..n)}, H)*mul(Z[jj]^degree(H,Z[jj]), jj=1..n); Hstar0 := simplify(subs({seq(Z[ii]=Z[ii]+1/(Y[ii][2]-Y[ii][1]), ii=1..n)}, Hstar)); const := coeftayl(Hstar0, Z=[0$n], [0$n]); #the constant term in Hstar0 Pos := IsPosNoPrint(Hstar0,Z); if Pos=true then if const=0 then if not [seq(Y[i][2],i=1..n)]=B then if not [seq(Y[i][1],i=1..n)]=[0$n] then Bad := Bad union {Y}; flag := false; #print(false); else #keep flag what it is #print(true); fi; else #keep flag what it is #print(true); fi; else #keep flag what it is #print(true); fi; elif Pos=false then Bad := Bad union {Y}; flag := false; #print(false); else #Pos=FAIL then #don't add Y to Bad since we don't know if Hstar0 is < 0 on Y flag := false; #print(FAIL); fi; od; if flag then RETURN(true); else if Bad={} then #flag was set to false but only because IsPos couldn't #determine the sign of Hstar0 for some Y RETURN(FAIL); else #flag was set to false becuase Hstar0 is <0 for some point #in each of the hyperboxes in Bad #print(Bad); RETURN(false); fi; fi; end: #==================================================================== #PolynomialPositiveAbsNoPrint(P,Z,xbar,N)); #Inputs:); # P - a polynomial # Z - variables in P # xbar - a positive real number which marks the first cut of Rplus^|Z| into 2^|Z| regions # N - an integer greater than 1. Will be used when calling SmallBoxMethod if necessary #Outputs: # Set containing elements from {true, false, FAIL}. # if {true} then # P(Z)>0 when Z[i]>=0 except possibly P(0,...,0)>=0 and P(xbar,...,xbar)>=0 # elif false in output then # there is a point where P(Z')<=0 when Z'[i]>=0 (or P(0,...,0)<0 or P(xbar,...,xbar)<0 ) # else (FAIL in output and false not in output) # there is some place where using IsPos is inconclusive #Method: # Cuts the positive orthant into 2^|Z| subsets and checks on each of # these whether P(Z)>0. Each subset, I, of {1,...,|Z|} corresponds to the # subset of the positive orthant where 0<=Z[i]<=xbar if i in I, # and xbar<=Z[j] if j not in I. #Also: # "AbsNoPrint" means that absolutely nothing is printed in this implementation # of PolynomialPositive PolynomialPositiveAbsNoPrint := proc(P,Z,xbar,N) local n,sets,J,f,g,C,neg,const,Result,h,compJ,B,i,K,j,g2,g3,flag,II, Os,C1,neg1,tDegNeg,look,tDeg,loc,checkPos,ii,disc,a,b,c,d,flag2,Cd,solNumer,Pos; n := nops(Z); sets := powerset(n); #print(sets); if xbar=0 then C := [coeffs(collect(P, Z, 'distributed'), Z)]; neg := select(i->evalb(i<=0) ,C); const := coeftayl(P, Z=[0$n], [0$n]); if nops(neg)=0 and const>0 then RETURN({true}); elif nops(neg)=0 and const=0 then RETURN({true}); else RETURN({false}); fi; fi; for J in sets do #print(PolPos,J); f[J] := subs({seq(Z[j]=1/Z[j] , j in J)}, P)*mul(Z[j]^degree(P,Z[j]), j in J); compJ := {seq(i, i=1..n)} minus J; #{1,...,n} minus J g[J] := simplify(subs({seq(Z[j]=Z[j]+1/xbar, j in J), seq(Z[cj]=Z[cj]+xbar, cj in compJ)}, f[J])); Pos := IsPosNoPrint(g[J],Z); if Pos=true then Result[J] := true; elif Pos=false then Result[J] := false; else #Pos=FAIL h[J] := simplify(subs({seq(Z[j]=1/Z[j], j in compJ)}, P) *mul(Z[j]^degree(P,Z[j]), j in compJ)); B[J] := []; for i from 1 to n do if i in J then B[J] := [op(B[J]), xbar]; else B[J] := [op(B[J]),1/xbar]; fi; od; Result[J] := SmallBoxMethodAbsNoPrint(h[J], Z, B[J], N); fi; od; for K in sets do #print(K,Result[K]); od; RETURN({seq(Result[K], K in sets)}); end: #==================================================================== #IsPosNoPrint(H,Z) #Inputs: # H - a polynomial # Z - variables in H #Outputs: # true if H(Z)>0 when all variables are >=0 (except H(Z)>=0 when all vars=0) # false if there is an assignment, Z', of the variables in the positive # orthant s.t. H(Z')<=0 (or H(0)<0) # FAIL if the sufficient conditions checked are inconclusive #Method: # Checks two sufficient conditions: # 1. All coefficients are positive, and the constant term is positive or 0. # If the constant term is 0 then check that there is no subset of the # variables s.t. setting that subset 0 yields the 0 polynomial. # 2. If the only negative coefficient is on a term from the set # {Z[i]*Z[j], Z[i], Z[i]^2} then checks that the sub-polynomial # consisting of the terms from this set (for all i,j) is still positive # using the quadratic formula. #Also: # "NoPrint" means that nothing is printed in this implementation of IsPos IsPosNoPrint := proc(H,Z) local Os,C1,neg1,const,n,sets,flag,II,h,LargestNegDeg,tDeg,look,ii,loc, checkPos,flag2,a,b,c,disc,Cd,j,sets2,flag3,JJ,disc0,CheckPos0,d; n := nops(Z); sets := powerset(n); Os := [op(sort(collect(H, Z, 'distributed'),order=tdeg(op(Z))))]; #list of terms in H C1 := [seq([coeffs(Os[i],Z),Os[i]/coeffs(Os[i],Z),i],i=1..nops(Os))]; #list of [coeff,vars,term#] neg1 := select(i->evalb(evalf(i[1])<=0),C1); #list of terms w/ negative coeff const := coeftayl(H, Z=[0$n], [0$n]); #constant term in H if nops(neg1)=0 and evalf(const)>0 then #all coeffs positive and constant term positive RETURN(true); elif nops(neg1)=0 and const=0 then #all coeffs positive and constant term = 0 for II in sets minus {{seq(i,i=1..n)}} do h[II] := subs({seq(Z[i]=0, i in II)},H); if h[II]=0 then RETURN(false); fi; od; RETURN(true); elif evalf(const)<0 then #constant term is negative RETURN(false); else #there are negative coefficients and constant term is >= 0 if max({seq(degree(neg1[1][2],Z[i]),i=1..n)})>2 then RETURN(FAIL); else LargestNegDeg := add(degree(neg1[1][2],Z[i]),i=1..n); look := neg1[1][3]; #the index in Os of the first term w/ neg. coeff loc := look; #if look=1 then this is needed because we won't enter into the for loop below tDeg := LargestNegDeg; #will provide ending condition for the for loop for ii from 1 to look-1 while tDeg = LargestNegDeg do tDeg := add(degree(C1[look-ii][2],Z[j]), j=1..n); loc := C1[look-ii][3]; #will be the position in Os from which to start checkPos od; if LargestNegDeg=tDeg then #we got all the way to look-1 with the same total degree checkPos := add(Os[j], j=loc..nops(Os)); else #we stopped before ii=look-1 so we must add 1 to loc checkPos := add(Os[j], j=loc+1..nops(Os)); fi; flag2 := true; #will serve as "flag" in the pseudocode algorithm for j from 1 to n while flag2 do if degree(checkPos,Z[j])=2 then a := coeff(checkPos,Z[j],2); b := coeff(checkPos,Z[j],1); c := coeff(checkPos,Z[j],0); disc := b^2-4*a*c; Cd := [coeffs(collect(disc,Z,'distributed'),Z)]; #list of coeffs in disc if {seq(sign(evalf(Cd[i])), i=1..nops(Cd))}={-1} then #all coeffs are negative sets2 := powerset({seq(i, i=1..j-1), seq(i,i=j+1..n)}) minus {{seq(i, i=1..j-1), seq(i,i=j+1..n)}}; flag3 := true; for JJ in sets2 while flag3 do d[JJ] := subs({seq(Z[i]=0, i in JJ)},disc); if d[JJ]=0 then flag3 := false; fi; od; if flag3 then #there are no subsets of {1,...,j-1,j+1,...,n} s.t. setting those vars to 0 yields disc=0 disc0 := subs({seq(Z[i]=0, i=1..j-1),seq(Z[i]=0,i=j+1..n)},disc); CheckPos0 := subs({seq(Z[i]=0, i=1..j-1),seq(Z[i]=0,i=j+1..n)},checkPos); if disc0 < 0 or (disc0 = 0 and coeff(CheckPos0,Z[j])*coeff(CheckPos0,Z[j],2)>=0) then #keep flag=true else flag2 := false; fi; else flag2 := false; fi; else #there are positive coefficients in disc flag2 := false; fi; else #degree(checkPos,Z[j]) <> 2 flag2 := false; fi; od; if flag2 then RETURN(true) fi; fi; RETURN(FAIL); fi; end: #==================================================================== #SmallBoxMethodParamsNoPrint(h,Z,params,B,N) #Inputs: # h - a polynomial # Z - variables in h # params - parameters in h # B - list of positive real numbers of length nops(Z) # N - an integer #Outputs: # true if h(Z)>0 for 0<=Z[i]<=B[i] for all i, and all params > 0, except possibly # h(0,...,0)>=0 and h(B)>=0 # false if there is a n assignment, Z', of variables s.t. h(Z')<0 # FAIL if the method is inconclusive #Method: # Splits the hyperbox 0<=Z[i]<=B[i] for all i into 2^N smaller hyperboxes # (split each axis into N intervals) and uses (the correct version of) # IsPos (with Z union params as variables) on each of the smaller hyperboxes. #Also: # "ParamsNoPrint" means that almost nothing is printed in this implementation # of SmallBoxMethod and that parameters are treated differently than variables SmallBoxMethodParamsNoPrint := proc(h,Z,params,B,N) local n,interval,i,j,L,Li,Y,H,Hstar,Hstar0,Os,C1,neg1,const,tryFlips,k,flag,flag2, disc,checkPos,tDegNeg,look,tDeg,ii,loc,a,b,c,Cd,solNumer,d,Bad,Pos,ZP; flag := true; n := nops(Z); ZP := [op(Z),op(params)]; for i from 1 to n do interval[i] := [seq([(j-1)*B[i]/N, j*B[i]/N], j=1..N)]; od; Li := Nlists(n,N); #All lists of length n consisting of integers between 1 and N Bad := {}; #Will be the set of all bad hyperboxes for L in Li do Y := [seq(interval[j][L[j]],j=1..n)]; print(SBM,Y); H := subs({seq(Z[ii]=Z[ii]+Y[ii][1], ii=1..n)} , h); Hstar := subs({seq(Z[ii]=1/Z[ii], ii=1..n)}, H)*mul(Z[jj]^degree(H,Z[jj]), jj=1..n); Hstar0 := simplify(subs({seq(Z[ii]=Z[ii]+1/(Y[ii][2]-Y[ii][1]), ii=1..n)}, Hstar)); const := coeftayl(Hstar0, Z=[0$n], [0$n]); #the constant term in Hstar0 Pos := IsPosNoPrint(Hstar0,ZP); if Pos=true then if const=0 then if not [seq(Y[i][2],i=1..n)]=B then if not [seq(Y[i][1],i=1..n)]=[0$n] then Bad := Bad union {Y}; flag := false; print(false); else #keep flag what it is print(true,1); fi; else #keep flag what it is print(true,2); fi; else #keep flag what it is print(true,3); fi; elif Pos=false then Bad := Bad union {Y}; flag := false; print(false); else #Pos=FAIL then #don't add Y to Bad since we don't know if Hstar0 is < 0 on Y flag := false; #print(Hstar0); print(FAIL); fi; od; if flag then RETURN(true); else if Bad={} then #flag was set to false but only because IsPos couldn't #determine the sign of Hstar0 for some Y RETURN(FAIL); else #flag was set to false becuase Hstar0 is <0 for some point #in each of the hyperboxes in Bad #print(Bad); RETURN(false); fi; fi; end: #==================================================================== #PolynomialPositiveParamsNoPrint(P,Z,params,xbar,N)); #Inputs:); # P - a polynomial # Z - variables in P # params - parameters in P # xbar - a positive real number which marks the first cut of Rplus^|Z| into 2^|Z| regions # N - an integer greater than 1. Will be used when calling SmallBoxMethod if necessary #Outputs: # Set containing elements from {true, false, FAIL}. # if {true} then # P(Z)>0 when Z[i]>=0 and params >=0 except possibly # P(0,...,0)>=0 and P(xbar,...,xbar)>=0 # elif false in output then # there is a point where P(Z')<=0 when Z'[i]>=0 # (or P(0,...,0)<0 or P(xbar,...,xbar)<0 ) # else (FAIL in output and false not in output) # there is some place where using IsPos is inconclusive #Method: # Cuts the positive orthant into 2^|Z| subsets and checks on each of # these whether P(Z)>0. Each subset, I, of {1,...,|Z|} corresponds to the # subset of the positive orthant where 0<=Z[i]<=xbar if i in I, # and xbar<=Z[j] if j not in I. Treats parameters differently. No cutting # of parameter space, always just assuming the params are <= 0. #Also: # "ParamsNoPrint" means that almost nothing is printed in this implementation # of PolynomialPositive and that parameters are treated differently than variables PolynomialPositiveParamsNoPrint := proc(P,Z,params,xbar,N) local n,sets,J,f,g,C,neg,const,Result,h,compJ,B,i,K,j,g2,g3,flag,II,ZP,np, Os,C1,neg1,tDegNeg,look,tDeg,loc,checkPos,ii,disc,a,b,c,d,flag2,Cd,solNumer,Pos; n := nops(Z); sets := powerset(n); ZP := [op(Z),op(params)]; np := nops(ZP); if xbar=0 then C := [coeffs(collect(P, ZP, 'distributed'), ZP)]; neg := select(i->evalb(i<=0) ,C); const := coeftayl(P, ZP=[0$np], [0$np]); if nops(neg)=0 and const>0 then RETURN({true}); elif nops(neg)=0 and const=0 then RETURN({true}); else RETURN({false}); fi; fi; for J in sets do print(PolPos,J); f[J] := subs({seq(Z[j]=1/Z[j] , j in J)}, P)*mul(Z[j]^degree(P,Z[j]), j in J); compJ := {seq(i, i=1..n)} minus J; #{1,...,n} minus J g[J] := simplify(subs({seq(Z[j]=Z[j]+1/xbar, j in J), seq(Z[cj]=Z[cj]+xbar, cj in compJ)}, f[J])); Pos := IsPosNoPrint(g[J],ZP); if Pos=true then Result[J] := true; elif Pos=false then Result[J] := false; else #Pos=FAIL h[J] := simplify(subs({seq(Z[j]=1/Z[j], j in compJ)}, P) *mul(Z[j]^degree(P,Z[j]), j in compJ)); B[J] := []; for i from 1 to n do if i in J then B[J] := [op(B[J]), xbar]; else B[J] := [op(B[J]),1/xbar]; fi; od; Result[J] := SmallBoxMethodParamsNoPrint(h[J], Z, params, B[J], N); fi; od; for K in sets do print(K,Result[K]); od; RETURN({seq(Result[K], K in sets)}); end: #==================================================================== #IsPos(H,Z) #Inputs: # H - a polynomial # Z - variables in H #Outputs: # true if H(Z)>0 when all variables are >=0 (except H(Z)>=0 when all vars=0) # false if there is an assignment, Z', of the variables in the positive # orthant s.t. H(Z')<=0 (or H(0)<0) # FAIL if the sufficient conditions checked are inconclusive #Method: # Checks two sufficient conditions: # 1. All coefficients are positive, and the constant term is positive or 0. # If the constant term is 0 then check that there is no subset of the # variables s.t. setting that subset 0 yields the 0 polynomial. # 2. If the only negative coefficient is on a term from the set # {Z[i]*Z[j], Z[i], Z[i]^2} then checks that the sub-polynomial # consisting of the terms from this set (for all i,j) is still positive # using the quadratic formula. #Also: # IsPos prints what it is doing as we go along. IsPos := proc(H,Z) local Os,C1,neg1,const,n,sets,flag,II,h,LargestNegDeg,tDeg,look,ii,loc, checkPos,flag2,a,b,c,disc,Cd,j,sets2,flag3,JJ,disc0,CheckPos0,d,Pos,pows; n := nops(Z); sets := powerset(n); Os := [op(sort(collect(H, Z, 'distributed'),order=tdeg(op(Z))))]; #list of terms in H C1 := [seq([coeffs(Os[i],Z),Os[i]/coeffs(Os[i],Z),i],i=1..nops(Os))]; #list of [coeff,vars,term#] neg1 := select(i->evalb(evalf(i[1])<=0),C1); #list of terms w/ negative coeff const := coeftayl(H, Z=[0$n], [0$n]); #constant term in H if nops(neg1)=0 and const>0 then #all coeffs positive and constant term positive ########## print(All coefficients are positive, and the constant term is positive,); print(so the polynomial is positive.); ########## RETURN(true); elif nops(neg1)=0 and const=0 then #all coeffs positive and constant term = 0 for II in sets minus {{seq(i,i=1..n)}} do h[II] := subs({seq(Z[i]=0, i in II)},H); if h[II]=0 then ########## print(If we set the variables , seq(Z[i],i in II), equal 0); print(then the polynomial equals 0 which is bad.); ########## RETURN(false); fi; od; ########## print(All coefficients are positive, and the constant term is 0.); print(Also, there is no proper subset of the variables); print(for which setting them all equal 0 yields the 0 polynomial.); print(Therefore, the polynomial is 0 only when all variables are 0.); ########## RETURN(true); elif evalf(const)<0 then #constant term is negative ########## print(The constant term in the polynomial is negative.); print(Therefore, the polynomial is negative when all variables are 0.); ########## RETURN(false); else #there are negative coefficients and constant term is >= 0 ########## print(There are negative coefficients, they are:); print(neg1); ########## if max({seq(degree(neg1[1][2],Z[i]),i=1..n)})>2 then ########## print(The degree of the largest negative term is greater than 2 in one of the variables,); print(so we can't use the discriminant method. We can't tell the sign of the polynomial for sure); ########## RETURN(FAIL); else ########## print(We must check that the portion of the polynomial with total degree); print(less than or equal to that of the "highest" negative coefficient term); print(is actually positive. Then, since the rest of the polynomial has no); print(negative coefficients, the whole polynomial is positive.); ########## ########## #print(The degree of the largest negative term is less than or equal to 2 in all of the); #print(variables. So we must check that the portion of the polynomial in which each variable has); #print(degree less than or equal to 2 is actually positive. Then, since the rest of); #print(the polynomial has no negative coefficients, the whole polynomial is positive.); ########## LargestNegDeg := add(degree(neg1[1][2],Z[i]),i=1..n); look := neg1[1][3]; #the index in Os of the first term w/ neg. coeff loc := look; #if look=1 then this is needed because we won't enter into the for loop below tDeg := LargestNegDeg; #will provide ending condition for the for loop for ii from 1 to look-1 while tDeg = LargestNegDeg do tDeg := add(degree(C1[look-ii][2],Z[j]), j=1..n); loc := C1[look-ii][3]; #will be the position in Os from which to start checkPos od; if LargestNegDeg=tDeg then #we got all the way to look-1 with the same total degree checkPos := add(Os[j], j=loc..nops(Os)); else #we stopped before ii=look-1 so we must add 1 to loc checkPos := add(Os[j], j=loc+1..nops(Os)); fi; #pows := Nlists0(n,2); #checkPos := add(coeftayl(H,Z=[0$n],l)*mul(Z[i]^l[i],i=1..n),l in pows); ########## print(The polynomial with negative coefficients that we want to check is positive is:); print(CP = , checkPos); ########## flag2 := true; #will serve as "flag" in the pseudocode algorithm for j from 1 to n while flag2 do if degree(checkPos,Z[j])=2 then ########## print(The degree of CP in the variable ,Z[j], is 2, so we check to see if the discriminant is negative.); ########## a := coeff(checkPos,Z[j],2); b := coeff(checkPos,Z[j],1); c := coeff(checkPos,Z[j],0); disc := b^2-4*a*c; ########## print(The discriminant is Disc=,disc); ########## Cd := [coeffs(collect(disc,Z,'distributed'),Z)]; #list of coeffs in disc if {seq(sign(evalf(Cd[i])), i=1..nops(Cd))}={-1} then #all coeffs are negative ########## #print(According to IsPos, -disc >0); print(All coefficients in the discriminant are negative.); ########## sets2 := powerset({seq(i, i=1..j-1), seq(i,i=j+1..n)}) minus {{seq(i, i=1..j-1), seq(i,i=j+1..n)}}; flag3 := true; for JJ in sets2 while flag3 do d[JJ] := subs({seq(Z[i]=0, i in JJ)},disc); if d[JJ]=0 then ########## print(If we set the variables , seq(Z[i],i in JJ), equal 0); print(then the discriminant equals 0 which is bad.); ########## flag3 := false; fi; od; if flag3 then #there are no subsets of {1,...,j-1,j+1,...,n} s.t. setting those vars to 0 yields disc=0 ########## print(There is no proper subset of the variables which, when set to 0, yields Disc=0); ########## disc0 := subs({seq(Z[i]=0, i=1..j-1),seq(Z[i]=0,i=j+1..n)},disc); CheckPos0 := subs({seq(Z[i]=0, i=1..j-1),seq(Z[i]=0,i=j+1..n)},checkPos); if disc0 < 0 or (disc0 = 0 and coeff(CheckPos0,Z[j])*coeff(CheckPos0,Z[j],2)>=0) then ########## print(Setting all variables except ,Z[j], equal 0 in Disc yields ,disc0); print(Setting all varialbes except ,Z[j], equal 0 in CP yields ,CheckPos0); ########## #keep flag=true else ########## print(Setting all variables except ,Z[j], equal 0 in Disc yields ,disc0); print(Setting all varialbes except ,Z[j], equal 0 in CP yields ,CheckPos0); ########## flag2 := false; fi; else flag2 := false; fi; else #there are positive coefficients in disc ########## print(There are some positive coefficients in Disc, so we can't tell if the discriminant is negative.); ########## flag2 := false; fi; else #degree(checkPos,Z[j]) <> 2 ########## print(The degree of CP is <>2 in the variable, Z[j], so we can't use); print(the discriminant method.); ########## flag2 := false; fi; od; if flag2 then ########## print(CP is positive when all variables are non-negative, and CP may equal 0 when all ); print(variables are 0. Therefore the polynomial is positive when all variables are ); print(non-negative and may equal 0 when all variables are 0.); ########## RETURN(true) fi; fi; RETURN(FAIL); fi; end: #==================================================================== #SmallBoxMethod(h,Z,B,N) #Inputs: # h - a polynomial # Z - variables in h # B - list of positive real numbers of length nops(Z) # N - an integer #Outputs: # true if h(Z)>0 for 0<=Z[i]<=B[i] for all i, except possibly # h(0,...,0)>=0 and h(B)>=0 # false if there is a n assignment, Z', of variables s.t. h(Z')<0 # FAIL if the method is inconclusive #Method: # Splits the hyperbox 0<=Z[i]<=B[i] for all i into 2^N smaller hyperboxes # (split each axis into N intervals) and uses (the correct version of) # IsPos on each of the smaller hyperboxes. #Also: # There are print statements throughout. SmallBoxMethod := proc(h,Z,B,N) local n,interval,i,j,L,Li,Y,H,Hstar,Hstar0,Os,C1,neg1,const,tryFlips,k,flag,flag2, disc,checkPos,tDegNeg,look,tDeg,ii,loc,a,b,c,Cd,solNumer,d,Bad,Pos; ########## print(Starting the procedure SmallBoxMethod with inputs:); #print(h=, h); print(Z=, Z); print(B=, B); print(N=, N); ########## flag := true; n := nops(Z); for i from 1 to n do interval[i] := [seq([(j-1)*B[i]/N, j*B[i]/N], j=1..N)]; od; Li := Nlists(n,N); #All lists of length n consisting of integers between 1 and N Bad := {}; #Will be the set of all bad hyperboxes for L in Li do Y := [seq(interval[j][L[j]],j=1..n)]; ########## print(Trying to prove h>0 for ,Z, in the region:); for k from 1 to n do print(Y[k][1], < , Z[k], <= , Y[k][2]); od; print(Make a new polynomial, H, by substituting into h:); for k from 1 to n do print(Z[k]=Z[k]+Y[k][1]); od; ########## H := simplify(subs({seq(Z[ii]=Z[ii]+Y[ii][1], ii=1..n)} , h)); ########## print(Make a new polynomial, Hstar, by substituting into H:); for k from 1 to n do print(Z[k]=1/Z[k]); od; print(and then multiplying by); for k from 1 to n do print(Z[k],^degree(H,Z[k],)=,Z[k]^degree(H,Z[k])); od; ########## Hstar := subs({seq(Z[ii]=1/Z[ii], ii=1..n)}, H)*mul(Z[jj]^degree(H,Z[jj]), jj=1..n); ########## print(Make a new polynomial, Hstar0, by substituting into Hstar:); for k from 1 to n do print(Z[k]=Z[k]+1/(Y[k][2]-Y[k][1])); od; ########## Hstar0 := simplify(subs({seq(Z[ii]=Z[ii]+1/(Y[ii][2]-Y[ii][1]), ii=1..n)}, Hstar)); #Hstar0 := subs({seq(Z[ii]=Z[ii]+1/(Y[ii][2]-Y[ii][1]), ii=1..n)}, Hstar); ########## print(Need to prove that Hstar0(,op(Z),) > 0 for all variables >= 0); ########## const := coeftayl(Hstar0, Z=[0$n], [0$n]); #the constant term in Hstar0 Pos := IsPos(Hstar0,Z); if Pos=true then if const=0 then if not [seq(Y[i][2],i=1..n)]=B then if not [seq(Y[i][1],i=1..n)]=[0$n] then ########## print(IsPos returned true in this hyperbox, but the constant term); print(is 0, and the box is not the furthest from the origin or closest to the origin.); print(Therefore, there is a point where h=0 that we don't want h=0.); print(Y,false); ########## Bad := Bad union {Y}; flag := false; RETURN(false); else ########## print(The polynomial is positive on the hyperbox except at the origin:); print(Y,true); ########## #keep flag what it is #continue to examine the hyperboxes fi; else ########## print(The polynomial is positive on the hyperbox:); print(Y,true); ########## #keep flag what it is #continue to examine the hyperboxes fi; else ########## print(The polynomial is positive on the hyperbox:); print(Y,true); ########## #keep flag what it is #continue to examine the hyperboxes fi; elif Pos=false then ########## print(IsPos returned false, so h is not positive on the hyperbox:); print(Y,false); ########## Bad := Bad union {Y}; flag := false; RETURN(false); else #Pos=FAIL then ########## print(IsPos was inconclusive, so the sign of h is not known on the hyperbox:); print(Y,FAIL); ########## #don't add Y to Bad since we don't know if Hstar0 is < 0 on Y flag := false; fi; print(-------------------------------------------------------------); od; if flag then RETURN(true); else if Bad={} then #flag was set to false but only because IsPos couldn't #determine the sign of Hstar0 for some Y RETURN(FAIL); else ########## print(Bad); ########## #flag was set to false becuase Hstar0 is <0 for some point #in each of the hyperboxes in Bad RETURN(false); fi; fi; end: #==================================================================== #PolynomialPositive(P,Z,xbar,N) #Inputs: # P - a polynomial # Z - variables in P # xbar - a positive real number which marks the first cut of Rplus^|Z| into 2^|Z| regions # N - an integer greater than 1. Will be used when calling SmallBoxMethod if necessary #Outputs: # Set containing elements from {true, false, FAIL}. # if {true} then # P(Z)>0 when Z[i]>=0 except possibly P(0,...,0)>=0 and P(xbar,...,xbar)>=0 # elif false in output then # there is a point where P(Z')<=0 when Z'[i]>=0 (or P(0,...,0)<0 or P(xbar,...,xbar)<0 ) # else (FAIL in output and false not in output) # there is some place where using IsPos is inconclusive #Method: # Cuts the positive orthant into 2^|Z| subsets and checks on each of # these whether P(Z)>0. Each subset, I, of {1,...,|Z|} corresponds to the # subset of the positive orthant where 0<=Z[i]<=xbar if i in I, # and xbar<=Z[j] if j not in I. #Also: # There are print statements in this implementation of PolynomialPositive #Try: # PolynomialPositive(x^2-x*y+y^2,[x,y],1,2); PolynomialPositive := proc(P,Z,xbar,N) local n,sets,J,f,g,C,neg,const,Result,h,compJ,B,i,K,j,g2,g3,flag,II, Os,C1,neg1,tDegNeg,look,tDeg,loc,checkPos,ii,disc,a,b,c,d,flag2,Cd,solNumer,Pos; ##### print(Starting the procedure PolynomialPositive with inputs:); print(P=,P); print(Z=,Z); print(xbar=,xbar); print(N=,N); ##### n := nops(Z); sets := powerset(n); if xbar=0 then ########## print(Since xbar=0 we check that all coefficients); print(and the constant term in P are positive); ########## C := [coeffs(collect(P, Z, 'distributed'), Z)]; neg := select(i->evalb(i<=0) ,C); const := coeftayl(P, Z=[0$n], [0$n]); if nops(neg)=0 and const>0 then ########## print(All coefficients in P and the constant term are positive,); print(so P>0 for all variables >=0); ########## RETURN({true}); elif nops(neg)=0 and const=0 then ########## print(All coefficients in P are positive, and the constant term is 0); print(So P is >=0 when all variables are >=0); ########## RETURN({true}); else ########## if nops(neg)>0 then print(The negative coefficients are ,neg); fi; if const<=0 then print(The constant term is ,const); fi; print(so this method doesn't prove that P>0); ########## RETURN({false}); fi; fi; for J in sets do print(J); ########## print(Proving P>0 in the region:); for j from 1 to n do if j in J then print(0, < , Z[j], <= , xbar); else print(xbar, <= , Z[j],  < , infinity); fi; od; if J<>{} then print(Make new polynomial f[,J,] by substituting ); for j from 1 to n do if j in J then print(Z[j]=1/Z[j]); fi; od; print(and multiplying by); for j from 1 to n do if j in J then print(Z[j],^degree(P,Z[j])=,Z[j]^degree(P,Z[j])); fi; od; fi; ########## f[J] := subs({seq(Z[j]=1/Z[j] , j in J)}, P)*mul(Z[j]^degree(P,Z[j]), j in J); ########## if J<>{} then #print(The resulting polynomial is:); #print(f[J]); print(Now need to prove that f[,J,]>0 in the region:); for j from 1 to n do if j in J then print(1/xbar, <= , Z[j], < , infinity); else print(xbar, <= , Z[j],  < , infinity); fi; od; fi; print(Make new polynomial g[,J,] by substituting); for j from 1 to n do if j in J then print(Z[j]=Z[j]+1/xbar); else print(Z[j]=Z[j]+xbar); fi; od; if J={} then print(into P); else print(into f[,J,] ); fi; ########## compJ := {seq(i, i=1..n)} minus J; #{1,...,n} minus J g[J] := simplify(subs({seq(Z[j]=Z[j]+1/xbar, j in J), seq(Z[cj]=Z[cj]+xbar, cj in compJ)}, f[J])); ########## #print(The resulting polynomial is ); #print(g[J]); print(Now need to prove that g[,J,]>0 in the region:); for j from 1 to n do print(0, <= , Z[j], < , infinity); od; print(except when all variables are simultaneously zero.); print(Entering IsPos with g[,J,] ); ########## Pos := IsPos(g[J],Z); if Pos=true then ########## print(According to IsPos, g[,J,] > 0 in the positive orthant!); ########## Result[J] := true; elif Pos=false then ########## print(According to IsPos, g[,J,] is not always positive in the positive orthant.); ########## Result[J] := false; else #Pos=FAIL ########## print(The procedure IsPos was not able to tell whether g[,J,] was positive or not.); if J<>{seq(i, i=1..n)} then print(Make new polynomial h[,J,] by substituting into P ); for j from 1 to n do if not member(j,J) then print(Z[j]=1/Z[j]); fi; od; print(and then multiplying by); for j from 1 to n do if not member(j,J) then print(Z[j],^degree(P,Z[j])=,Z[j]^degree(P,Z[j])); fi; od; fi; ########## h[J] := simplify(subs({seq(Z[j]=1/Z[j], j in compJ)}, P) *mul(Z[j]^degree(P,Z[j]), j in compJ)); B[J] := []; for i from 1 to n do if i in J then B[J] := [op(B[J]), xbar]; else B[J] := [op(B[J]),1/xbar]; fi; od; ########## if J<>{seq(i, i=1..n)} then #print(The resulting polynomial is); #print(h[J]); print(Now need to prove that h[,J,]>0 in the region:); else print(Now need to prove that P>0 in the region:); fi; for j from 1 to n do if j in J then print(0, < , Z[j], <= , xbar); else print(0, < , Z[j],  <= , 1/xbar); fi; od; print(using the procedure SmallBoxMethod); ########## Result[J] := SmallBoxMethod(h[J], Z, B[J], N); ########## if Result[J]=true then if J<>{seq(i, i=1..n)} then print(SmallBoxMethod was able to prove that h[,J,] is positive in the region:); for j from 1 to n do if j in J then print(0, < , Z[j], <= , xbar); else print(0, < , Z[j],  <= , 1/xbar); fi; od; print(Therefore, P>0 in the region:); else print(SmallBoxMethod was able to prove that P is positive in the region:); fi; for j from 1 to n do if j in J then print(0, < , Z[j], <= , xbar); else print(xbar, <= , Z[j],  < , infinity); fi; od; elif Result[J]=FAIL then if J<>{seq(i, i=1..n)} then print(SmallBoxMethod was not able to conclude that h[,J,] is positive in the region); for j from 1 to n do if j in J then print(0, < , Z[j], <= , xbar); else print(0, < , Z[j],  <= , 1/xbar); fi; od; else print(SmallBoxMethod was not able to conclude that P is positive in the region:); for j from 1 to n do if j in J then print(0, < , Z[j], <= , xbar); else print(xbar, <= , Z[j],  < , infinity); fi; od; fi; else if J<>{seq(i, i=1..n)} then print(SmallBoxMethod proved that h[,J,] is negative at some points in the region); for j from 1 to n do if j in J then print(0, < , Z[j], <= , xbar); else print(0, < , Z[j],  <= , 1/xbar); fi; od; else print(SmallBoxMethod proved that P is negative at some points in the region:); for j from 1 to n do if j in J then print(0, < , Z[j], <= , xbar); else print(xbar, <= , Z[j],  < , infinity); fi; od; fi; fi; ########## fi; ########## print(============================================================); ########## od; for K in sets do print(K,Result[K]); od; RETURN({seq(Result[K], K in sets)}); end: #==================================================================== #ProveK(R,vars,K,delta,N) #Inputs: # R - a rational difference equation # vars - variables in R # K - K value to test # delta - delta value to test # N - positive integer to be used when calling PolynomialPositive #Outputs: # List of lists of length 2. First element is equilibrium point # second element is output of PolynomialPositive for that equil. #Also: # There are print statements throughout. ProveK := proc(R,vars,K,delta,N) local P,p,GAS,newvars; newvars := [seq(z[i], i=1..nops(vars))]; P := OrigPoly(R,vars,newvars,K,delta); for p in P do GAS[p[1]] := PolynomialPositive(p[2],newvars,p[1],N); od; RETURN([seq([p[1], GAS[p[1]]], p in P)]); end: #==================================================================== #Prove(R,vars,MinK,MaxK,delta,N) #Inputs: # R - A rational difference equation # vars - variables in R # MinK - the minimum K value to try # MaxK - the maximum K value to try # delta - the delta to try # N - positive integer to be used in PolynomialPositive #Outputs: # A set of lists of length 2. For each list the first element is # the equilibrium value, the second element is the K value for that equ. # If K=-1 then that equilibrium is not LAS. If K=0 then the MaxK # value is not high enough. #Method: # For each equilibrium calls PolynomialPositive with OrigPoly for # that equilibrium. Prove := proc(R,vars,MinK,MaxK,delta,N) local newvars,P,P2,LS1,p,i,flag,Ret,K,GAS,l; newvars := [seq(z[i], i=1..nops(vars))]; LS1 := LinStab(R,vars); Ret := {}; for l in LS1 do flag := false; if l[2] then for K from MinK to MaxK while not flag do print(Testing K=,K, for the equilibrium xbar=,l[1]); P := OrigPoly(R,vars,newvars,K,delta); P2 := select(i->evalb(i[1]=l[1]),P)[1][2]; #print(P2); GAS[l] := PolynomialPositiveAbsNoPrint(P2,newvars,l[1],N); print(GAS[l]); if GAS[l]={true} then flag := true; Ret := Ret union {[l[1],K]}; PolynomialPositive(P2,newvars,l[1],N); fi; od; if not flag then print(The MaxK for the equilibrium xbar=,l[1], is not large enough.); Ret := Ret union {[l[1],0]}; fi; else print(The equilibrium xbar=,l[1], is not LAS); Ret := Ret union {[l[1],-1]}; fi; od; RETURN(Ret); end: #==================================================================== #WebBook(R,vars,params,paramPoss,numParams,MinK,MaxK,delta,N) #Inputs: # R - a rational difference equation # vars - variables in R # params - parameters in R # paramPoss - possible values for parameters # numParams - number of parameter sets to test # MinK - the minimum K value to test # MaxK - the maximum K value to test # delta - the delta to test # N - a positive integer to be used when calling PolynomialPositive #Outputs: # List of lists of length 2. First element is set of parameters, second # element is output of Prove for subs(params,R). #Method: # Chooses parameters from paramPoss for which the equilibrium is # rational if possible. If "not" possible (tries 10*numParams) then # checks other parameter sets. #Also: # There are print statements throughout. WebBook := proc(R,vars,params,paramPoss,numParams,MinK,MaxK,delta,N) local Ps,P,i,j,Ret,R1,k,K,r,r2,x,xK,C; Ps := {}; Ret := {}; for i from 1 to 10*numParams while nops(Ps)FAIL then print(P); Ps := Ps union {P}; R1 := subs(P,R); C := {coeffs(numer(R1),vars),coeffs(denom(R1),vars)}; if not member(-1,{seq(sign(C[i]),i=1..nops(C))}) then Ret := Ret union {[P,Prove(R1,vars,MinK,MaxK,delta,N)]}; else Ps := Ps minus {P}; fi; else fi; od; if nops(Ps)0 when all variables are >=0 (except H(Z)>=0 when all vars=0) # false if there is an assignment, Z', of the variables in the positive # orthant s.t. H(Z')<=0 (or H(0)<0) # FAIL if the sufficient conditions checked are inconclusive #Method: # Checks two sufficient conditions: # 1. All coefficients are positive, and the constant term is positive or 0. # If the constant term is 0 then check that there is no subset of the # variables s.t. setting that subset 0 yields the 0 polynomial. # 2. If the only negative coefficient is on a term from the set # {Z[i]*Z[j], Z[i], Z[i]^2} then checks that the sub-polynomial # consisting of the terms from this set (for all i,j) is still positive # using the quadratic formula. #Also: # "Tex" means that it prints what it is doing as we go along in # such a way that the resulting text can be TeXed IsPosTex := proc(H,Z) local Os,C1,neg1,const,n,sets,flag,II,h,LargestNegDeg,tDeg,look,ii,loc, checkPos,flag2,a,b,c,disc,Cd,j,sets2,flag3,JJ,disc0,CheckPos0,d,Pos,pows; n := nops(Z); sets := powerset(n); Os := [op(sort(collect(H, Z, 'distributed'),order=tdeg(op(Z))))]; #list of terms in H C1 := [seq([coeffs(Os[i],Z),Os[i]/coeffs(Os[i],Z),i],i=1..nops(Os))]; #list of [coeff,vars,term#] neg1 := select(i->evalb(evalf(i[1])<=0),C1); #list of terms w/ negative coeff const := coeftayl(H, Z=[0$n], [0$n]); #constant term in H if nops(neg1)=0 and const>0 then #all coeffs positive and constant term positive ########## print(All coefficients are positive, and the constant term is positive,); print(so the polynomial is positive.); ########## RETURN(true); elif nops(neg1)=0 and const=0 then #all coeffs positive and constant term = 0 for II in sets minus {{seq(i,i=1..n)}} do h[II] := subs({seq(Z[i]=0, i in II)},H); if h[II]=0 then ########## print(If we set the variables $\\{, seq(Z[i],i in II),\\}$ equal 0); print(then the polynomial equals 0 which is bad.); ########## RETURN(false); fi; od; ########## print(All coefficients are positive, and the constant term is 0.); print(Also, there is no proper subset of the variables); print(for which setting them all equal 0 yields the 0 polynomial.); print(Therefore, the polynomial is 0 only when all variables are 0.); ########## RETURN(true); elif evalf(const)<0 then #constant term is negative ########## print(The constant term in the polynomial is negative.); print(Therefore, the polynomial is negative when all variables are 0.); ########## RETURN(false); else #there are negative coefficients and constant term is >= 0 ########## print(There are terms with negative coefficients, they are:); print(\\begin{center}); print( $,seq(Os[neg1[i][3]], i=1..nops(neg1)),$ ); print(\\end{center}); ########## if max({seq(degree(neg1[1][2],Z[i]),i=1..n)})>2 then ########## print(The degree of the largest negative term is greater than 2 in one of the variables,); print(so we can't use the discriminant method. We can't tell the sign of the polynomial for sure); ########## RETURN(FAIL); else ########## print(We must check that the portion of the polynomial with total degree); print(less than or equal to that of the "highest" negative coefficient term); print(is actually positive. Then, since the rest of the polynomial has no); print(negative coefficients, the whole polynomial is positive.); ########## ########## #print(The degree of the largest negative term is less than or equal to 2 in all of the); #print(variables. So we must check that the portion of the polynomial in which each variable has); #print(degree less than or equal to 2 is actually positive. Then, since the rest of); #print(the polynomial has no negative coefficients, the whole polynomial is positive.); ########## LargestNegDeg := add(degree(neg1[1][2],Z[i]),i=1..n); look := neg1[1][3]; #the index in Os of the first term w/ neg. coeff loc := look; #if look=1 then this is needed because we won't enter into the for loop below tDeg := LargestNegDeg; #will provide ending condition for the for loop for ii from 1 to look-1 while tDeg = LargestNegDeg do tDeg := add(degree(C1[look-ii][2],Z[j]), j=1..n); loc := C1[look-ii][3]; #will be the position in Os from which to start checkPos od; if LargestNegDeg=tDeg then #we got all the way to look-1 with the same total degree checkPos := add(Os[j], j=loc..nops(Os)); else #we stopped before ii=look-1 so we must add 1 to loc checkPos := add(Os[j], j=loc+1..nops(Os)); fi; #pows := Nlists0(n,2); #checkPos := add(coeftayl(H,Z=[0$n],l)*mul(Z[i]^l[i],i=1..n),l in pows); ########## print(The polynomial with negative coefficients that we want to check is positive is:); print(\\begin{center}); print($CP = , checkPos,$); print(\\end{center}); ########## flag2 := true; #will serve as "flag" in the pseudocode algorithm for j from 1 to n while flag2 do if degree(checkPos,Z[j])=2 then ########## print(The degree of$CP$in the variable$,Z[j],$is 2, so we check to see if the discriminant is negative.); ########## a := coeff(checkPos,Z[j],2); b := coeff(checkPos,Z[j],1); c := coeff(checkPos,Z[j],0); disc := b^2-4*a*c; ########## print(The discriminant is); print(\\begin{center}); print($Disc=,disc,$); print(\\end{center} ); ########## Cd := [coeffs(collect(disc,Z,'distributed'),Z)]; #list of coeffs in disc if {seq(sign(evalf(Cd[i])), i=1..nops(Cd))}={-1} then #all coeffs are negative ########## #print(According to IsPos, -disc >0); print(All coefficients in the discriminant are negative.); ########## sets2 := powerset({seq(i, i=1..j-1), seq(i,i=j+1..n)}) minus {{seq(i, i=1..j-1), seq(i,i=j+1..n)}}; flag3 := true; for JJ in sets2 while flag3 do d[JJ] := subs({seq(Z[i]=0, i in JJ)},disc); if d[JJ]=0 then ########## print(If we set the variables$\\{, seq(Z[i],i in JJ),\\}$equal 0); print(then the discriminant equals 0 which is bad.); ########## flag3 := false; fi; od; if flag3 then #there are no subsets of {1,...,j-1,j+1,...,n} s.t. setting those vars to 0 yields disc=0 ########## print(There is no proper subset of the variables which, when set to 0, yields$Disc=0$); ########## disc0 := subs({seq(Z[i]=0, i=1..j-1),seq(Z[i]=0,i=j+1..n)},disc); CheckPos0 := subs({seq(Z[i]=0, i=1..j-1),seq(Z[i]=0,i=j+1..n)},checkPos); if disc0 < 0 or (disc0 = 0 and coeff(CheckPos0,Z[j])*coeff(CheckPos0,Z[j],2)>=0) then ########## print(Setting all variables except$,Z[j],$equal 0 in$Disc$yields \$,disc0,\$ ); print(Setting all varialbes except$,Z[j],$equal 0 in$CP$yields \$,CheckPos0,\$ ) ; ########## #keep flag=true else ########## print(Setting all variables except$,Z[j],$equal 0 in$Disc$yields \$,disc0,\$ ); print(Setting all variables except$,Z[j],$equal 0 in$CP$yields \$,CheckPos0,\$ ); ########## flag2 := false; fi; else flag2 := false; fi; else #there are positive coefficients in disc ########## print(There are some positive coefficients in$Disc$, so we can't tell if the discriminant is negative.); ########## flag2 := false; fi; else #degree(checkPos,Z[j]) <> 2 ########## print(The degree of$CP$is$\\neq 2$in the variable$, Z[j],$so we can't use); print(the discriminant method.); ########## flag2 := false; fi; od; if flag2 then ########## print($CP$is positive when all variables are non-negative, and$CPmay equal 0 when all ); print(variables are 0. Therefore the polynomial is positive when all variables are ); print(non-negative and may equal 0 when all variables are 0.); ########## RETURN(true) fi; fi; RETURN(FAIL); fi; end: #==================================================================== #SmallBoxMethodTex(h,Z,B,N) #Inputs: # h - a polynomial # Z - variables in h # B - list of positive real numbers of length nops(Z) # N - an integer #Outputs: # true if h(Z)>0 for 0<=Z[i]<=B[i] for all i, except possibly # h(0,...,0)>=0 and h(B)>=0 # false if there is a n assignment, Z', of variables s.t. h(Z')<0 # FAIL if the method is inconclusive #Method: # Splits the hyperbox 0<=Z[i]<=B[i] for all i into 2^N smaller hyperboxes # (split each axis into N intervals) and uses (the correct version of) # IsPos on each of the smaller hyperboxes. #Also: # "Tex" means that it prints what it is doing as we go along in # such a way that the resulting text can be TeXed SmallBoxMethodTex := proc(h,Z,B,N) local n,interval,i,j,L,Li,Y,H,Hstar,Hstar0,Os,C1,neg1,const,tryFlips,k,flag,flag2, disc,checkPos,tDegNeg,look,tDeg,ii,loc,a,b,c,Cd,solNumer,d,Bad,Pos; ########## print(Starting the procedure SmallBoxMethod with inputs:); print(\\begin{align*}); #print(h &=, h,\\\\ ); print(Z &=, Z,\\\\ ); print(B &=, B,\\\\ ); print(N&=, N); print(\\end{align*}); ########## flag := true; n := nops(Z); for i from 1 to n do interval[i] := [seq([(j-1)*B[i]/N, j*B[i]/N], j=1..N)]; od; Li := Nlists(n,N); #All lists of length n consisting of integers between 1 and N Bad := {}; #Will be the set of all bad hyperboxes for L in Li do Y := [seq(interval[j][L[j]],j=1..n)]; ########## print(Trying to proveh>0$for$,Z,in the region:); print(\\begin{align*}); for k from 1 to n-1 do print(Y[k][1], &< , Z[k], \\leq , Y[k][2],\\\\ ); od; print(Y[n][1], &< , Z[n], \\leq , Y[n][2]); print(\\end{align*}); print(Make a new polynomial,H$, by substituting into$h:); print(\\begin{align*}); for k from 1 to n-1 do print(Z[k],&=,Z[k]+Y[k][1],\\\\ ); od; print(Z[n],&=,Z[n]+Y[n][1]); print(\\end{align*}); ########## H := simplify(subs({seq(Z[ii]=Z[ii]+Y[ii][1], ii=1..n)} , h)); ########## print(Make a new polynomial,H^\\star$, by substituting into$H:); print(\\begin{align*}); for k from 1 to n-1 do print(Z[k],&=,1/Z[k],\\\\ ); od; print(Z[n],&=,1/Z[n]); print(\\end{align*}); print(and then multiplying by); print(\\begin{align*}); for k from 1 to n-1 do print(Z[k],^{degree(H, , Z[k],)}&=,Z[k]^degree(H,Z[k]),\\\\ ); od; print(Z[n],^{degree(H, , Z[n],)}&=,Z[n]^degree(H,Z[n])); print(\\end{align*}); ########## Hstar := subs({seq(Z[ii]=1/Z[ii], ii=1..n)}, H)*mul(Z[jj]^degree(H,Z[jj]), jj=1..n); ########## print(Make a new polynomial,H^{\\star}_0$, by substituting into$H^{\\star}:); print(\\begin{align*}); for k from 1 to n-1 do print(Z[k],&=,Z[k]+1/(Y[k][2]-Y[k][1]),\\\\ ); od; print(Z[n],&=,Z[n]+1/(Y[n][2]-Y[n][1])); print(\\end{align*}); ########## Hstar0 := simplify(subs({seq(Z[ii]=Z[ii]+1/(Y[ii][2]-Y[ii][1]), ii=1..n)}, Hstar)); #Hstar0 := subs({seq(Z[ii]=Z[ii]+1/(Y[ii][2]-Y[ii][1]), ii=1..n)}, Hstar); ########## print(Need to prove thatH^{\\star}_0(,op(Z),) > 0$for all variables$\\geq 0$using IsPos.); ########## const := coeftayl(Hstar0, Z=[0$n], [0$n]); #the constant term in Hstar0 Pos := IsPosTex(Hstar0,Z); if Pos=true then if const=0 then if not [seq(Y[i][2],i=1..n)]=B then if not [seq(Y[i][1],i=1..n)]=[0$n] then ########## print(IsPos returned true in this hyperbox, but the constant term); print(is 0, and the box is not the furthest from the origin or closest to the origin.); print(Therefore, there is a point where $h=0$ that we don't want $h=0$.); #print(Y,false); ########## Bad := Bad union {Y}; flag := false; RETURN(false); else ########## print(The polynomial, $h$, is positive on the hyperbox ); print(\\begin{align*}); for k from 1 to n-1 do print(Y[k][1], &< , Z[k], \\leq , Y[k][2],\\\\ ); od; print(Y[n][1], &< , Z[n], \\leq , Y[n][2]); print(\\end{align*}); print(except at the origin:); #print(Y,true); ########## #keep flag what it is #continue to examine the hyperboxes fi; else ########## print(The polynomial, $h$, is positive on the hyperbox:); print(\\begin{align*}); for k from 1 to n-1 do print(Y[k][1], &< , Z[k], \\leq , Y[k][2],\\\\ ); od; print(Y[n][1], &< , Z[n], \\leq , Y[n][2]); print(\\end{align*}); #print(Y,true); ########## #keep flag what it is #continue to examine the hyperboxes fi; else ########## print(The polynomial, $h$, is positive on the hyperbox:); print(\\begin{align*}); for k from 1 to n-1 do print(Y[k][1], &< , Z[k], \\leq , Y[k][2],\\\\ ); od; print(Y[n][1], &< , Z[n], \\leq , Y[n][2]); print(\\end{align*}); #print(Y,true); ########## #keep flag what it is #continue to examine the hyperboxes fi; elif Pos=false then ########## print(IsPos returned false, so $h$ is not positive on the hyperbox:); print(\\begin{align*}); for k from 1 to n-1 do print(Y[k][1], &< , Z[k], \\leq , Y[k][2],\\\\ ); od; print(Y[n][1], &< , Z[n], \\leq , Y[n][2]); print(\\end{align*}); #print(Y,false); ########## Bad := Bad union {Y}; flag := false; RETURN(false); else #Pos=FAIL then ########## print(IsPos was inconclusive, so the sign of $h$ is not known on the hyperbox:); print(\\begin{align*}); for k from 1 to n-1 do print(Y[k][1], &< , Z[k], \\leq , Y[k][2],\\\\ ); od; print(Y[n][1], &< , Z[n], \\leq , Y[n][2]); print(\\end{align*}); #print(Y,FAIL); ########## #don't add Y to Bad since we don't know if Hstar0 is < 0 on Y flag := false; fi; print(\\\\-------------------------------------------------------------\\\\); od; if flag then RETURN(true); else if Bad={} then #flag was set to false but only because IsPos couldn't #determine the sign of Hstar0 for some Y RETURN(FAIL); else ########## print(The hyperboxes for which $h$ is negative is:); print( \$,Bad,\$ ); ########## #flag was set to false becuase Hstar0 is <0 for some point #in each of the hyperboxes in Bad RETURN(false); fi; fi; end: #==================================================================== #PolynomialPositiveTex(P,Z,xbar,N)); #Inputs:); # P - a polynomial # Z - variables in P # xbar - a positive real number which marks the first cut of Rplus^|Z| into 2^|Z| regions # N - an integer greater than 1. Will be used when calling SmallBoxMethod if necessary #Outputs: # Set containing elements from {true, false, FAIL}. # if {true} then # P(Z)>0 when Z[i]>=0 except possibly P(0,...,0)>=0 and P(xbar,...,xbar)>=0 # elif false in output then # there is a point where P(Z')<=0 when Z'[i]>=0 (or P(0,...,0)<0 or P(xbar,...,xbar)<0 ) # else (FAIL in output and false not in output) # there is some place where using IsPos is inconclusive #Method: # Cuts the positive orthant into 2^|Z| subsets and checks on each of # these whether P(Z)>0. Each subset, I, of {1,...,|Z|} corresponds to the # subset of the positive orthant where 0<=Z[i]<=xbar if i in I, # and xbar<=Z[j] if j not in I. #Also: # "Tex" means that it prints what it is doing as we go along in # such a way that the resulting text can be TeXed PolynomialPositiveTex := proc(P,Z,xbar,N) local n,sets,J,f,g,C,neg,const,Result,h,compJ,B,i,K,j,g2,g3,flag,II, Os,C1,neg1,tDegNeg,look,tDeg,loc,checkPos,ii,disc,a,b,c,d,flag2,Cd,solNumer,Pos; ##### print(Starting the procedure PolynomialPositive with inputs:); print(\\begin{align*}); print(P&=\\text{ from above}\\\\ ); print(Z&=,Z,\\\\ ); print(\\xbar&=,xbar,\\\\ ); print(N&=,N); print(\\end{align*}); ##### n := nops(Z); sets := powerset(n); if xbar=0 then ########## print(Since $\\xbar=0$ we check that all coefficients); print(and the constant term in P are positive); ########## C := [coeffs(collect(P, Z, 'distributed'), Z)]; neg := select(i->evalb(i<=0) ,C); const := coeftayl(P, Z=[0$n], [0$n]); if nops(neg)=0 and const>0 then ########## print(All coefficients in $P$ and the constant term are positive,); print(so $P>0$ for all variables $\\geq 0$); ########## RETURN({true}); elif nops(neg)=0 and const=0 then ########## print(All coefficients in $P$ are positive, and the constant term is 0); print(So $P\\geq 0$ when all variables are $\\geq 0$); ########## RETURN({true}); else ########## if nops(neg)>0 then print(The negative coefficients are $,neg,$ ); fi; if const<=0 then print(The constant term is $,const,$ ); fi; print(so this method doesn't prove that $P>0$); ########## RETURN({false}); fi; fi; for J in sets do #print($\\{, op(J),\\}$); #print(\\\\ ); ########## print(Proving $P>0$ in the region:); print(\\begin{align*}); for j from 1 to n-1 do if j in J then print(0, <& , Z[j], \\leq , xbar,\\\\ ); else print(xbar, \\leq & , Z[j],  < \\infty \\\\); fi; od; if n in J then print(0, <& , Z[n], \\leq , xbar); else print(xbar, \\leq & , Z[n],  < \\infty); fi; print(\\end{align*}); if J<>{} then print(Make new polynomial $f[\\{,op(J),\\}]$ by substituting ); print(\\begin{align*}); for j from 1 to n-1 do if j in J then print(Z[j],&=,1/Z[j],\\\\ ); fi; od; if n in J then print(Z[n],&=,1/Z[n]); fi; print(\\end{align*}); print(and multiplying by); print(\\begin{align*}); for j from 1 to n-1 do if j in J then print(Z[j],^{degree(P, , Z[j],)}&=,Z[j]^degree(P,Z[j]),\\\\ ); fi; od; if n in J then print(Z[n],^{degree(P, , Z[n],)}&=,Z[n]^degree(P,Z[n])); fi; print(\\end{align*}); fi; ########## f[J] := subs({seq(Z[j]=1/Z[j] , j in J)}, P)*mul(Z[j]^degree(P,Z[j]), j in J); ########## if J<>{} then #print(The resulting polynomial is:); #print(f[J]); print(Now need to prove that $f[\\{,op(J),\\}]>0$ in the region:); print(\\begin{align*}); for j from 1 to n-1 do if j in J then print(1/xbar, \\leq & , Z[j], < \\infty \\\\); else print(xbar, \\leq & , Z[j],  < \\infty \\\\); fi; od; if n in J then print(1/xbar, \\leq & , Z[n], < \\infty); else print(xbar, \\leq & , Z[n],  < \\infty); fi; print(\\end{align*}); fi; print(Make new polynomial $g[\\{,op(J),\\}]$ by substituting); print(\\begin{align*} ); for j from 1 to n-1 do if j in J then print(Z[j],=& ,Z[j]+1/xbar,\\\\ ); else print(Z[j],=& ,Z[j]+xbar,\\\\ ); fi; od; if n in J then print(Z[n],=& ,Z[n]+1/xbar); else print(Z[n],=& ,Z[n]+xbar); fi; print(\\end{align*}); if J={} then print(into $P$); else print(into $f[\\{,op(J),\\}]$ ); fi; ########## compJ := {seq(i, i=1..n)} minus J; #{1,...,n} minus J g[J] := simplify(subs({seq(Z[j]=Z[j]+1/xbar, j in J), seq(Z[cj]=Z[cj]+xbar, cj in compJ)}, f[J])); ########## #print(The resulting polynomial is ); #print(g[J]); print(Now need to prove that $g[\\{,op(J),\\}]>0$ in the region:); print(\\begin{align*}); for j from 1 to n-1 do print(0, \\leq & , Z[j], < \\infty \\\\); od; print(0, \\leq & , Z[n], < \\infty \\\\); print(\\end{align*}); print(except when all variables are simultaneously zero.); print(\\\\ ); print(Entering IsPos with $g[\\{,op(J),\\}]$ ); print(\\\\ ); ########## Pos := IsPosTex(g[J],Z); if Pos=true then ########## print(\\\\ ); print(According to IsPos, $g[\\{,op(J),\\}] > 0$ in the positive orthant!); ########## Result[J] := true; elif Pos=false then ########## print(\\\\ ); print(According to IsPos, $g[\\{,op(J),\\}]$ is not always positive in the positive orthant.); ########## Result[J] := false; else #Pos=FAIL ########## print(\\\\ ); print(The procedure IsPos was not able to tell whether $g[\\{,op(J),\\}]$ was positive or not.); if J<>{seq(i, i=1..n)} then print(Make new polynomial $h[\\{,op(J),\\}]$ by substituting into $P$ ); print(\\begin{align*}); for j from 1 to max(J)-1 do if not member(j,J) then print(Z[j],=& ,1/Z[j],\\\\ ); fi; od; print(Z[max(J)],=& ,1/Z[max(J)]); print(\\end{align*}); print(and then multiplying by); print(\\begin{align*}); for j from 1 to n-1 do if not member(j,J) then print(Z[j],^{degree(P, , Z[j],)}&=,Z[j]^degree(P,Z[j]),\\\\ ); fi; od; if not member(n,J) then print(Z[n],^{degree(P, , Z[n],)}&=,Z[n]^degree(P,Z[n])); fi; print(\\end{align*}); fi; ########## h[J] := simplify(subs({seq(Z[j]=1/Z[j], j in compJ)}, P) *mul(Z[j]^degree(P,Z[j]), j in compJ)); B[J] := []; for i from 1 to n do if i in J then B[J] := [op(B[J]), xbar]; else B[J] := [op(B[J]),1/xbar]; fi; od; ########## if J<>{seq(i, i=1..n)} then #print(The resulting polynomial is); #print(h[J]); print(Now need to prove that $h[\\{,op(J),\\}]>0$ in the region:); else print(Now need to prove that $P>0$ in the region:); fi; print(\\begin{align*}); for j from 1 to n-1 do if j in J then print(0, < & , Z[j], \\leq , xbar, \\\\); else print(0, < &, Z[j],  \\leq , 1/xbar, \\\\); fi; od; if n in J then print(0, < & , Z[n], \\leq , xbar); else print(0, < &, Z[n],  \\leq , 1/xbar); fi; print(\\end{align*}); print(using the procedure SmallBoxMethod); ########## Result[J] := SmallBoxMethodTex(h[J], Z, B[J], N); ########## if Result[J]=true then if J<>{seq(i, i=1..n)} then print(SmallBoxMethod was able to prove that $h[\\{,op(J),\\}]$ is positive in the region:); print(\\begin{align*}); for j from 1 to n-1 do if j in J then print(0, < & , Z[j], \\leq , xbar,\\\\ ); else print(0, < & , Z[j],  \\leq , 1/xbar,\\\\ ); fi; od; if n in J then print(0, < & , Z[n], \\leq , xbar); else print(0, < & , Z[n],  \\leq , 1/xbar); fi; print(\\end{align*}); print(Therefore, $P>0$ in the region:); else print(SmallBoxMethod was able to prove that $P$ is positive in the region:); fi; print(\\begin{align*}); for j from 1 to n-1 do if j in J then print(0, < & , Z[j], \\leq , xbar,\\\\ ); else print(xbar, \\leq & , Z[j],  < \\infty \\\\); fi; od; if n in J then print(0, < & , Z[n], \\leq , xbar); else print(xbar, \\leq & , Z[n],  < \\infty); fi; print(\\end{align*}); elif Result[J]=FAIL then if J<>{seq(i, i=1..n)} then print(SmallBoxMethod was not able to conclude that $h[\\{,op(J),\\}]$ is positive in the region); print(\\begin{align*}); for j from 1 to n-1 do if j in J then print(0, < & , Z[j], \\leq , xbar,\\\\ ); else print(0, < & , Z[j],  \\leq , 1/xbar,\\\\ ); fi; od; if n in J then print(0, < & , Z[n], \\leq , xbar); else print(0, < & , Z[n],  \\leq , 1/xbar); fi; print(\\begin{align*}); else print(SmallBoxMethod was not able to conclude that $P$ is positive in the region:); print(\\begin{align*}); for j from 1 to n-1 do if j in J then print(0, < &, Z[j], \\leq , xbar,\\\\ ); else print(xbar, \\leq & , Z[j],  < \\infty \\\\); fi; od; if n in J then print(0, < &, Z[n], \\leq , xba); else print(xbar, \\leq & , Z[n],  < \\infty); fi; print(\\end{align*}); fi; else if J<>{seq(i, i=1..n)} then print(SmallBoxMethod proved that $h[\\{,op(J),\\}]$ is negative at some points in the region); print(\\begin{align*}); for j from 1 to n-1 do if j in J then print(0, < & , Z[j], \\leq , xbar,\\\\ ); else print(0, < & , Z[j],  \\leq , 1/xbar,\\\\ ); fi; od; if n in J then print(0, < & , Z[n], \\leq , xbar); else print(0, < & , Z[n],  \\leq , 1/xbar); fi; print(\\end{align*}); else print(SmallBoxMethod proved that $P$ is negative at some points in the region:); print(\\begin{align*}); for j from 1 to n-1 do if j in J then print(0, < & , Z[j], \\leq , xbar,\\\\ ); else print(xbar, \\leq & , Z[j],  < \\infty \\\\); fi; od; if n in J then print(0, < & , Z[n], \\leq , xbar); else print(xbar, \\leq & , Z[n],  < \\infty); fi; print(\\end{align*}); fi; fi; ########## fi; ########## print(\\\\-------------------------------------------------------------\\\\); ########## od; #for K in sets do # print(K,Result[K]); #od; RETURN({seq(Result[K], K in sets)}); end: #==================================================================== #ProveTex(R,vars,MinK,MaxK,delta,N) #Inputs: # R - A rational difference equation # vars - variables in R # MinK - the minimum K value to try # MaxK - the maximum K value to try # delta - the delta to try # N - positive integer to be used in PolynomialPositive #Outputs: # A set of lists of length 2. For each list the first element is # the equilibrium value, the second element is the K value for that equ. # If K=-1 then that equilibrium is not LAS. If K=0 then the MaxK # value is not high enough. #Method: # For each equilibrium calls PolynomialPositive with OrigPoly for # that equilibrium. #Also: # "Tex" means that it prints what it is doing as we go along in # such a way that the resulting text can be TeXed ProveTex := proc(R,vars,MinK,MaxK,delta,N) local newvars,P,P2,LS1,p,i,flag,Ret,K,GAS,l,Zmap; newvars := [seq(z[i], i=1..nops(vars))]; LS1 := LinStab(R,vars); Ret := {}; for l in LS1 do flag := false; print(First we check that the equilibrium ,l[1], is LAS.\\\\); if l[2] then print(It is LAS, so we continue to test $K$ values.); for K from MinK to MaxK while not flag do ########## print(\\\\ ); print(Testing $K=,K,$ for the equilibrium $\\xbar=,l[1],$ \\\\); ########## P := OrigPoly(R,vars,newvars,K,delta); P2 := select(i->evalb(i[1]=l[1]),P)[1][2]; #print(P2); GAS[l] := PolynomialPositiveAbsNoPrint(P2,newvars,l[1],N); ########## print(For $K=,K,$ we get $\\{,op(GAS[l]),\\}$ output from PolynomialPositive. ); ########## if GAS[l]={true} then flag := true; Ret := Ret union {[l[1],K]}; ########## print(\\begin{thm}); print(The equilibrium ,l[1], for the rational difference equation \$x[n+1] = ,R,\$ is GAS.); print(\\end{thm}); print(\\begin{proof}); print(From the rational difference equation, the $K$ value $,K,$, and the $\\delta$ value, delta, we get a polynomial:\\\\); print( $P = ,P2,$ \\\\); print(The goal is to prove that this polynomial is positive when all variables are positive. We create this polynomial by:); print(\\begin{enumerate}); print(\\item Create a new rational difference equation, by setting $z[n]=x[n]-,l[1],$ in the original difference equation. This new difference equation has 0 as its equilibrium, and we now wish to prove that the equilibrium 0 is GAS when initial conditions are greater than $,-l[1],$. The new difference equation is:); Zmap := normal(subs({seq(vars[i]=z[op(vars[i])]+l[1],i=1..nops(vars))},R)-l[1]); print(\$z[n+1]=, Zmap,\$ ); print(\\item Consider the vector valued mapping, $Q:\\R^{, nops(vars),}\\rightarrow \\R^{, nops(vars),}$, which represents the new difference equation: ); if nops(vars)=1 then print(\$Q(\\tup{z[n]}) = \\tup{,Zmap,}\$); elif nops(vars)=2 then print(\$Q(\\tup{z[n-1], z[n]}) = \\tup{z[n], ,Zmap,}\$); elif nops(vars)=3 then print(\$Q(\\tup{z[n-2],\\ldots,z[n]}) = \\tup{z[n-1], z[n], ,Zmap,}\$); else print(\$Q(\\tup{z[n-,nops(vars)-1,],\\ldots,z[n]}) = \\tup{z[n-,nops(vars)-2,],\\ldots,z[n], ,Zmap,}\$); fi; print(\\item Notice that if); print(\$\\frac{\\left| Q^K(\\tup{,op(newvars),})\\right|^2}{\\left|\\tup{,op(newvars),}\\right|^2}<,delta,\$ ); print(for all $,op(newvars), > ,-l[1],$ then we are done (see applicable Theorem in Emilie Hogan's PhD thesis). Simplify this so that it is of the form something $>0$:); print(\$0< ,delta,*\\left|\\tup{,op(newvars),}\\right|^2 - \\left| Q^K(\\tup{,op(newvars),}) \\right|^2 \$); print(and then take the numerator. The denominator will always be a product of squares and sums of squares, so if the numerator is positive when all variables are greater than ,-l[1], then the whole expression is positive.); print(\\item Finally, replace $z[i]$ by $z[i]-,l[1],$ so that we have a polynomial that we wish to prove is positive when all variables are positive (rather than when all variables are greater than some negative number).); print(\\end{enumerate}); print(Now we run the algorithm PolynomialPositive on this polynomial. If the polynomial is positive when all variables are positive then the equilibrium ,l[1], is GAS for the original difference equation.); ########## PolynomialPositiveTex(P2,newvars,l[1],N); ########## print(Since $P>0$ when all variables are positive, the $K$ value , K, is proven to work for the rational difference equation \$x[n+1]=,R,\$ ); print(\\end{proof}); print(========================================\\\\); ########## fi; od; if not flag then print(The $MaxK$ for the equilibrium $\\xbar=,l[1],$ is not large enough.); Ret := Ret union {[l[1],0]}; print(\\\\========================================\\\\); fi; else print(The equilibrium $\\xbar=,l[1],$ is not LAS); Ret := Ret union {[l[1],-1]}; print(\\\\========================================\\\\); fi; od; RETURN(Ret); end: #==================================================================== #WebBookTex(R,vars,params,paramPoss,numParams,MinK,MaxK,delta,N) #Inputs: # R - a rational difference equation # vars - variables in R # params - parameters in R # paramPoss - possible values for parameters # numParams - number of parameter sets to test # MinK - the minimum K value to test # MaxK - the maximum K value to test # delta - the delta to test # N - a positive integer to be used when calling PolynomialPositive #Outputs: # List of lists of length 2. First element is set of parameters, second # element is output of Prove for subs(params,R). #Method: # Chooses parameters from paramPoss for which the equilibrium is # rational if possible. If "not" possible (tries 10*numParams) then # checks other parameter sets. #Also: # "Tex" means that it prints what it is doing as we go along in # such a way that the resulting text can be TeXed WebBookTex := proc(R,vars,params,paramPoss,numParams,MinK,MaxK,delta,N) local Ps,P,i,j,Ret,R1,k,K,r,r2,x,xK,C,jj; Ps := {}; Ret := {}; ########## print(For the rational difference equation); print(\$x[n+1]=,R,\$); print(We will try to prove that the equilibrium is GAS for various values of the parameters $\\{,op(params),\\}$. ); ########## for i from 1 to 10*numParams while nops(Ps)FAIL then ########## print(\\\\For the parameters $\\{,op(P),\\}$:); ########## Ps := Ps union {P}; R1 := subs(P,R); C := {coeffs(numer(R1),vars),coeffs(denom(R1),vars)}; if not member(-1,{seq(sign(C[i]),i=1..nops(C))}) then Ret := Ret union {[P,ProveTex(R1,vars,MinK,MaxK,delta,N)]}; else Ps := Ps minus {P}; fi; else fi; od; if nops(Ps)0 then print(The parameter values for which the equilibrium is not LAS are:); print(\\begin{center}); print($); ######### #for jj from 1 to nops(K[-1])-1 do # print([\\{,seq(params[i]=K[-1][jj][1][i],i=1..nops(params)),\\}, ,K[-1][jj][2] ,], ); #od; #print([\\{,seq(params[i]=K[-1][nops(K[-1])][1][i],i=1..nops(params)),\\}, ,K[-1][nops(K[-1])][2] ,] ); ######### for jj from 1 to nops(K[-1])-1 do print([\\{,op(K[-1][jj][1]),\\}, ,K[-1][jj][2] ,], ); od; print([\\{,op(K[-1][nops(K[-1])][1]),\\}, ,K[-1][nops(K[-1])][2] ,], ); print($ ); print(\\end{center}); fi; if nops(K[0])>0 then print(The parameter values for which the MaxK value, ,MaxK, , is not high enough are:); print(\\begin{center}); print($); ######### #for jj from 1 to nops(K[0])-1 do # print([\\{,seq(params[i]=K[0][jj][1][i],i=1..nops(params)),\\}, ,K[0][jj][2] ,], ); #od; #print([\\{,seq(params[i]=K[0][nops(K[0])][1][i],i=1..nops(params)),\\}, ,K[0][nops(K[0])][2],] ); ######### for jj from 1 to nops(K[0])-1 do print([\\{,op(K[0][jj][1]),\\}, ,K[0][jj][2] ,], ); od; print([\\{,op(K[0][nops(K[0])][1]),\\}, ,K[0][nops(K[0])][2] ,], ); print($ ); print(\\end{center}); fi; for k from MinK to MaxK do if nops(K[k])>0 then print(The parameter values for which the $K=,k,$ are:); print(\\begin{center}); print( $); ######### #for jj from 1 to nops(K[k])-1 do # print([\\{,seq(params[i]=K[k][jj][1][i],i=1..nops(params)),\\}, ,K[k][jj][2] ,], ); #od; #print([\\{,seq(params[i]=K[k][nops(K[k])][1][i],i=1..nops(params)),\\}, ,K[k][nops(K[k])][2],] ); ######### for jj from 1 to nops(K[k])-1 do print([\\{,op(K[k][jj][1]),\\}, ,K[k][jj][2] ,], ); od; print([\\{,op(K[k][nops(K[k])][1]),\\}, ,K[k][nops(K[k])][2] ,], ); print($ ); print(\\end{center} ); fi; od; RETURN(Ret); end: #==================================================================== #ProveAlmostNoPrintTex(R,vars,MinK,MaxK,delta,N) #Inputs: # R - A rational difference equation # vars - variables in R # MinK - the minimum K value to try # MaxK - the maximum K value to try # delta - the delta to try # N - positive integer to be used in PolynomialPositive #Outputs: # A set of lists of length 2. For each list the first element is # the equilibrium value, the second element is the K value for that equ. # If K=-1 then that equilibrium is not LAS. If K=0 then the MaxK # value is not high enough. #Method: # For each equilibrium calls PolynomialPositive with OrigPoly for # that equilibrium. #Also: # "AlmostNoPrintTex" means that it prints some of what it is doing # as we go along in such a way that the resulting text can be TeXed ProveAlmostNoPrintTex := proc(R,vars,MinK,MaxK,delta,N) local newvars,P,P2,LS1,p,i,flag,Ret,K,GAS,l; newvars := [seq(z[i], i=1..nops(vars))]; LS1 := LinStab(R,vars); Ret := {}; for l in LS1 do flag := false; if l[2] then for K from MinK to MaxK while not flag do ########## print(\\\\ ); print(Testing $K=,K,$ for the equilibrium $\\xbar=,l[1],$ \\\\); ########## P := OrigPoly(R,vars,newvars,K,delta); P2 := select(i->evalb(i[1]=l[1]),P)[1][2]; #print(P2); GAS[l] := PolynomialPositiveAbsNoPrint(P2,newvars,l[1],N); ########## print(For $K=,K,$ we get $\\{,op(GAS[l]),\\}$ output from PolynomialPositive. ); ########## if GAS[l]={true} then flag := true; Ret := Ret union {[l[1],K]}; print(For the equilibrium $\\xbar=, l[1],$ the $K$ value , K, is proven to work.); print(\\\\========================================\\\\); fi; od; if not flag then print(The $MaxK$ for the equilibrium $\\xbar=,l[1],$ is not large enough.); Ret := Ret union {[l[1],0]}; print(\\\\========================================\\\\); fi; else print(The equilibrium $\\xbar=,l[1],$ is not LAS); Ret := Ret union {[l[1],-1]}; print(\\\\========================================\\\\); fi; od; RETURN(Ret); end: #==================================================================== #WebBookAlmostNoPrintTex(R,vars,params,paramPoss,numParams,MinK,MaxK,delta,N) #Inputs: # R - a rational difference equation # vars - variables in R # params - parameters in R # paramPoss - possible values for parameters # numParams - number of parameter sets to test # MinK - the minimum K value to test # MaxK - the maximum K value to test # delta - the delta to test # N - a positive integer to be used when calling PolynomialPositive #Outputs: # List of lists of length 2. First element is set of parameters, second # element is output of Prove for subs(params,R). #Method: # Chooses parameters from paramPoss for which the equilibrium is # rational if possible. If "not" possible (tries 10*numParams) then # checks other parameter sets. #Also: # "AlmostNoPrintTex" means that it prints some of what it is doing as we go along in # such a way that the resulting text can be TeXed WebBookAlmostNoPrintTex := proc(R,vars,params,paramPoss,numParams,MinK,MaxK,delta,N) local Ps,P,i,j,Ret,R1,k,K,r,r2,x,xK,C,jj; Ps := {}; Ret := {}; ########## print(For the rational difference equation); print(\$x[n+1]=,R,\$); print(We will try to prove that the equilibrium is GAS for various ); print(values of the parameters $\\{,op(params),\\}$. ); ########## for i from 1 to 10*numParams while nops(Ps)FAIL then ########## print(\\\\For the parameters $\\{,op(P),\\}$:); ########## Ps := Ps union {P}; R1 := subs(P,R); C := {coeffs(numer(R1),vars),coeffs(denom(R1),vars)}; if not member(-1,{seq(sign(C[i]),i=1..nops(C))}) then Ret := Ret union {[P,ProveAlmostNoPrintTex(R1,vars,MinK,MaxK,delta,N)]}; else Ps := Ps minus {P}; fi; else fi; od; if nops(Ps)0 then print(The parameter values and equilibrium for which the equilibrium is not LAS are:); print(\\begin{center}); print( $); ######### #for jj from 1 to nops(K[-1])-1 do # print([\\{,seq(params[i]=K[-1][jj][1][i],i=1..nops(params)),\\}, ,K[-1][jj][2] ,], ); #od; #print([\\{,seq(params[i]=K[-1][nops(K[-1])][1][i],i=1..nops(params)),\\}, ,K[-1][nops(K[-1])][2] ,] ); ######### for jj from 1 to nops(K[-1])-1 do print([\\{,op(K[-1][jj][1]),\\}, ,K[-1][jj][2] ,], ); od; print([\\{,op(K[-1][nops(K[-1])][1]),\\}, ,K[-1][nops(K[-1])][2] ,], ); print($ ); print(\\end{center}); fi; if nops(K[0])>0 then print(The parameter values and equilibrium for which the MaxK value, ,MaxK, , is not high enough are:); print(\\begin{center}); print( $); ######### #for jj from 1 to nops(K[0])-1 do # print([\\{,seq(params[i]=K[0][jj][1][i],i=1..nops(params)),\\}, ,K[0][jj][2] ,], ); #od; #print([\\{,seq(params[i]=K[0][nops(K[0])][1][i],i=1..nops(params)),\\}, ,K[0][nops(K[0])][2],] ); ######### for jj from 1 to nops(K[0])-1 do print([\\{,op(K[0][jj][1]),\\}, ,K[0][jj][2] ,], ); od; print([\\{,op(K[0][nops(K[0])][1]),\\}, ,K[0][nops(K[0])][2] ,], ); print($ ); print(\\end{center}); fi; for k from MinK to MaxK do if nops(K[k])>0 then print(The parameter values and equilibrium for which the $K=,k,$ are:); print(\begin{center}); print( $); ######### #for jj from 1 to nops(K[k])-1 do # print([\\{,seq(params[i]=K[k][jj][1][i],i=1..nops(params)),\\}, ,K[k][jj][2] ,], ); #od; #print([\\{,seq(params[i]=K[k][nops(K[k])][1][i],i=1..nops(params)),\\}, ,K[k][nops(K[k])][2],] ); ######### for jj from 1 to nops(K[k])-1 do print([\\{,op(K[k][jj][1]),\\}, ,K[k][jj][2] ,], ); od; print([\\{,op(K[k][nops(K[k])][1]),\\}, ,K[k][nops(K[k])][2] ,], ); print($ ); print(\end{center} ); fi; od; RETURN(Ret); end: #========================================================== #GoodParams(R,vars,params,possParams) #Inputs: # R - a rational difference equation # vars - variables in R # params - parameters in R # possParams - possible values for elements of params #Outputs: # A random parameter assignment such that the positive # equilibrium of R with this parameter assignment is rational. GoodParams := proc(R,vars,params,possParams,n) local ra,flag,P,E,i; ra := rand(1..nops(possParams)); flag := false; for i from 1 to n while not flag do P := {seq(params[i]=possParams[ra()],i=1..nops(params))}; #print(GP,P); if not type(subs(P,R),constant) then E := FindPosEquil(subs(P,R),vars); if {seq(type(convert(E[i],fraction),rational),i=1..nops(E))}={true} then flag:=true; fi; fi; od; if not flag then RETURN(FAIL); else RETURN(P); fi; end: #---------------------------------------------------------- #RandParams(R,vars,params,possParams) #Inputs: # R - a rational difference equation # vars - variables in R # params - parameters in R # possParams - possible values for elements of params #Outputs: # A random parameter assignment from the set possParams RandParams := proc(R,vars,params,possParams) local ra,flag,P,E; ra := rand(1..nops(possParams)); flag := false; while not flag do P := {seq(params[i]=possParams[ra()],i=1..nops(params))}; #print(RP,P); if not type(subs(P,R),constant) then flag:=true; fi; od; RETURN(P); end: #==================================================================== #Nlists(n,N) #Inputs: # n - a positive integer # N - a positive integer #Outputs: # Lists of length n consisting of integers between 1 and N (inclusively) Nlists := proc(n,N) local pre,Ret,p; if n=0 then RETURN({}); fi; if n=1 then RETURN({seq([i],i=1..N)}); fi; pre := Nlists(n-1,N); Ret := {}; for p in pre do Ret := Ret union {seq([op(p),j],j=1..N)}; od; RETURN(Ret); end: #==================================================================== #Nlists0(n,N) #Inputs: # n - a positive integer # N - a positive integer #Outputs: # Lists of length n consisting of integers between 0 and N (inclusively) Nlists0 := proc(n,N) local pre,Ret,p; if n=0 then RETURN({}); fi; if n=1 then RETURN({seq([i],i=0..N)}); fi; pre := Nlists0(n-1,N); Ret := {}; for p in pre do Ret := Ret union {seq([op(p),j],j=0..N)}; od; RETURN(Ret); end: #==================================================================== #OrigPoly(R,vars,newvars,k,delta) #Inputs: # -R A rational difference equation with equilibrium xBar # -vars List of variables in R in order [x[n-k],...,x[n]] # -newvars List of variables to replace vars when the equilibrium is moved to 0 # -k The k in |Q^k(newvars)|^2/|newvars|^2 # -delta We hope to prove |Q^k(newvars)|^2/|newvars|^2 <= delta<1 #Outputs: # The polynomial P=numer(delta)*|newvars|^2-|Q^k(newvars)|^2*denom(delta) # shifted to the right in each variable by xBar. We want to prove that # P>=0 for all newvars>=0 OrigPoly := proc(R,vars,newvars,k,delta) option remember; local F1,Ret,f1,xBar,f; F1 := ToMin(R,vars,newvars,k,delta); Ret := []; for f1 in F1 do xBar := f1[1]; f := subs({seq(newvars[i]=newvars[i]-xBar,i=1..nops(newvars))},f1[2]); Ret := [op(Ret),[xBar,f]]; od; RETURN(Ret); end: #========================================================== #ToMin(R1,vars,newvars,k,delta) #Inputs: # -R A rational difference equation with equilibrium xBar # -vars List of variables in R in order [x[n-k],...,x[n]] # -newvars List of variables to replace vars when the equilibrium is moved to 0 # -k The k in |Q^k(newvars)|^2/|newvars|^2 # -delta We hope to prove |Q^k(newvars)|^2/|newvars|^2 <= delta<1 #Outputs: # The polynomial P=numer(delta)*|newvars|^2-|Q^k(newvars)|^2*denom(delta). We want # to prove that P>=0 for all newvars>=-xbar. Coefficients must be numerical # since FindPosEquil is called (within MakeEquil0). ToMin := proc(R1,vars,newvars,k,delta) local iterK,n,R,Ret,r,Rat; R := MakeEquil0(R1,vars,newvars); Ret := {}; for r in R do iterK := Iter(r[2],newvars,k); #print(iterK); n := nops(vars); Rat := add(iterK[i]^2,i=1..nops(iterK))/add(newvars[i]^2,i=1..n); Ret := Ret union {[r[1], -numer(Rat)*denom(delta)+denom(Rat)*numer(delta)]}; od; RETURN(Ret); end: #========================================================== #Iter(R,vars,k): #Inputs: # R - a rational difference equation # vars - variables in R # k - number of times to iterate R #Outputs: # Iterates the map # T:(vars)->(vars[2..nops(vars)],R(vars)) # k times. The output is the vector, T^k(vars), of length nops(vars). Iter:=proc(R,vars,k) local gu,i: gu:=vars: for i from 1 to k do gu:=[op(gu[2..nops(gu)]),normal(subs({seq(vars[i]=gu[i],i=1..nops(vars))},R))]: od: RETURN(gu): end: #========================================================== #MakeEquil0(R,vars,newvars): #Inputs: # R - A rational difference equation # vars - variables in R # newvars - variables for the new rational difference eqn #Outputs: # If X is the positive equilibrium of R then the output is # the rational difference equation R with x[n] replaced by # z[n]+X (if z is the base of newvars) # The output format is a set of lists of length 2 in which # the first element is the equilibrium point of R and the # second element is the new rat. diff. eqn. which now # has equilibrium 0. MakeEquil0 := proc(R,vars,newvars) local x,Ret,X; X := FindPosEquil(R,vars); Ret := {}; for x in X do Ret := Ret union {[x,normal(subs({seq(vars[i]=newvars[i]+x,i=1..nops(vars))},R)-x)]}; od; RETURN(Ret); end: #========================================================== #FindPosEquil(R,vars) #Inputs: # R - a rational difference equation # vars - variables in R #Outputs: # The set of positive equilibrium points of R. This is # obtained by setting v=x for all v in vars and solving # for x. # The output format is a list of all positive equilibria FindPosEquil := proc(R,vars) local x,Rx,equil; if subs({seq(v=x,v in vars)},denom(R))=0 then RETURN([]); fi; Rx := subs({seq(v=x,v in vars)},R)-x; equil := [solve(Rx,x)]; equil := select(i->evalb(evalf(i)>=0),equil); RETURN([seq(simplify(e),e in equil)]); end: #========================================================== #LinStab(R,vars) #Inputs: # R - rational difference equation # vars - variables in R of the form [x[n-k],...,x[n-1],x[n]] #Outputs: # For each non-negative equilibrium checks the criteria given in # Camouzis and Ladas' book (Theorem 1.2.1) for an equilibria of # R to be locally asymptotically stable (LAS). The output format is # a set of lists of length two in which the first element is # the equilibrium and the second element is true if this equil # is LAS and false if this equil is not LAS. LinStab := proc(R,vars) local equil,xBar,k,lambda,Ret,i,q,charEq,ans; k := nops(vars)-1; #assume vars = [x[n-k],...,x[n-1],x[n]] Ret := {}; #Find xBar equil := FindPosEquil(R,vars); #print(equil); for xBar in equil do #Find q[i] = dF/dx[n-i](xBar,...,xBar) for i from 0 to k do q[i] := subs({seq(vars[i]=xBar,i=1..k+1)},diff(R,vars[k+1-i])); od; #print(op(q)); #Make characteristic equation # lambda^(k+1)-q[0]lambda^(k)-q[1]lambda^(k-1)-...-q[i]lambda^(k-i)-...-q[k-1]lambda-q[k] = 0 charEq := lambda^(k+1)-add(q[i]*lambda^(k-i),i=0..k); #print(xBar,charEq); #solve char. eqn. for lambda, check if all roots lie inside the open unit disk #ans := {solve(charEq,lambda)}; #Ret := Ret union {[xBar,ans]}; if k=0 then #print(ans); ans := {solve(charEq,lambda)}; Ret := Ret union {[xBar,seq(evalb(evalf(abs(ans[i]))<1),i=1..nops(ans))]}; elif k=1 then #use thm 1.2.2 Ret := Ret union {[xBar,Thm122(charEq,lambda)]}; elif k=2 then #use thm 1.2.3 Ret := Ret union {[xBar,Thm123(charEq,lambda)]}; elif k=3 then #use thm 1.2.4 Ret := Ret union {[xBar,Thm124(charEq,lambda)]}; else #use thm 1.2.5 Ret := Ret union {[xBar,Thm125(charEq,lambda)]}; fi; od; #RETURN({seq([r[1],{seq(evalb(evalf(abs(r[2][j]))<1),j=1..nops(r[2]))}],r in Ret)}); RETURN(Ret); end: #========================================================== #For order 2 equations Thm122 := proc(charEq,lambda) local a1,a0; a1 := coeff(charEq,lambda,1); a0 := coeff(charEq,lambda,0); if evalf(abs(a1))=0 for all newvars>=0. Input is allowed to have parameters # (FindPosEquil not called). OrigPolySymb := proc(R,vars,newvars,k,delta) option remember; local F1,Ret,f1,xBar,f; F1 := ToMinSymb(R,vars,newvars,k,delta); Ret := []; for f1 in F1 do xBar := f1[1]; f := subs({seq(newvars[i]=newvars[i]-xBar,i=1..nops(newvars))},f1[2]); Ret := [op(Ret), [xBar,f]]; od; RETURN(Ret); end: #========================================================== #ToMinSymb(R1,vars,newvars,k,delta) #Inputs: # -R A rational difference equation with equilibrium xBar # -vars List of variables in R in order [x[n-k],...,x[n]] # -newvars List of variables to replace vars when the equilibrium is moved to 0 # -k The k in |Q^k(newvars)|^2/|newvars|^2 # -delta We hope to prove |Q^k(newvars)|^2/|newvars|^2 <= delta<1 #Outputs: # The polynomial P=numer(delta)*|newvars|^2-|Q^k(newvars)|^2*denom(delta). We want # to prove that P>=0 for all newvars>=-xbar. Coefficients are allowed to be # parameters since FindPosEquil is not called. ToMinSymb := proc(R1,vars,newvars,k,delta) local iterK,n,R,Ret,r,Rat; R := MakeEquil0Symb(R1,vars,newvars); Ret := {}; for r in R do iterK := Iter(r[2],newvars,k); #print(iterK); n := nops(vars); Rat := add(iterK[i]^2,i=1..nops(iterK))/add(newvars[i]^2,i=1..n); Ret := Ret union {[r[1], -numer(Rat)*denom(delta)+denom(Rat)*numer(delta)]}; od; RETURN(Ret); end: #---------------------------------------------------------- #MakeEquil0Symb(R,vars,newvars): #Inputs: # R - A rational difference equation # vars - variables in R # newvars - variables for the new rational difference eqn #Outputs: # If X is the positive equilibrium of R then the output is # the rational difference equation R with x[n] replaced by # z[n]+X (if z is the base of newvars) # The output format is a set of lists of length 2 in which # the first element is the equilibrium point of R and the # second element is the new rat. diff. eqn. which now # has equilibrium 0. MakeEquil0Symb := proc(R,vars,newvars) local x,Ret,X; X := FindEquilSymb(R,vars); Ret := []; for x in X do #print(x); Ret := [op(Ret),[x,normal(subs({seq(vars[i]=newvars[i]+x,i=1..nops(vars))},R)-x)]]; od; RETURN(Ret); end: #---------------------------------------------------------- #FindEquilSymb(R,vars) #Inputs: # R - a rational difference equation # vars - variables in R #Outputs: # The set of equilibrium points of R. This is # obtained by setting v=x for all v in vars and solving # for x. # The output format is the list of all symbolic equilibria FindEquilSymb := proc(R,vars) local x,Rx,equil; if subs({seq(v=x,v in vars)},denom(R))=0 then RETURN([]); fi; Rx := subs({seq(v=x,v in vars)},R)-x; equil := [solve(Rx,x)]; RETURN([seq(simplify(e),e in equil)]); end: #========================================================== #LinStabSymb(R,vars,params) #Inputs: # R - rational difference equation # vars - variables in R of the form [x[n-k],...,x[n-1],x[n]] # params - the parameters that occur in R #Outputs: # For each non-negative equilibrium checks the criteria given in # Camouzis and Ladas' book (Theorem 1.2.1) for an equilibria of # R to be locally asymptotically stable (LAS). The output format is # a set of lists of length two in which the first element is # the equilibrium and the second element is the restrictions on # params for this equilibrium to be LAS. LinStabSymb := proc(R,vars,params) local equil,xBar,k,lambda,Ret,i,q,charEq,ans,toCheck,r,a0; k := nops(vars)-1; #assume vars = [x[n-k],...,x[n-1],x[n]] #toCheck := {}; Ret := {}; #Find xBar equil := FindEquilSymb(R,vars); #print(equil); for xBar in equil do #Find q[i] = dF/dx[n-i](xBar,...,xBar) for i from 0 to k do q[i] := subs({seq(vars[i]=xBar,i=1..k+1)},diff(R,vars[k+1-i])); od; #print(op(q)); #Make characteristic equation # lambda^(k+1)-q[0]lambda^(k)-q[1]lambda^(k-1)-...-q[i]lambda^(k-i)-...-q[k-1]lambda-q[k] = 0 charEq := lambda^(k+1)-add(q[i]*lambda^(k-i),i=0..k); charEq := simplify(charEq); #print(charEq); #solve char. eqn. for lambda #ans := {solve(charEq,lambda)}; #toCheck := toCheck union {[xBar,simplify(ans)]}; if k=0 then a0 := coeff(charEq,lambda,0); #print(a0); Ret := Ret union {[xBar,solve({abs(a0)<1},params)]}; elif k=1 then #use thm 1.2.2 Ret := Ret union {[xBar,Thm122Symb(charEq,lambda,params)]}; elif k=2 then #use thm 1.2.3 Ret := Ret union {[xBar,Thm123Symb(charEq,lambda,params)]}; elif k=3 then #use thm 1.2.4 Ret := Ret union {[xBar,Thm124Symb(charEq,lambda,params)]}; else #use thm 1.2.5 Ret := Ret union {[xBar,Thm125Symb(charEq,lambda,params)]}; fi; od; #Ret := {}; #for r in toCheck do # Ret := Ret union {[r[1],{solve({seq(abs(r[2][j])<1,j=1..nops(r[2]))},params)}]}; #od; RETURN(Ret); #RETURN({seq([r[1],{seq(evalb(evalf(abs(r[2][j]))<1),j=1..nops(r[2]))}],r in Ret)}); end: #========================================================== #For order 2 symbolic equations Thm122Symb := proc(charEq,lambda,params) local a1,a0,eq,ans,a,ansS; a1 := simplify(coeff(charEq,lambda,1)); a0 := simplify(coeff(charEq,lambda,0)); #print(a1,a0); eq := {simplify(abs(a1)<1+a0), simplify(1+a0<2)}; #print(eq); ans:= solve(eq,params); #print(ans); ansS := {}; for a in ans do ansS := ansS union {simplify(a)}; od; RETURN(ansS); end: #========================================================== #For order 3 symbolic equations Thm123Symb := proc(charEq,lambda,params) local a2,a1,a0,eq,ans,a,ansS; a2 := coeff(charEq,lambda,2); a1 := coeff(charEq,lambda,1); a0 := coeff(charEq,lambda,0); #print(a2,a1,a0); #print(abs(a2+a0)-(1+a1),abs(a2-3*a0)-(3-a1),a0^2+a1-a0*a2-1); eq := {simplify(abs(a2+a0)<(1+a1)), simplify(abs(a2-3*a0)<(3-a1)), simplify(a0^2+a1-a0*a2<1)}; #print(eq); ans := solve(eq,params); ansS := {}; for a in ans do ansS := ansS union {simplify(a)}; od; RETURN(ansS); end: #========================================================== #For order 4 symbolic equations Thm124Symb := proc(charEq,lambda,params) local a3,a2,a1,a0,eq,ans,a,ansS; a3 := coeff(charEq,lambda,3); a2 := coeff(charEq,lambda,2); a1 := coeff(charEq,lambda,1); a0 := coeff(charEq,lambda,0); eq :={abs(a1+a3)<(1+a0+a2), abs(a1-a3)<2*(1-a0),a2-3*a0<3, a0+a2+a0^2+a1^2+a0^2*a2+a0*a3^2<(1+2*a0*a2+a1*a3+a0*a1*a3+a0^3)}; #print(eq); ans := solve(eq,params); ansS := {}; for a in ans do ansS := ansS union {simplify(a)}; od; RETURN(ansS); end: #==================================================================== Help := proc() if args=NULL then print(The main procedures in this package are:); print(ProveNoPrint, ProveKNoPrint, Prove, ProveK, ProveTex, ProveAlmostNoPrintTex); print(); print(The other procedures for proving GAS are:); print(WebBookNoPrint ); print(PolynomialPositiveNoPrint ); print(PolynomialPositiveAbsNoPrint ); print(PolynomialPositiveParamsNoPrint ); print(WebBook ); print(PolynomialPositive ); print(WebBookTex ); print(PolynomialPositiveTex ); print(WebBookAlmostNoPrintTex ); print(); print(There are also supplementary procedures:); print(FindPosEquil, MakeEquil0, LinStab, Iter, OrigPoly, ToMin); print(FindEquilSymb, MakeEquil0Symb, OrigPolySymb, ToMinSymb); print(MakeSeq, Nlists, Nlists0, RandParams, GoodParams); print(); print(For help on a specific procedure type); print(Help();); else if args = WebBookNoPrint then lprint(WebBookNoPrint(R,vars,params,paramPoss,numParams,MinK,MaxK,delta,N)); lprint(Inputs:); lprint( R - a rational difference equation); lprint( vars - variables in R); lprint( params - parameters in R); lprint( paramPoss - possible values for parameters); lprint( numParams - number of parameter sets to test); lprint( MinK - the minimum K value to test); lprint( MaxK - the maximum K value to test); lprint( delta - the delta to test); lprint( N - a positive integer to be used when calling PolynomialPositive); lprint(Outputs:); lprint( List of lists of length 2. First element is set of parameters, second); lprint( element is output of Prove for subs(params,R).); lprint(Method:); lprint( Chooses parameters from paramPoss for which the equilibrium is ); lprint( rational if possible. If "not" possible (tries 10*numParams) then); lprint( checks other parameter sets.); lprint(Also:); lprint( "NoPrint" means that almost nothing is printed.); elif args = ProveKNoPrint then lprint(ProveKNoPrint(R,vars,K,delta,N)); lprint(Inputs:); lprint( R - a rational difference equation); lprint( vars - variables in R); lprint( K - K value to test); lprint( delta - delta value to test); lprint( N - positive integer to be used when calling PolynomialPositive); lprint(Outputs:); lprint( List of lists of length 2. First element is equilibrium point); lprint( second element is output of PolynomialPositive for that equil.); lprint(Also:); lprint( "NoPrint" means that almost nothing is printed.); lprint(Try:); lprint( ProveKNoPrint(1/(9/20+x[n]),[x[n]],2,1,2);); elif args = ProveNoPrint then lprint(ProveNoPrint(R,vars,MinK,MaxK,delta,N)); lprint(Inputs:); lprint( R - A rational difference equation); lprint( vars - variables in R); lprint( MinK - the minimum K value to try); lprint( MaxK - the maximum K value to try); lprint( delta - the delta to try); lprint( N - positive integer to be used in PolynomialPositive); lprint(Outputs:); lprint( A set of lists of length 2. For each list the first element is); lprint( the equilibrium value, the second element is the K value for that equ.); lprint( If K=-1 then that equilibrium is not LAS. If K=0 then the MaxK); lprint( value is not high enough.); lprint(Method:); lprint( For each equilibrium calls PolynomialPositive with OrigPoly for); lprint( that equilibrium.); lprint(Also:); lprint( "NoPrint" means that almost nothing is printed); lprint(Try:); lprint( ProveNoPrint(1/(9/20+x[n]),[x[n]],1,3,1,2);); elif args = PolynomialPositiveNoPrint then lprint(PolynomialPositiveNoPrint(P,Z,xbar,N)); lprint(Inputs:); lprint( P - a polynomial); lprint( Z - variables in P); lprint( xbar - a positive real number which marks the first cut of Rplus^|Z| into 2^|Z| regions); lprint( N - an integer greater than 1. Will be used when calling SmallBoxMethod if necessary); lprint(Outputs: ); lprint( Set containing elements from {true, false, FAIL}. ); lprint( if {true} then ); lprint( P(Z)>0 when Z[i]>=0 except possibly P(0,...,0)>=0 and P(xbar,...,xbar)>=0); lprint( elif false in output then); lprint( there is a point where P(Z')<=0 when Z'[i]>=0 (or P(0,...,0)<0 or P(xbar,...,xbar)<0 )); lprint( else (FAIL in output and false not in output)); lprint( there is some place where using IsPos is inconclusive); lprint(Method:); lprint( Cuts the positive orthant into 2^|Z| subsets and checks on each of ); lprint( these whether P(Z)>0. Each subset, I, of {1,...,|Z|} corresponds to the); lprint( subset of the positive orthant where 0<=Z[i]<=xbar if i in I, ); lprint( and xbar<=Z[j] if j not in I. ); lprint(Also:); lprint( "NoPrint" means that almost nothing is printed in this implementation ); lprint( of PolynomialPositive); elif args = PolynomialPositiveAbsNoPrint then lprint(PolynomialPositiveAbsNoPrint(P,Z,xbar,N)); lprint(Inputs:); lprint( P - a polynomial); lprint( Z - variables in P); lprint( xbar - a positive real number which marks the first cut of Rplus^|Z| into 2^|Z| regions); lprint( N - an integer greater than 1. Will be used when calling SmallBoxMethod if necessary); lprint(Outputs: ); lprint( Set containing elements from {true, false, FAIL}. ); lprint( if {true} then ); lprint( P(Z)>0 when Z[i]>=0 except possibly P(0,...,0)>=0 and P(xbar,...,xbar)>=0); lprint( elif false in output then); lprint( there is a point where P(Z')<=0 when Z'[i]>=0 (or P(0,...,0)<0 or P(xbar,...,xbar)<0 )); lprint( else (FAIL in output and false not in output)); lprint( there is some place where using IsPos is inconclusive); lprint(Method:); lprint( Cuts the positive orthant into 2^|Z| subsets and checks on each of ); lprint( these whether P(Z)>0. Each subset, I, of {1,...,|Z|} corresponds to the); lprint( subset of the positive orthant where 0<=Z[i]<=xbar if i in I, ); lprint( and xbar<=Z[j] if j not in I. ); lprint(Also:); lprint( "AbsNoPrint" means that absolutely nothing is printed in this implementation ); lprint( of PolynomialPositive); elif args = PolynomialPositiveParamsNoPrint then lprint(PolynomialPositiveParamsNoPrint(P,Z,params,xbar,N)); lprint(Inputs:); lprint( P - a polynomial); lprint( Z - variables in P); lprint( params - parameters in P); lprint( xbar - a positive real number which marks the first cut of Rplus^|Z| into 2^|Z| regions); lprint( N - an integer greater than 1. Will be used when calling SmallBoxMethod if necessary); lprint(Outputs: ); lprint( Set containing elements from {true, false, FAIL}. ); lprint( if {true} then ); lprint( P(Z)>0 when Z[i]>=0 and params >=0 except possibly ); lprint( P(0,...,0)>=0 and P(xbar,...,xbar)>=0); lprint( elif false in output then); lprint( there is a point where P(Z')<=0 when Z'[i]>=0 ); lprint( (or P(0,...,0)<0 or P(xbar,...,xbar)<0 )); lprint( else (FAIL in output and false not in output)); lprint( there is some place where using IsPos is inconclusive); lprint(Method:); lprint( Cuts the positive orthant into 2^|Z| subsets and checks on each of ); lprint( these whether P(Z)>0. Each subset, I, of {1,...,|Z|} corresponds to the); lprint( subset of the positive orthant where 0<=Z[i]<=xbar if i in I, ); lprint( and xbar<=Z[j] if j not in I. Treats parameters differently. No cutting); lprint( of parameter space, always just assuming the params are <= 0.); lprint(Also:); lprint( "ParamsNoPrint" means that almost nothing is printed in this implementation); lprint( of PolynomialPositive and that parameters are treated differently than variables); elif args = PolynomialPositive then lprint(PolynomialPositive(P,Z,xbar,N)); lprint(Inputs:); lprint( P - a polynomial); lprint( Z - variables in P); lprint( xbar - a positive real number which marks the first cut of Rplus^|Z| into 2^|Z| regions); lprint( N - an integer greater than 1. Will be used when calling SmallBoxMethod if necessary); lprint(Outputs: ); lprint( Set containing elements from {true, false, FAIL}. ); lprint( if {true} then ); lprint( P(Z)>0 when Z[i]>=0 except possibly P(0,...,0)>=0 and P(xbar,...,xbar)>=0); lprint( elif false in output then); lprint( there is a point where P(Z')<=0 when Z'[i]>=0 (or P(0,...,0)<0 or P(xbar,...,xbar)<0 )); lprint( else (FAIL in output and false not in output)); lprint( there is some place where using IsPos is inconclusive); lprint(Method:); lprint( Cuts the positive orthant into 2^|Z| subsets and checks on each of ); lprint( these whether P(Z)>0. Each subset, I, of {1,...,|Z|} corresponds to the); lprint( subset of the positive orthant where 0<=Z[i]<=xbar if i in I, ); lprint( and xbar<=Z[j] if j not in I. ); lprint(Also:); lprint( There are print statements in this implementation of PolynomialPositive); lprint(Try: ); lprint( PolynomialPositive(x^2-x*y+y^2,[x,y],1,2);); elif args = WebBook then lprint(WebBook(R,vars,params,paramPoss,numParams,MinK,MaxK,delta,N)); lprint(Inputs:); lprint( R - a rational difference equation); lprint( vars - variables in R); lprint( params - parameters in R); lprint( paramPoss - possible values for parameters); lprint( numParams - number of parameter sets to test); lprint( MinK - the minimum K value to test); lprint( MaxK - the maximum K value to test); lprint( delta - the delta to test); lprint( N - a positive integer to be used when calling PolynomialPositive); lprint(Outputs:); lprint( List of lists of length 2. First element is set of parameters, second); lprint( element is output of Prove for subs(params,R).); lprint(Method:); lprint( Chooses parameters from paramPoss for which the equilibrium is ); lprint( rational if possible. If "not" possible (tries 10*numParams) then); lprint( checks other parameter sets.); lprint(Also:); lprint( There are print statements throughout.); elif args = ProveK then lprint(ProveK(R,vars,K,delta,N)); lprint(Inputs:); lprint( R - a rational difference equation); lprint( vars - variables in R); lprint( K - K value to test); lprint( delta - delta value to test); lprint( N - positive integer to be used when calling PolynomialPositive); lprint(Outputs:); lprint( List of lists of length 2. First element is equilibrium point); lprint( second element is output of PolynomialPositive for that equil.); lprint(Also:); lprint( There are print statements throughout.); lprint(Try:); lprint( ProveK((4+x[n])/(1+x[n-1]),[x[n-1],x[n]],5,1,2);); elif args = Prove then lprint(Prove(R,vars,MinK,MaxK,delta,N)); lprint(Inputs:); lprint( R - A rational difference equation); lprint( vars - variables in R); lprint( MinK - the minimum K value to try); lprint( MaxK - the maximum K value to try); lprint( delta - the delta to try); lprint( N - positive integer to be used in PolynomialPositive); lprint(Outputs:); lprint( A set of lists of length 2. For each list the first element is); lprint( the equilibrium value, the second element is the K value for that equ.); lprint( If K=-1 then that equilibrium is not LAS. If K=0 then the MaxK); lprint( value is not high enough.); lprint(Method:); lprint( For each equilibrium calls PolynomialPositive with OrigPoly for); lprint( that equilibrium.); lprint(Try:); lprint( Prove((4+x[n])/(1+x[n-1]),[x[n-1],x[n]],1,5,1,2);); elif args = PolynomialPositiveTex then lprint(PolynomialPositiveTex(P,Z,xbar,N)); lprint(Inputs:); lprint( P - a polynomial); lprint( Z - variables in P); lprint( xbar - a positive real number which marks the first cut of Rplus^|Z| into 2^|Z| regions); lprint( N - an integer greater than 1. Will be used when calling SmallBoxMethod if necessary); lprint(Outputs: ); lprint( Set containing elements from {true, false, FAIL}. ); lprint( if {true} then ); lprint( P(Z)>0 when Z[i]>=0 except possibly P(0,...,0)>=0 and P(xbar,...,xbar)>=0); lprint( elif false in output then); lprint( there is a point where P(Z')<=0 when Z'[i]>=0 (or P(0,...,0)<0 or P(xbar,...,xbar)<0 )); lprint( else (FAIL in output and false not in output)); lprint( there is some place where using IsPos is inconclusive); lprint(Method:); lprint( Cuts the positive orthant into 2^|Z| subsets and checks on each of ); lprint( these whether P(Z)>0. Each subset, I, of {1,...,|Z|} corresponds to the); lprint( subset of the positive orthant where 0<=Z[i]<=xbar if i in I, ); lprint( and xbar<=Z[j] if j not in I. ); lprint(Also:); lprint( "Tex" means that it prints what it is doing as we go along in ); lprint( such a way that the resulting text can be TeXed); elif args = WebBookTex then lprint(WebBookTex(R,vars,params,paramPoss,numParams,MinK,MaxK,delta,N)); lprint(Inputs:); lprint( R - a rational difference equation); lprint( vars - variables in R); lprint( params - parameters in R); lprint( paramPoss - possible values for parameters); lprint( numParams - number of parameter sets to test); lprint( MinK - the minimum K value to test); lprint( MaxK - the maximum K value to test); lprint( delta - the delta to test); lprint( N - a positive integer to be used when calling PolynomialPositive); lprint(Outputs:); lprint( List of lists of length 2. First element is set of parameters, second); lprint( element is output of Prove for subs(params,R).); lprint(Method:); lprint( Chooses parameters from paramPoss for which the equilibrium is ); lprint( rational if possible. If "not" possible (tries 10*numParams) then); lprint( checks other parameter sets.); lprint(Also:); lprint( "Tex" means that it prints what it is doing as we go along in ); lprint( such a way that the resulting text can be TeXed); elif args = ProveTex then lprint(ProveTex(R,vars,MinK,MaxK,delta,N)); lprint(Inputs:); lprint( R - A rational difference equation); lprint( vars - variables in R); lprint( MinK - the minimum K value to try); lprint( MaxK - the maximum K value to try); lprint( delta - the delta to try); lprint( N - positive integer to be used in PolynomialPositive); lprint(Outputs:); lprint( A set of lists of length 2. For each list the first element is); lprint( the equilibrium value, the second element is the K value for that equ.); lprint( If K=-1 then that equilibrium is not LAS. If K=0 then the MaxK); lprint( value is not high enough.); lprint(Method:); lprint( For each equilibrium calls PolynomialPositive with OrigPoly for); lprint( that equilibrium.); lprint(Also:); lprint( "Tex" means that it prints what it is doing as we go along in ); lprint( such a way that the resulting text can be TeXed); lprint(Try:); lprint( ProveTex((4+x[n])/(1+x[n-1]),[x[n-1],x[n]],1,5,1,2);); elif args = WebBookAlmostNoPrintTex then lprint(WebBookAlmostNoPrintTex(R,vars,params,paramPoss,numParams,MinK,MaxK,delta,N)); lprint(Inputs:); lprint( R - a rational difference equation); lprint( vars - variables in R); lprint( params - parameters in R); lprint( paramPoss - possible values for parameters); lprint( numParams - number of parameter sets to test); lprint( MinK - the minimum K value to test); lprint( MaxK - the maximum K value to test); lprint( delta - the delta to test); lprint( N - a positive integer to be used when calling PolynomialPositive); lprint(Outputs:); lprint( List of lists of length 2. First element is set of parameters, second); lprint( element is output of Prove for subs(params,R).); lprint(Method:); lprint( Chooses parameters from paramPoss for which the equilibrium is ); lprint( rational if possible. If "not" possible (tries 10*numParams) then); lprint( checks other parameter sets.); lprint(Also:); lprint( "AlmostNoPrintTex" means that it prints some of what it is doing as we go along in ); lprint( such a way that the resulting text can be TeXed); elif args = ProveAlmostNoPrintTex then lprint(ProveAlmostNoPrintTex(R,vars,MinK,MaxK,delta,N)); lprint(Inputs:); lprint( R - A rational difference equation); lprint( vars - variables in R); lprint( MinK - the minimum K value to try); lprint( MaxK - the maximum K value to try); lprint( delta - the delta to try); lprint( N - positive integer to be used in PolynomialPositive); lprint(Outputs:); lprint( A set of lists of length 2. For each list the first element is); lprint( the equilibrium value, the second element is the K value for that equ.); lprint( If K=-1 then that equilibrium is not LAS. If K=0 then the MaxK); lprint( value is not high enough.); lprint(Method:); lprint( For each equilibrium calls PolynomialPositive with OrigPoly for); lprint( that equilibrium.); lprint(Also:); lprint( "AlmostNoPrintTex" means that it prints some of what it is doing); lprint( as we go along in such a way that the resulting text can be TeXed); lprint(Try:); lprint( ProveAlmostNoPrintTex((4+x[n])/(1+x[n-1]),[x[n-1],x[n]],1,5,1,2);); elif args = GoodParams then lprint(GoodParams(R,vars,params,possParams)); lprint(Inputs:); lprint( R - a rational difference equation); lprint( vars - variables in R); lprint( params - parameters in R); lprint( possParams - possible values for elements of params); lprint(Outputs:); lprint( A random parameter assignment such that the positive ); lprint( equilibrium of R with this parameter assignment is rational.); elif args = RandParams then lprint(RandParams(R,vars,params,possParams)); lprint(Inputs:); lprint( R - a rational difference equation); lprint( vars - variables in R); lprint( params - parameters in R); lprint( possParams - possible values for elements of params); lprint(Outputs:); lprint( A random parameter assignment from the set possParams); elif args = Nlists then lprint(Nlists(n,N)); lprint(Inputs:); lprint( n - a positive integer); lprint( N - a positive integer); lprint(Outputs:); lprint( Lists of length n consisting of integers between 1 and N (inclusively)); elif args = Nlists0 then lprint(Nlists0(n,N)); lprint(Inputs:); lprint( n - a positive integer); lprint( N - a positive integer); lprint(Outputs:); lprint( Lists of length n consisting of integers between 0 and N (inclusively)); elif args = OrigPoly then lprint(OrigPoly(R,vars,newvars,k,delta)); lprint(Inputs:); lprint( -R A rational difference equation with equilibrium xBar); lprint( -vars List of variables in R in order [x[n-k],...,x[n]]); lprint( -newvars List of variables to replace vars when the equilibrium is moved to 0); lprint( -k The k in |Q^k(newvars)|^2/|newvars|^2); lprint( -delta We hope to prove |Q^k(newvars)|^2/|newvars|^2 <= delta<1); lprint(Outputs:); lprint( The polynomial P=numer(delta)*|newvars|^2-|Q^k(newvars)|^2*denom(delta)); lprint( shifted to the right in each variable by xBar. We want to prove that ); lprint( P>=0 for all newvars>=0); lprint(Try:); lprint( OrigPoly((4+x[n])/(1+x[n-1]),[x[n-1],x[n]],[z[n-1],z[n]],5,1);); elif args = ToMin then lprint(ToMin(R1,vars,newvars,k,delta)); lprint(Inputs:); lprint( -R A rational difference equation with equilibrium xBar); lprint( -vars List of variables in R in order [x[n-k],...,x[n]]); lprint( -newvars List of variables to replace vars when the equilibrium is moved to 0); lprint( -k The k in |Q^k(newvars)|^2/|newvars|^2); lprint( -delta We hope to prove |Q^k(newvars)|^2/|newvars|^2 <= delta<1); lprint(Outputs:); lprint( The polynomial P=numer(delta)*|newvars|^2-|Q^k(newvars)|^2*denom(delta). We want); lprint( to prove that P>=0 for all newvars>=-xbar. Coefficients must be numerical); lprint( since FindPosEquil is called (within MakeEquil0).); lprint(Try:); lprint( ToMin((4+x[n])/(1+x[n-1]),[x[n-1],x[n]],[z[n-1],z[n]],5,1);); elif args = Iter then lprint(Iter(R,vars,k):); lprint(Inputs:); lprint( R - a rational difference equation); lprint( vars - variables in R); lprint( k - number of times to iterate R); lprint(Outputs:); lprint( Iterates the map); lprint( T:(vars)->(vars[2..nops(vars)],R(vars))); lprint( k times. The output is the vector, T^k(vars), of length nops(vars).); lprint(Try:); lprint( Iter((4+x[n])/(1+x[n-1]),[x[n-1],x[n]],2);); elif args = MakeEquil0 then lprint(MakeEquil0(R,vars,newvars):); lprint(Inputs:); lprint( R - A rational difference equation); lprint( vars - variables in R); lprint( newvars - variables for the new rational difference eqn); lprint(Outputs:); lprint( If X is the positive equilibrium of R then the output is); lprint( the rational difference equation R with x[n] replaced by); lprint( z[n]+X (if z is the base of newvars)); lprint( The output format is a set of lists of length 2 in which ); lprint( the first element is the equilibrium point of R and the); lprint( second element is the new rat. diff. eqn. which now ); lprint( has equilibrium 0.); lprint(Try:); lprint( MakeEquil0((4+x[n])/(1+x[n-1]),[x[n-1],x[n]],[z[n-1],z[n]]);); elif args = FindPosEquil then lprint(FindPosEquil(R,vars)); lprint(Inputs:); lprint( R - a rational difference equation); lprint( vars - variables in R); lprint(Outputs:); lprint( The set of positive equilibrium points of R. This is ); lprint( obtained by setting v=x for all v in vars and solving); lprint( for x.); lprint( The output format is a list of all positive equilibria); lprint(Try:); lprint( FindPosEquil((4+x[n])/(1+x[n-1]),[x[n-1],x[n]]);); elif args = LinStab then lprint(LinStab(R,vars)); lprint(Inputs:); lprint( R - rational difference equation); lprint( vars - variables in R of the form [x[n-k],...,x[n-1],x[n]]); lprint(Outputs:); lprint( For each non-negative equilibrium checks the criteria given in ); lprint( Camouzis and Ladas' book (Theorem 1.2.1) for an equilibria of ); lprint( R to be locally asymptotically stable (LAS). The output format is); lprint( a set of lists of length two in which the first element is); lprint( the equilibrium and the second element is true if this equil); lprint( is LAS and false if this equil is not LAS.); lprint(Try:); lprint( LinStab((4+x[n])/(1+x[n-1]),[x[n-1],x[n]]);); elif args = OrigPolySymb then lprint(OrigPolySymb(R,vars,newvars,k,delta)); lprint(Inputs:); lprint( -R A rational difference equation with equilibrium xBar); lprint( -vars List of variables in R in order [x[n-k],...,x[n]]); lprint( -newvars List of variables to replace vars when the equilibrium is moved to 0); lprint( -k The k in |Q^k(newvars)|^2/|newvars|^2); lprint( -delta We hope to prove |Q^k(newvars)|^2/|newvars|^2 <= delta<1); lprint(Outputs:); lprint( The polynomial P=numer(delta)*|newvars|^2-|Q^k(newvars)|^2*denom(delta)); lprint( shifted to the right in each variable by xBar. We want to prove that ); lprint( P>=0 for all newvars>=0. Input is allowed to have parameters ); lprint( (FindPosEquil not called).); lprint(Try:); lprint( OrigPolySymb(1/(A+x[n-1]),[x[n-1],x[n]],[z[n-1],z[n]],2,1);); elif args = ToMinSymb then lprint(ToMinSymb(R1,vars,newvars,k,delta)); lprint(Inputs:); lprint( -R A rational difference equation with equilibrium xBar); lprint( -vars List of variables in R in order [x[n-k],...,x[n]]); lprint( -newvars List of variables to replace vars when the equilibrium is moved to 0); lprint( -k The k in |Q^k(newvars)|^2/|newvars|^2); lprint( -delta We hope to prove |Q^k(newvars)|^2/|newvars|^2 <= delta<1); lprint(Outputs:); lprint( The polynomial P=numer(delta)*|newvars|^2-|Q^k(newvars)|^2*denom(delta). We want); lprint( to prove that P>=0 for all newvars>=-xbar. Coefficients are allowed to be ); lprint( parameters since FindPosEquil is not called.); lprint(Try:); lprint( ToMinSymb(1/(A+x[n-1]),[x[n-1],x[n]],[z[n-1],z[n]],2,1);); elif args = MakeEquil0Symb then lprint(MakeEquil0Symb(R,vars,newvars):); lprint(Inputs:); lprint( R - A rational difference equation); lprint( vars - variables in R); lprint( newvars - variables for the new rational difference eqn); lprint(Outputs:); lprint( If X is the positive equilibrium of R then the output is); lprint( the rational difference equation R with x[n] replaced by); lprint( z[n]+X (if z is the base of newvars)); lprint( The output format is a set of lists of length 2 in which ); lprint( the first element is the equilibrium point of R and the); lprint( second element is the new rat. diff. eqn. which now ); lprint( has equilibrium 0.); lprint(Try:); lprint( MakeEquil0Symb((a+x[n])/(A+x[n-1]),[x[n-1],x[n]],[z[n-1],z[n]]);); elif args = FindEquilSymb then lprint(FindEquilSymb(R,vars)); lprint(Inputs:); lprint( R - a rational difference equation); lprint( vars - variables in R); lprint(Outputs:); lprint( The set of equilibrium points of R. This is ); lprint( obtained by setting v=x for all v in vars and solving); lprint( for x.); lprint( The output format is the list of all symbolic equilibria); lprint(Try:); lprint( FindEquilSymb((a+x[n])/(A+x[n-1]),[x[n-1],x[n]]);); elif args = MakeSeq then lprint(MakeSeq(R,vars,init,K)); lprint(Inputs:); lprint( R - a rational difference equation); lprint( vars - variables in R); lprint( init - initial values for the vars); lprint( K - number of iterations); lprint(Outputs:); lprint( The sequence of length K generated by x[n+1]=R (assuming); lprint( vars=[x[n-k],...,x[n]]) with initial conditions init.); lprint(Try:); lprint( MakeSeq(1/(9/20+x[n-1]),[x[n-1],x[n]],[1,2],5);); fi; fi; end: