######################################################################
##AmalGerry.txt: Save this file as  AmalGerry.txt                    #
## To use it, stay in the                                            #
##same directory, get into Maple (by typing: maple <Enter> )         #
##and then type:  read `AmalGerry.txt`<Enter>                        #
##Then follow the instructions given there                           #
##                                                                   #
##Written by George Spahn and                                        #
#Doron Zeilberger, Rutgers University ,                              #
#DoronZeil at gmail dot com                                          #
######################################################################
 
#Created: Spring 2023
 
Digits:=50:

print(`Created:  Spring  2023`):
print(` This is  AmalGerry.txt `):
print(`It is one of the packages that accompany the article `):
print(`Experimenting with Discrete Rational Dynamical Systems`):
print(`by George Spahn and Doron Zeilberger`):

print(`and also available from Zeilberger's website`):
print(``):
print(`Please report bugs to DoronZeil at gmail dot com `):
print(``):
 print(`The most current version of this  package and paper`):
 print(` are  available from`):
 print(`http://sites.math.rutgers.edu/~zeilberg/  .`):
print(``):

print(`The main procedures are run_nlp and run_symb`):
print(`For more info, type ezra(run_nlp) or ezra(run_symb)`):

ezra := proc()
if op(1,[args])=run_nlp then
printf("run_nlp(T,v,m,n,gran) runs a series of tests to determine if T will always converge to points in R3 with small v-norm.\n\
T is a transformation on 3 variables. For example try print(T4), which gives the sequence corresponding to equation 4 in Ladas's paper.\n\
v is a norm on 3 variables, should be always positive. Try print(v5).\n\
The tests compare m iterations of T to n iterations of T. The objective function is given by F(x) = v(T^m(x)) - v(T^n(x))\n\
Maple attempts to minimize this objective function using the numerical procedure 'NLPSolve'.\n\
If the minimum is no less than -10^-6, the program will output that 'All tests passed'\n\
The argument gran specifies how many different initial points will be generated for NLPSolve. gran=3 corresponds to 27 initial points, spaced in a cube.\n\
Results so far verify conjectures 1-4. Try:\n\
run_nlp(T4,v4,3,7,3)\n\
run_nlp(T5,v5,3,7,3)\n\
run_nlp(T6,v6,3,8,3)\n\
run_nlp(T7,v7,3,7,3)\n"):
fi:

if op(1,[args])=run_symb then
printf("run_symb(T,v,n) attempts to prove rigorously that T will always converge to points in R3 with small v-norm.\n\
T is a transformation on 3 variables. For example try print(T4), which gives the sequence corresponding to equation 4 in Ladas's paper.\n\
v is a norm on 3 variables, should be always positive. Try print(v5).\n\
The objective function is given by F(x) = v(x) - v(T^n(x))\n\
Maple attempts to minimize this objective function using the symbolic procedure 'minimize'.\n\
The program will output the denominator and the minimum of the numerator.\n\
In most cases the denominator is nonnegative."):
print(`In order to get a fully rigorous proof of Conjecture 1 in the paper. Type:`):
print(`run_symb(T4,v4A,4);`):
print(`run_symb(T4,v4B,4);`):
print(`Note that v4A is what is called v_{xy} in the paper, and v4B is alled v_{yz}`):
fi:
end:

T4 := proc(x,y,z)
	y,z,evalf(y/(x+y)):
end:

v4 := proc(x,y,z)
	1+x*x+y*y+z*z-x-2*y-z+x*y-x*z+y*z:
end:
#run_nlp01(T4,v4,4,3)
v4A := proc(x,y,z)
	(y+x-1)*(y+x-1):
end:
v4B := proc(x,y,z)
	(y+z-1)*(y+z-1):
end:


T5 := proc(x,y,z)
	y,z,evalf((z+x)/y):
end:


v5 := proc(x,y,z)
	(z - (x+y*y)/(x*y-1))^2:
end:

#run_nlp(T5,v5,3,7,3)

T6 := proc(x,y,z)
	y,z,evalf((1+x)/z):
end:

v6 := proc(x,y,z)
	(z - (x+1)/(x*y-1))^2:
end:
#run_nlp(T6,v6,3,8,3)


T7 := proc(x,y,z)
	y,z,evalf((1+z)/(x+y)):
end:

v7 := proc(x,y,z)
	(z - y/x)^2:
end:

#run_nlp(T7,v7,3,7,3)

T10 := proc(x,y,z)
	y,z,evalf((1+y)/(1+z)):
end:

v10 := proc(x,y,z)
	(z - 1/y)^2:
end:
#run_nlp(T10,v10,3,5,3)

T11 := proc(x,y,z)
	y,z,evalf((y+z)/(x+z)):
end:

v11 := proc(x,y,z)
	(z - y/(2*y-1))^2:
end:

T12 := proc(x,y,z)
	y,z,evalf((x+y)/(x+z)):
end:

v12 := proc(x,y,z)
	(z - y/(2*y-1))^2:
end:


T13 := proc(x,y,z)
	y,z,evalf((x+y)/(x)):
end:

v13 := proc(x,y,z)
	(z - y/(y-1))^2:
end:
#run_nlp(T13,v13,3,7,3)

T14 := proc(x,y,z)
	y,z,evalf((1+y)/(1+x+z)):
end:

v14 := proc(x,y,z)
	(z - .5/y)^2:
end:

T15 := proc(x,y,z)
	y,z,evalf((1+x+y)/(x)):
end:

v15 := proc(x,y,z)
	(z - (y+1)/(y-1))^2:
end:



run_test01 := proc(T,v,n,gran)
local i,j,k,l,a,b,c:

for i from 1 to gran do
for j from 1 to gran do
for k from 1 to gran do

	a,b,c := i/gran,j/gran,k/gran:
	for l from 1 to n do
		a,b,c := T(a,b,c):
	od:
	
	if v(i/gran,j/gran,k/gran) < v(a,b,c) then
		print	(i/gran,j/gran,k/gran,v(i/gran,j/gran,k/gran) - v(a,b,c)):
	fi:
od:
od:
od:
end:

#compare m iterations to n iterations
run_test := proc(T,v,m,n,gran)
local i,j,k,l,a,b,c,a1,b1,c1:

for i from 2 to gran do
for j from 2 to gran do
for k from 2 to gran do

	a,b,c := i,j,k:
	for l from 1 to m do
		a,b,c := T(a,b,c):
	od:
	a1,b1,c1 := a,b,c:
	for l from 1 to n-m do
		a,b,c := T(a,b,c):
	od:
	
	if v(a1,b1,c1) < v(a,b,c) then
		print	(a1,b1,c1,v(i,j,k) - v(a,b,c)):
	fi:
od:
od:
od:
end:


run_nlp01 := proc(T,v,n,gran)
local i,j,k,l,a,b,c,r,fs,i1,j1,k1:


a,b,c := x,y,z:
for l from 1 to n do
	a,b,c := T(a,b,c):
od:
fs := v(x,y,z)-v(a,b,c):
print(fs):
for i from 1 to gran do
for j from 1 to gran do
for k from 1 to gran do


	i1 := (i-.51)/gran:
	j1 := (j-.52)/gran:
	k1 := (k-.53)/gran:
	
	r:=Optimization:-NLPSolve(fs,x=0..1,y=0..1,z=0..1,initialpoint=[x=i1,y=j1,z=k1]):

	print((i1,j1,k1),r):
od:
od:
od:
end:

run_nlp := proc(T,v,m,n,gran)
local i,j,k,l,a,b,c,r,fs,i1,j1,k1,a1,b1,c1,fl:	

fl := true:
a,b,c := x,y,z:
for l from 1 to n do
	a,b,c := T(a,b,c):
od:
a1,b1,c1 := a,b,c:
	for l from 1 to n-m do
		a,b,c := T(a,b,c):
	od:
fs := v(a1,b1,c1)-v(a,b,c):
#print(fs):
for i from 1 to gran do
for j from 1 to gran do
for k from 1 to gran do


	i1 := (i-.51):
	j1 := (j-.52):
	k1 := (k-.53):
	
	r:=Optimization:-NLPSolve(fs,x=0.1..10,y=0.1..10,z=0.1..10,initialpoint=[x=i1,y=j1,z=k1]):
	
	if r[1] < -10^(-6) then
		print((i1,j1,k1),r):
		fl := false:
	fi:
od:
od:
od:
if fl then
	print("All tests passed"):
fi:

end:


appl := proc(T,x,y,z,n) 
local i,a,b,c:
a,b,c := x,y,z:
for i from 1 to n do
  a,b,c := T(a,b,c)
od:
normal(a),normal(b),normal(c):
end:

run_symb := proc(T,v,n)
local F1,m:
F1 := normal( v(x,y,z) - v(appl(T,x,y,z,n)) ):
print("denominator:",denom(F1)):
m := minimize(numer(F1),x=0..infinity,y=0..infinity,z=0..infinity):
print("min of numerator:",m):
end: