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 $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: #==================================================================== #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 prove $h>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 that $H^{\\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: