""" Computations concerning the solutions of the Calogero-Moser model using the root-type Lax pair. Revision 1, 30 Nov 2011 This code was used for [1]. It can be interpreted by the Sage mathematical system, which can be found at [2]. It is not neccessary to have a copy installed; Sage can be used online at [3] -- but please do not do heavy calculations on that public server. GETTING STARTED To verify the results in [1], first load this file in sage by load("/path/to/calogero-moser-root-type-lax-pair.sage") (substitute the path on your computer for /path/to). Or, if you are not running Sage on your own computer (but online on [3], say), then copy-and-paste the contents of this file into a cell and click 'evaluate'. Then, run the following commands: To see a plot of the example system mentioned in section 5.2: image = plot_example(0,-2,2,100,"blue",*examples[0]) image.show(ymin=-3, ymax=3) Note that there are dents where the solution crosses the boundary of a fundamental domain, because we haven't calculated the indicator function yet. The discriminant that is mentioned there is obtained by: f = factor(get_indicator_discriminant(*examples[0])) f We can use the indicator function thus obtained by: c = f[0][0] image = plot_example(0,-2,2,100,"blue",*examples[0], indicator_function = c) image.show(ymin=-3, ymax=3) (Here, f[0] gives the first factor, and f[0][0] gives the factor itself. f[0][1] contains the multiplicity). You can tick the "Typeset" checkbox at the top of the page to make it look pretty. The search for counterexamples to the generic statement of proposition 3.1 was done by find_problems(RootSystem(['A',5]), 0, 6) Of course, the point in publishing this code is that everything is verifiable, so you are invited to look at the implementation of all these functions. I have tried to make it as modular and descriptive as possible so that you can 'play' with it. REFERENCES [1] Timo Kluck, On the Calogero-Moser solution by root-type Lax pair, http://arxiv.org/{to-be-submitted} [2] www.sagemath.org [3] www.sagenb.org COPYING Copyright (C) Universiteit Utrecht 2011 Copyright (C) Timo Kluck 2011 This is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. It is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with this program. If not, see . """ """ Some data describing examples The data contains a root system together with initial values for p,q, a value for g, and a vector `chi` whose orthoplement separates two fundamental domains in a Weyl chamber. """ import numpy import numpy.linalg examples = [ (RootSystem(['A',2]), [-10, 10, 0], [6/10,-1/10,-1/2], 1, (0,-1,0)), (RootSystem(['A',5]), [0,0,0,0,0,0], [-28,-22,-16,8,20,38], 1, (0,0,1,1,0,0)), (RootSystem(['A',5]), [0,0,0,0,0,0], [-34,-28,2,8,20,32], 1, (0,0,1,1,0,0)) ] def do_example(initial_value_time, minimal_time, maximal_time, number_of_intervals, root_system,p0_value,q0_value,g_value,chi, indicator_function = lambda t: 1): """ Calculate the time evolution of a C-M system described by the parameters USAGE: q = do_example(initial_value_time, minimal_time, maximal_time, time_step, root_system, p0, q0, g, chi) or q = do_example(initial_value_time, minimal_time, maximal_time, *example[i]) to do one of the pre-defined examples (currently, i can be between 0 and 3 (excl). It is not necessary that the initial_value_time is in between of minimal_time and maximal_time. The return value is a dict matching times to values of q. """ number_of_particles = len(root_system.index_set()) if root_system.cartan_type()[0] == 'A': number_of_particles += 1 p,q,p0,q0,g,t = define_variables(number_of_particles, False) L,M,HH = lax_pair(root_system, p, q, g) W = lax_W(root_system,L,M,p,p0_value,q,q0_value,g,g_value,t) delta_t = (maximal_time - minimal_time) / number_of_intervals q_t = {} for time in [minimal_time + delta_t * n for n in range(number_of_intervals+1)]: q_t[time] = calculate_q_at_time(root_system, W, time - initial_value_time, q, chi, indicator_function(time), t, 10^(-3)) return q_t def plot_example(initial_value_time, minimal_time, maximal_time, number_of_intervals, color, root_system, p0_value, q0_value, g_value, chi, indicator_function = lambda t: 1): qs = do_example(initial_value_time, minimal_time, maximal_time, number_of_intervals, root_system, p0_value, q0_value, g_value, chi, indicator_function) times = sorted(k for k in qs.keys() if qs[k] != []) times_with_no_data = sorted(k for k in qs.keys() if qs[k] == []) if(len(times_with_no_data) > 0): print("Warning: no data at times:") print(times_with_no_data) return worldline_plot(times,[qs[t][0] for t in times], color) def worldline_plot(times,qs,color): if(len(times) != len(qs)): raise ValueError("should be as many times as qs") if(len(set(len(q) for q in qs)) != 1): raise ValueError("all qs should have the same number of particle coordinates") lines = plot([]) # for every particle for i in range(len(qs[0])): coords = [qs[t][i] for t in range(len(times))] lines += line(zip(times,coords),color=color) return lines def define_variables(number_of_particles, use_symbolic_ring = false): """ Define the dynamical variables for the C-M system USAGE: p,q,p0,q0,g,t = define_variables(number_of_particles, use_symbolic_ring) The return value p is a list of variables p_1 up to p_n, where n is `number_of_particles`, and similarly for q, p0, q0. The variable `g` is the coupling constant and `t` is time. The variables can either be defined in the symbolic ring or, trading of readability of expressions for speed, in a polynomial ring. In the latter case, we inject the variables into the global namespace so that they can be used regardless. This includes the formal variable I for the square root of -1. """ # canonical variables variable_names=['p%i'%i for i in [1..number_of_particles]] + ['q%i'%i for i in [1..number_of_particles]] # initial values variable_names+=['p0%i'%i for i in [1..number_of_particles]] + ['q0%i'%i for i in [1..number_of_particles]] # coupling constant and time variable_names+=['g','t'] # we can either do our calculations in the symbolic ring (so that tVhe equations look nicer) # or in a polynomial ring (which is much faster) if not use_symbolic_ring: QQi. = QuadraticField(-1) RR = PolynomialRing(QQi, names=variable_names) RR.inject_variables() p = RR.gens()[0:number_of_particles] q = RR.gens()[number_of_particles:2*number_of_particles] p0 = RR.gens()[2*number_of_particles:3*number_of_particles] q0 = RR.gens()[3*number_of_particles:4*number_of_particles] g,t = RR.gens()[4*number_of_particles:4*number_of_particles+2] else: QQi = QQ variable_names = join(variable_names,', ') p = var(variable_names)[0:number_of_particles] q = var(variable_names)[number_of_particles:2*number_of_particles] p0 = var(variable_names)[2*number_of_particles:3*number_of_particles] q0 = var(variable_names)[3*number_of_particles:4*number_of_particles] g,t =var(variable_names)[4*number_of_particles:4*number_of_particles+2] return (p,q,p0,q0,g,t) def bracket(a,b): return a*b-b*a def Poisson(f,g,p,q): return sum([diff(f,p[k])*diff(g,q[k]) - diff(g,p[k])*diff(f,q[k]) for k in range(len(p))]) def dot(a,b): return sum(x*y for (x,y) in zip(a,b)) def lax_pair(root_system, p, q, g): """ Calculate the root type lax pair for the C-M system associated to `root_system` USAGE: L,M,HH = lax_pair(root_system, p, q, g) The return values `L` and `M` form the Lax pair. The return value `HH` is the Hamiltonian. All of these are functions of the variables in the lists `p` and `q` and of the coupling constant `g`. """ Phi = root_system.ambient_space().roots() # define matrices def E(eta): return matrix(ZZ, [[int((alpha - beta)==eta) for beta in Phi] for alpha in Phi]) def x(t): return 1/t def y(t): return -1/t^2 def z(t): return -1/t^2 H = diagonal_matrix([dot(alpha.to_vector(), p) for alpha in Phi]) D = diagonal_matrix([z(dot(alpha.to_vector(),q)) + sum([z(dot(beta.to_vector(),q)) for beta in Phi if dot(beta.to_vector(),alpha.to_vector())==1]) for alpha in Phi]) X = sum([x(dot(eta.to_vector(),q))*(E(eta)+2*E(2*eta)) for eta in Phi]) Y = sum([y(dot(eta.to_vector(),q))*(E(eta)+E(2*eta)) for eta in Phi]) L = H + I*g*X M = I*g*(Y - D) HH = 1/2*dot(p,p) - g^2/2*sum([z(dot(eta.to_vector(),q)) for eta in Phi]) return (L,M,HH) def multisets_equal_approx(a,b, epsilon=10^(-10)): """ Compare two lists of floating-points values and returns `True` if both lists have the same elements with the same multiplicities The parameter `epsilon` controls how much the floating point values may differ to be considered equal. """ def multisets_equal_destructive(a,b): if a == [] and b == []: return True if a == []: return False x = a[0] minimal_diff = min(abs(y - x) for y in b) if minimal_diff >= epsilon: return False for y in b: if abs(y - x) <= minimal_diff: b.remove(y) a = a[1:] return multisets_equal_destructive(a,b) return False return multisets_equal_destructive(a[:], b[:]) def partial_permutations(l, length): """ Return all permutations of all choices of `length` elements from `l`. If `l` has double values, we return each distinct partial permutation only once. """ if l == [] or length == 0: yield [] return for x in set(l): m = l[:] m.remove(x) for p in partial_permutations(m,length-1): yield [x] + p def get_simple_root_weights(root_system): """ Return the list of coefficients of the simple roots in the sum of all roots. """ roots = root_system.root_lattice().roots() indices = root_system.root_lattice().index_set() pos_roots = [r for r in roots if not any(c < 0 for c in r.coefficients())] weights = {} for ix in indices: weights[ix] = 0 for root in pos_roots: for ix in indices: weights[ix] += root.coefficient(ix) return weights.values() def maximal_simple_root_value(root_system, possible_values): """ If `possible_values` are arranged additively over the positive roots in `root_system`, we can compute an upper bound on how great values for the simple roots can be. `maximal_simple_root_value` returns this upper bound. """ w = get_simple_root_weights(root_system) w.sort() p = sorted(possible_values) min_weight = w[0] w = w[1:] w.reverse() min_total_rest = sum(a*b for (a,b) in zip(w,p)) return (sum(possible_values) - min_total_rest)/min_weight def get_possible_values(root_system, simple_values): """ For a root system with `len(simple_values)` simple roots, and assuming that the values of the simple roots are given by `simple_values`, returns all values. These values are defined by additivity. """ roots = root_system.root_lattice().roots() indices = root_system.root_lattice().index_set() pos_roots = [r for r in roots if not any(c < 0 for c in r.coefficients())] possible_values = [] v = dict(zip(indices,simple_values)) for root in pos_roots: possible_values.append(sum(root.coefficient(ix) * v[ix] for ix in indices)) return possible_values def find_alpha_q_from_q(root_system, q_value): """ Find the values of the inner products (alpha, q) """ Phi = root_system.ambient_space().roots() return [dot(alpha.to_vector(), q_value) for alpha in Phi] def find_q_from_alpha_q(root_system, alpha_q, q, epsilon = 10^(-10), debug=False): """ Recover q from the values of the inner products (alpha, q) USAGE: possible_qs = find_q_from_alpha_q(root_system, alpha_q, q, epsilon) Here, `possible_qs` is a list of all possibilities assuming that `q` is in the Weyl chamber corresponding to the base choice of `root_system.ambient_space().simple_roots()` We have proved that generically, we find as many solutions as there are Dynkin diagram automorphisms, but we have also given an example (examples 2,3) where we find more. `epsilon` is a parameter controlling how much floating point values may differ to be considered equal. This is an implementation detail; we have tried to rely on it as litle as possible. """ rank = len(root_system.index_set()) alpha_q = [x for x in alpha_q if x>0] if(len(alpha_q)!=len(root_system.ambient_space().positive_roots())): raise ValueError("Error: number of positive values is unexpected") weights = get_simple_root_weights(root_system) max_simple_value = maximal_simple_root_value(root_system, alpha_q) total_of_values = sum(v for v in alpha_q) possible_simple_values = [v for v in alpha_q if v <= max_simple_value + epsilon] if(len(possible_simple_values) < rank): raise ValueError("Error: something went wrong while trying to calculate simple values") possible_matchings = [] delta_perm = [] for perm in partial_permutations(possible_simple_values, rank): delta = abs(sum(v*w for (v,w) in zip(perm,weights)) - total_of_values) delta_perm.append((delta,perm)) if delta >= epsilon: continue if(debug): print("found partial permutation:") print(perm) print("with total weight:") print(sum(v*w for (v,w) in zip(perm,weights))) print("should be") print(total_of_values) print("and with possible values:") print(sorted(get_possible_values(root_system, perm))) print("should be") print(sorted(alpha_q)) if(multisets_equal_approx(get_possible_values(root_system, perm), alpha_q, epsilon)): if(debug): print("so these are equal") possible_matchings.append(dict(zip(root_system.ambient_space().simple_roots(),perm))) else: if(debug): print("so these are unequal (epsilon is: %f" % epsilon) if(debug): min_delta_perm = sorted(delta_perm, key=lambda a: a[0])[0:2] print("Minimal delta values are") print([d for (d,p) in min_delta_perm]) print("occurring with possible_values:") print(sorted(get_possible_values(root_system,min_delta_perm[0][1]))) print("resp") print(sorted(get_possible_values(root_system,min_delta_perm[1][1]))) print("should be:") print(sorted(alpha_q)) retval = [] for m in possible_matchings: eqs_for_q = [SR(dot(alpha.to_vector(),q)) == m[alpha] for alpha in root_system.ambient_space().simple_roots()] if(root_system.cartan_type()[0] == 'A'): eqs_for_q += [SR(sum([q_i for q_i in q])) == 0] solutions_for_q = solve(eqs_for_q,*(map(SR,q)),solution_dict=True)[0] retval.append([SR(q_i).subs(solutions_for_q) for q_i in q]) return retval def select_q_in_domain(possible_qs, chi, sign_of_chi_q): """ Choose the value of q that is in the right fundamental domain The fundamental domains are distinguished by the sign of (chi, q). Generically, there is only one solution in a fundamental domain, but we have shown an example (2,3) where there is more than one. """ return [q for q in possible_qs if sign_of_chi_q * dot(q,chi)>=0] def lax_W(root_system,L,M,p,p0_value,q,q0_value,g,g_value,t): """ Return the `W` matrix associated to a Lax pair USAGE: W = lax_W(root_system,L,M,p,p0_value,q,q0_value,g,g_value) The `W` matrix is obtained by substituting values `p0,q0,g_value` for `p,q,g` into the matrix diag(q) + t*L """ Phi = root_system.ambient_space().roots() ev0 = [dot(alpha.to_vector(),q0_value) for alpha in Phi] Q0 = diagonal_matrix(ev0) # substitute q0 and p0 in L L0 = L for (a,b) in zip(p,p0_value): L0 = L0.subs({a:b}) for (a,b) in zip(q,q0_value): L0 = L0.subs({a:b}) L0=L0.subs({g:g_value}) W = Q0 + t*L0 return W def find_unique_q(q_values, epsilon): """ Because of numerical precision issues, we may find the same solution for q a number of times, but with small differences. We try to identify those. """ unique_q_values = [] for q in q_values: found_one = False for existing_q in unique_q_values: if max(abs(v-w) for (v,w) in zip(q,existing_q)) < epsilon: found_one = True if not found_one: unique_q_values.append(q) return unique_q_values def calculate_q_at_time(root_system, W, time, q, chi, sign_of_chi_q, t, epsilon = 10^(-10)): W_t = W.subs(t = time).change_ring(SR).change_ring(CC) # get the eigenvalues of W_t eigenvalues = list(numpy.linalg.eigh(numpy.array(W_t))[0]) q_t = find_q_from_alpha_q(root_system,eigenvalues,q, epsilon) if(len(q_t)==0): print("warning: no solutions for q when alpha_q is") print(sorted(eigenvalues)) print("now running with debug flag") q_t = find_q_from_alpha_q(root_system,eigenvalues,q,epsilon,debug=True) q_t = find_unique_q(q_t, epsilon) q_t = select_q_in_domain(q_t, chi, sign_of_chi_q) return q_t def num_approx(q,p,q0_value,p0_value,g,HH,times): dHdq = [diff(HH,q_i) for q_i in q] dHdp = [diff(HH,p_i) for p_i in p] q_t=q0 p_t=p0 q_t = [q0] p_t = [p0] last_time = times[0] for t in times[1:]: d = dict([(q[i], q_t[n][i]) for i in range(len(q))] + [(p[i],p_t[n][i]) for i in range(len(p))] ) q_t[n+1] = [CC(q_t[n][i] + delta_t * dHdp[i].subs(d)) for i in range(len(q))] p_t[n+1] = [CC(p_t[n][i] - delta_t * dHdq[i].subs(d)) for i in range(len(p))] t = t +delta_t return q_t,p_t def get_indicator_discriminant(root_system,p0_value,q0_value,g_value,chi): """ Calculate the discrimant, a factor of which is the indicator function we are looking for. USAGE: d = get_indicator_discriminant(root_system,p0_value,q0_value,g_value,chi) or d = get_indicator_discriminant(*examples[i]) """ p,q,p0,q0,g,t = define_variables(3, False) L,M,HH = lax_pair(root_system,p,q,g) W = lax_W(root_system,L,M,p,p0_value,q,q0_value,g,g_value,t) charpoly = W.characteristic_polynomial() x = charpoly.variables()[0] # coerce to make the base ring a fraction field, otherwise Sage # cannot calculate the discriminant charpoly = (FractionField(QQ[t])[x])(charpoly) return charpoly.discriminant() def some_values(length,possible_values): if length <= 0: yield [] return for vs in some_values(length - 1, possible_values): for v in possible_values: yield [v] + vs def multisets_equal_hashable(a,b): """ Compares two lists and returns `True` if both lists have the same elements with the same multiplicities. `multisets_equal_hashable` differs from `multisets_equal` in that it assumes that the list elements are hashable. """ if len(a) != len(b): return False a_multiplicities = {} b_multiplicities = {} for x in a: if x in a_multiplicities.keys(): a_multiplicities[x] += 1 else: a_multiplicities[x] = 1 for x in b: if x in b_multiplicities.keys(): b_multiplicities[x] += 1 else: b_multiplicities[x] = 1 return a_multiplicities == b_multiplicities def find_problems(root_system,min_value,max_value): """ Try to find an additive assignment of positive values to the positive roots in `root_system`, which have more than the expected number (2 in most cases) of ways of assigning these numbers. """ rank = len(root_system.index_set()) weights = get_simple_root_weights(root_system) expected_number = root_system.dynkin_diagram().automorphism_group().order() for simple_values in some_values(rank, range(1,max_value)): if (gcd(simple_values) > 1): continue if (max(simple_values) < min_value): continue possible_values = get_possible_values(root_system, simple_values) total_of_values = sum(possible_values) max_simple_value = maximal_simple_root_value(root_system, possible_values) possible_simple_values = [v for v in possible_values if v <= max_simple_value] possible_permutations = [] for a in partial_permutations(possible_simple_values,rank): if sum(v*w for (v,w) in zip(a,weights)) != total_of_values: continue new_possible_values = get_possible_values(root_system, a) if multisets_equal_hashable(possible_values,new_possible_values): possible_permutations.append(a) print possible_permutations if len(possible_permutations) > expected_number: return possible_permutations