|\^/| Maple 15 (IBM INTEL LINUX) ._|\| |/|_. Copyright (c) Maplesoft, a division of Waterloo Maple Inc. 2011 \ MAPLE / All rights reserved. Maple is a trademark of <____ ____> Waterloo Maple Inc. | Type ? for help. # ======================================================================== # # # Copyright 2012, Gabor Hegedus, Josef Schicho, and Hans-Peter Schroecker # # This program 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. # # This program 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 . # # ======================================================================== # > > read("bonds.mpl"): > # define planar four-bar linkage # ============================== > h1 := [0, e, 0, 1]: > h2 := [0, 0, e, 1]: > h3 := [0, 0, 0, 1]: > h4 := [0, e, 2*e, 1]: > # sanity checks of input data > for i from 1 to 4 do > if Mul(h||i, Cj(h||i)) <> [1, 0, 0, 0] then > print("Error: Input quaternions are not normalized!"): > stop: > end if: > if h||i + Cj(h||i) <> [0, 0, 0, 0] then > print("Error: Input quaternions are not rotations quaternions!"): > stop: > end if: > od: > # Ideal of configuration curve # ============================ > > for i from 1 to 4 do: > u||i := [t||i, 0, 0, 0] - h||i: > od: > eqs := RVec(Mulall(u1, u2, u3, u4)): > > es := {op(eqs[2..-1]), eqs[1]*u-1}: > vars := [t1, t2, u, t3, t4]: > Gb1 := Groebner[Basis]([op(es)], tdeg(t1, t2, u, t3, t4)): > egb := select(T->not (u in indets(T)), Gb1): > egb := Groebner[Basis]({op(egb)}, plex(t4, t3, t2, t1)): > configcurve := allvalues(solve(egb, [t1, t2, t3, t4])): > simplify(configcurve[1]); 2 2 2 t4 - 2 t4 + %1 - 5 t4 + 1 + %1 t4 + 1 - 4 t4 + %1 [[t1 = - -------------------, t2 = - ------------, t3 = -------------------, 2 (t4 + 3) 4 (t4 - 2) 4 (-1 + t4) t4 = t4]] 4 2 3 1/2 %1 := (t4 + 2 t4 - 8 t4 - 47 + 56 t4) > simplify(configcurve[2]); 2 2 2 -t4 + 2 t4 + %1 + 5 -t4 - 1 + %1 -t4 - 1 + 4 t4 + %1 [[t1 = --------------------, t2 = -------------, t3 = - --------------------, 2 (t4 + 3) 4 (t4 - 2) 4 (-1 + t4) t4 = t4]] 4 2 3 1/2 %1 := (t4 + 2 t4 - 8 t4 - 47 + 56 t4) > # parametrized equation of (parts of) configuration curve > C1 := map(simplify, subs(t4 = t, [rhs(configcurve[1][1][1]), > rhs(configcurve[1][1][2]), > rhs(configcurve[1][1][3]), > rhs(configcurve[1][1][4])])); 2 2 2 t - 2 t + %1 - 5 t + 1 + %1 t + 1 - 4 t + %1 C1 := [- -----------------, - -----------, -----------------, t] 2 (t + 3) 4 (t - 2) 4 (-1 + t) 4 2 3 1/2 %1 := (t + 2 t - 8 t - 47 + 56 t) > C2 := map(simplify, subs(t4 = t, [rhs(configcurve[2][1][1]), > rhs(configcurve[2][1][2]), > rhs(configcurve[2][1][3]), > rhs(configcurve[2][1][4])])); 2 2 2 -t + 2 t + %1 + 5 -t - 1 + %1 -t - 1 + 4 t + %1 C2 := [------------------, ------------, - ------------------, t] 2 (t + 3) 4 (t - 2) 4 (-1 + t) 4 2 3 1/2 %1 := (t + 2 t - 8 t - 47 + 56 t) > vars:=[t1, t2, t3, t4]: > eqss := [op(egb), eqs[1]]: > bonds := solve(eqss, vars); bonds := [[t1 = -%2, t2 = %2, t3 = -%2, t4 = %2], [t1 = 2 + %2, t2 = -%2, t3 = -1 + 2 %2, t4 = %2], [t1 = 4 - %1, t2 = 2 - %1, t3 = -4 + %1, t4 = %1]] 2 %1 := RootOf(_Z - 8 _Z + 17) 2 %2 := RootOf(_Z + 1) > # The next step looses conjugate complex bonds. This is no problem as # the bond structure of conjugate complex bonds is equal. > bonds := simplify(convert(bonds, radical)): > > for b in bonds do > print(b): > od: [t1 = -I, t2 = I, t3 = -I, t4 = I] [t1 = 2 + I, t2 = -I, t3 = -1 + 2 I, t4 = I] [t1 = -I, t2 = -2 - I, t3 = I, t4 = 4 + I] > # Verify vanishing of certain coupler polynomials: > u5 := u1: u6 := u2: u7 := u3: u8 := u4: > > for j from 1 to nops(bonds) do > assign(bonds[j]): > printf(" beta = (%a, %a, %a, %a)\n\n", t1, t2, t3, t4): > for i from 1 to 4 do > for k from i+1 to i+3 do > F||i||k := Mulall(seq(u||l, l=i+1..k)): > if F||i||k = [0, 0, 0, 0] then > printf(" F_{%d,%d}(beta) = %a\n\n", i, cc(k, 4), F||i||k): > end if: > od: > od: > unassign('t1', 't2', 't3', 't4'): > od: beta = (-I, I, -I, I) F_{1,4}(beta) = [0, 0, 0, 0] F_{2,1}(beta) = [0, 0, 0, 0] F_{3,2}(beta) = [0, 0, 0, 0] F_{4,3}(beta) = [0, 0, 0, 0] beta = (2+I, -I, -1+2*I, I) F_{1,4}(beta) = [0, 0, 0, 0] F_{3,2}(beta) = [0, 0, 0, 0] beta = (-I, -2-I, I, 4+I) F_{2,1}(beta) = [0, 0, 0, 0] F_{4,3}(beta) = [0, 0, 0, 0] > > C1 := subs(t4 = t, [rhs(configcurve[1][1][1]), > rhs(configcurve[1][1][2]), > rhs(configcurve[1][1][3]), > rhs(configcurve[1][1][4])]): > C2 := subs(t4 = t, [rhs(configcurve[2][1][1]), > rhs(configcurve[2][1][2]), > rhs(configcurve[2][1][3]), > rhs(configcurve[2][1][4])]): > C1 := map(simplify, C1); 2 2 2 t - 2 t + %1 - 5 t + 1 + %1 t + 1 - 4 t + %1 C1 := [- -----------------, - -----------, -----------------, t] 2 (t + 3) 4 (t - 2) 4 (-1 + t) 4 2 3 1/2 %1 := (t + 2 t - 8 t - 47 + 56 t) > C2 := map(simplify, C2); 2 2 2 -t + 2 t + %1 + 5 -t - 1 + %1 -t - 1 + 4 t + %1 C2 := [------------------, ------------, - ------------------, t] 2 (t + 3) 4 (t - 2) 4 (-1 + t) 4 2 3 1/2 %1 := (t + 2 t - 8 t - 47 + 56 t) > # Compute all local distances # =========================== > for j from 1 to nops(bonds) do > assign(bonds[j]): > t0 := t4: > C := subs(t = t0, C1): > if (map(simplify, C - [t1, t2, t3, t4]) <> [0, 0, 0, 0]) then > C := C2: > else > C := C1: > end if: > printf("\n %a\n", map(simplify, subs(t = t0, C))): > for k from 1 to 4 do > for l from k+1 to 4 do > f||k||l := Mulall(seq([C[m], 0, 0, 0] - h||m, m=k+1..l)): > v||k||l := min(map(vanishingOrder, f||k||l, t, t0)): > d||k||l := 1/2*vanishingOrder(pNorm(f||k||l), t, t0) - v||k||l: > printf(" d(%d, %d) = %a (v_{%d,%d}=%a)\n", k, l, d||k||l, k, l, v||k||l): > od: > od: > unassign('t0', 't1', 't2', 't3', 't4'): > od: [-I, I, -I, I] d(1, 2) = 1/2 (v_{1,2}=0) d(1, 3) = 1 (v_{1,3}=0) d(1, 4) = 1/2 (v_{1,4}=1) d(2, 3) = 1/2 (v_{2,3}=0) d(2, 4) = 1 (v_{2,4}=0) d(3, 4) = 1/2 (v_{3,4}=0) [2+I, -I, -1+2*I, I] d(1, 2) = 1/2 (v_{1,2}=0) d(1, 3) = 1/2 (v_{1,3}=0) d(1, 4) = 0 (v_{1,4}=1) d(2, 3) = 0 (v_{2,3}=0) d(2, 4) = 1/2 (v_{2,4}=0) d(3, 4) = 1/2 (v_{3,4}=0) [-I, -2-I, I, 4+I] d(1, 2) = 0 (v_{1,2}=0) d(1, 3) = 1/2 (v_{1,3}=0) d(1, 4) = 1/2 (v_{1,4}=0) d(2, 3) = 1/2 (v_{2,3}=0) d(2, 4) = 1/2 (v_{2,4}=0) d(3, 4) = 0 (v_{3,4}=0) > # Details of computation of d14 and bond (-I, I, -I, I) # ===================================================== > t0 := I: > tmp := simplify(subs(t = t0, C1)): > if (tmp[1]^2 + 1 = 0) and (tmp[2]^2 + 1 = 0) and (tmp[2]^2 + 1 = 0) and (tmp[3]^2 + 1 = 0) then > C := C1: > else > C := C2: > end if: > C0 := simplify(subs(t = t0, C)); C0 := [-I, I, -I, I] > for i from 1 to 4 do > t||i := subs(t = t0, C0[i]): > od: > > F14 := map(factor, Mulall([C[2], 0, 0, 0] - h2, [C[3], 0, 0, 0] - h3, [C[4], 0, 0, 0] - h4)): > v14 := min(map(vanishingOrder, F14, t, t0)); v14 := 1 > d14 := 1/2*vanishingOrder(simplify(pNorm(F14)), t, t0); d14 := 3/2 > d14 := d14 - v14; d14 := 1/2 > quit memory used=25.6MB, alloc=21.1MB, time=0.61