with(linalg): alias( om = RootOf(Z_^2 + Z_ + 1 = 0) ); eqL[1] := (y = z): eqL[4] := (z = x): eqL[7] := (x = y): eqL[2] := (y = om*z): eqL[5] := (z = om*x): eqL[8] := (x = om*y): eqL[3] := (y = om^2*z): eqL[6] := (z = om^2*x): eqL[9] := (x = om^2*y): p1 := [1,1,1]: p4 := [om^2,1,1]: p7 := [1,om,1]: p2 := [1,om^2,om]: p5 := [1,1,om^2]: p8 := [om,1,1]: p3 := [1,om,om^2]: p6 := [1,om^2,1]: p9 := [1,1,om]: ii1 := [1,4,7]: ii4 := [1,5,9]: ii7 := [2,4,9]: ii2 := [2,5,8]: ii5 := [2,6,7]: ii8 := [1,6,8]: ii3 := [3,6,9]: ii6 := [3,4,8]: ii9 := [3,5,7]: # # ii1[j],...,ii9[j] are denoted in the paper by $i_j(1),...,i_j(9)$ # pp := [p1,p2,p3,p4,p5,p6,p7,p8,p9]: ii := [ii1,ii2,ii3,ii4,ii5,ii6,ii7,ii8,ii9]: ## ii[k,j] is $i_j(k)$ p := seq( [pp[k,1]/pp[k,3],pp[k,2]/pp[k,3]], k=1..9 ); ## p[k] is $p_k$ is the affine coordinates # # Check that the numbering of the triple points is correct: # factor(subs(solve([x = 1, eqL[1], eqL[4], eqL[7]]), [x,y,z]) - p1), factor(subs(solve([x = 1, eqL[2], eqL[5], eqL[8]]), [x,y,z]) - p2), factor(subs(solve([x = 1, eqL[3], eqL[6], eqL[9]]), [x,y,z]) - p3), factor(subs(solve([z = 1, eqL[1], eqL[5], eqL[9]]), [x,y,z]) - p4), factor(subs(solve([y = 1, eqL[2], eqL[6], eqL[7]]), [x,y,z]) - p5), factor(subs(solve([x = 1, eqL[3], eqL[4], eqL[8]]), [x,y,z]) - p6), factor(subs(solve([z = 1, eqL[2], eqL[4], eqL[9]]), [x,y,z]) - p7), factor(subs(solve([y = 1, eqL[1], eqL[6], eqL[8]]), [x,y,z]) - p8), factor(subs(solve([x = 1, eqL[3], eqL[5], eqL[7]]), [x,y,z]) - p9); # # Check that the numbers ii[k,j] are correct: # for k from 1 to 9 do solve([eqL[ii[k,1]],eqL[ii[k,2]],eqL[ii[k,3]]]); print(factor(subs(%,[x/z,y/z])-p[k])) od: L0x:=T: L0y:=1-T: L0z:=2: parL0 := (x=L0x, y=L0y, z=L0z): ## parametrization of the line $L_0$ ## ## Check that the ponts $p_i$ do not lye on the line $L_0$. ## seq([solve([p[i,1]=L0x/L0z,p[i,2]=L0y/L0z],t)],i=1..9); ## should be [],[],...,[] for i from 1 to 9 do t[i] := simplify(solve(subs(parL0,eqL[i]),T)); od; t1:=t[1]: t2:=t[2]: t3:=t[3]: t4:=t[4]: t5:=t[5]: t6:=t[6]: t7:=t[7]: t8:=t[8]: t9:=t[9]: a[1]:=a1: a[2]:=a2: a[3]:=a3: a[4]:=a4: a[5]:=a5: a[6]:=a6: a[7]:=a7: a[8]:=a8: a[9]:=a9: b[1]:=b1: b[2]:=b2: b[3]:=b3: b[4]:=b4: b[5]:=b5: b[6]:=b6: b[7]:=b7: b[8]:=b8: b[9]:=b9: xx := L0x *(T-t1-u*(a1+b1))*(T-t2-u*(a2+b2))*(T-t3-u*(a3+b3)) *(T-t4-u*a4)*(T-t5-u*a5)*(T-t6-u*a6) *(T-t7-u*a7)*(T-t8-u*a8)*(T-t9-u*a9): yy := L0y *(T-t1-u*a1)*(T-t2-u*a2)*(T-t3-u*a3) *(T-t4-u*(a4+b4))*(T-t5-u*(a5+b5))*(T-t6-u*(a6+b6)) *(T-t7-u*a7)*(T-t8-u*a8)*(T-t9-u*a9): zz := L0z *(T-t1-u*a1)*(T-t2-u*a2)*(T-t3-u*a3) *(T-t4-u*a4)*(T-t5-u*a5)*(T-t6-u*a6) *(T-t7-u*(a7+b7))*(T-t8-u*(a8+b8))*(T-t9-u*(a9+b9)): X := factor(xx/zz); Y := factor(yy/zz); ## Parametrization of the curve C in the affine coordinates xx := 'xx': yy := 'yy': zz := 'zz': ## Forget about the projective parametrization of C #================================================================================================= # # \S 2.3. Perturbation with triple points preserved up to O(u^2) # #================================================================================================= for k from 1 to 9 do # compute tau[k,j] for j from 1 to 3 do i:=ii[k,j]; so := solve( {subs(u=0,factor(subs( T = t[i] + Tau*u, X ))) = p[k,1], subs(u=0,factor(subs( T = t[i] + Tau*u, Y ))) = p[k,2]}, Tau ); tau[k,j] := factor(subs(so,Tau)); od; print( k,tau[k,1],tau[k,2],tau[k,3] ) od: # check the expression for tau[1,2] in the paper: "Check tau[1,2]: ", [[simplify( a4 + 1/3*b4 - tau[1,2] )]]; #================================================================================================= w1 := [1,0]: w2 := [0,1]: for k from 1 to 9 do for j from 1 to 3 do i:=ii[k,j]; print([k,i]); se := series(subs(T=t[i]+tau[k,j]*u,X),u,3); vX[j] := factor(op(3,se)); # the coefficient of u in X(...) se := series(subs(T=t[i]+tau[k,j]*u,Y),u,3); vY[j] := factor(op(3,se)); # the coefficient of u in Y(...) print(vX[j]); print(vY[j]); # # check the formula for the X-component of $p'_{1,2}(0)$ in the paper: # if k=1 and j=2 then print("check p'_{1,2}(0): ", [[[simplify( -2*b1 + (2+4*om)*(b2-b3) + 3*a4 + b4 + 4*b7 + (4+2*om)*b8 + (2-2*om)*b9 - 6*vX[j] )]]] ) fi; od; v1 := [vX[1],vY[1]]; v2 := [vX[2],vY[2]]; v3 := [vX[3],vY[3]]; w3 := p[k]; M:=matrix(4,4, [ w1[1], -w2[1], 0, v2[1]-v1[1], w1[2], -w2[2], 0, v2[2]-v1[2], w1[1], 0, -w3[1], v3[1]-v1[1], w1[2], 0, -w3[2], v3[2]-v1[2] ] ): ## the matrix in Equation (3) e[k] := expand(simplify(factor(det(M)))); print(e[k]); od: # # check the formula for $e_1$ in the paper: # "check e_1: ",simplify( 1/2*(a1+b1+a4+b4) - 4*(a7+b7) + b8 + b9 + (2*om-1)/7*(b2+b6) + (2*om^2-1)/7*(b3+b5) - e[1] ); # ====================================================================================== so := solve({a9=0,a4=a1,b8=1,b9=1,b1=1,b2=1,b3=1,b5=1,b6=1,seq(e[k],k=1..9)}); ## ## subs(so,a[i]) is the value of the i-th component of $\hat\aa$ (see \S 2.3), ## which is also $a_{i,0}$ in the notation of \S 2.5 ## # # Check the expression for $\hat\aa$ in the paper (Equation (4)): # Alpha := (38*om - 5)/21: ConjAlpha := (38*om^2 - 5)/21: simplify( [-1/3,Alpha,ConjAlpha, -1/3,ConjAlpha,Alpha, 37/42,0,0] - subs(so,[a1,a2,a3,a4,a5,a6,a7,a8,a9]) ); simplify( [1,1,1, 1,1,1, -1/2,1,1] - subs(so,[b1,b2,b3,b4,b5,b6,b7,b8,b9]) ); # ====================================================================================== # # \S 2.5; Initialization of the recursion # # ====================================================================================== X0 := factor(subs(so,X)); Y0 := factor(subs(so,Y)); for k from 1 to 9 do for j from 1 to 3 do i:=ii[k,j]; tau0[k,j] := simplify(subs(so,tau[k,j])) od: od: # # [X0,Y0] is $f_{u,\hat\aa}$ (here and below \aa is the boldface a ) # tau0[k,j] is $\tau_{k,j}(\hat\aa)$ # # ====================================================================================== # # compute S0[k,j], which is $S_{k,j}^{[0]}$ # for k from 1 to 9 do i1 := ii[k,1]; i2 := ii[k,2]; i3 := ii[k,3]; print(k,ii[k]); se := series( subs( T = t[i1] + tau0[k,1]*u + s1*u^2, X0 ),u,3): xx1 := factor(op(3,se)); # the coef of u^1 se := series( subs( T = t[i2] + tau0[k,2]*u + s2*u^2, X0 ),u,3): xx2 := factor(op(3,se)); se := series( subs( T = t[i3] + tau0[k,3]*u + s3*u^2, X0 ),u,3): xx3 := factor(op(3,se)); se := series( subs( T = t[i1] + tau0[k,1]*u + s1*u^2, Y0 ),u,3): yy1 := factor(op(3,se)); se := series( subs( T = t[i2] + tau0[k,2]*u + s2*u^2, Y0 ),u,3): yy2 := factor(op(3,se)); se := series( subs( T = t[i3] + tau0[k,3]*u + s3*u^2, Y0 ),u,3): yy3 := factor(op(3,se)); so123 := solve( {xx1=xx3, xx2=xx3, yy1=yy3, yy2=yy3} ); print([xx1,yy1],[xx2,yy2],[xx3,yy3]); S0[k,1] := subs( so123, s1 ); S0[k,2] := subs( so123, s2 ); S0[k,3] := subs( so123, s3 ); print( S0[k,1],S0[k,2],S0[k,3] ); od: # ====================================================================================== # # The recursion step [m-1] ---> [m] for m = 1 # # ====================================================================================== su := seq( a[k]=subs(so,a[k])+u*a[k], k=1..9 ), seq( b[k]=subs(so,b[k])+u*b[k], k=1..9 ); # # subs(su,a[k]) is the k-th component of $\aa^{[1]}$ where # the coefficient of $u^0$ is already computed and the coefficient # of $u^1$ is yet indeterminate. The same for subs(su,b[k]). # # Since $a_{i,n}$ with different $n$ are never used simultaneously, # they are denoted in the maple code by a[i] # X1a := factor(subs(su,X)); Y1a := factor(subs(su,Y)); for k from 1 to 9 do for j from 1 to 3 do tau1a[k,j] := simplify(subs(su,tau[k,j])) od: od: # # [X1a,Y1a] and tau1a[k,j] are $f_{u,a^{[1]}(u)}$ and $\tau(a^{[1]}$ respectively # with indeterminate coefficients $a_{k,1}$ and $b_{k,1}$ # # ====================================================================================== # Compute $s_{k,j}$, $s_{k,j}^*$, and Equation (7). for k from 1 to 9 do i1 := ii[k,1]; i2 := ii[k,2]; i3 := ii[k,3]; se := series( subs( T = t[i1] + tau1a[k,1]*u + (S0[k,1] + s1*u)*u^2, X1a ),u,4): xx1 := simplify(expand(op(5,se))); # coef of u^2 se := series( subs( T = t[i2] + tau1a[k,2]*u + (S0[k,2] + s2*u)*u^2, X1a ),u,4): xx2 := simplify(expand(op(5,se))); se := series( subs( T = t[i3] + tau1a[k,3]*u + (S0[k,3] + s3*u)*u^2, X1a ),u,4): xx3 := simplify(expand(op(5,se))); se := series( subs( T = t[i1] + tau1a[k,1]*u + (S0[k,1] + s1*u)*u^2, Y1a ),u,4): yy1 := simplify(expand(op(5,se))); se := series( subs( T = t[i2] + tau1a[k,2]*u + (S0[k,2] + s2*u)*u^2, Y1a ),u,4): yy2 := simplify(expand(op(5,se))); se := series( subs( T = t[i3] + tau1a[k,3]*u + (S0[k,3] + s3*u)*u^2, Y1a ),u,4): yy3 := simplify(expand(op(5,se))); so13 := solve({xx1=xx3,yy1=yy3},{s1,s3}); so23 := solve({xx2=xx3,yy2=yy3},{s2,s3}); s1star[k] := factor(subs(so13,s1)); ### $s_{k,1}^*$ s2star[k] := factor(subs(so23,s2)); ### $s_{k,2}^*$ s1a[k] := factor(subs(so13,s3)); ### $s_{k,1}$ s2a[k] := factor(subs(so23,s3)); ### $s_{k,2}$ e1[k] := factor(s1a[k]-s2a[k]); ### the equation (7) print([k],expand(simplify(e1[k]))); od: # ====================================================================================== # Solve Equation (7): so1 := solve({a4=a1, a7=0, b1=0, seq(b[k]=0,k = 4..9), seq(e1[k]=0,k=1..9)}); # # subs(so1,a[i]) is $a_{i,1}$ # subs(so1,b[i]) is $b_{i,1}$ # var := [a1,a2,a3, a5,a6,a8, a9,b2,b3]; li:=[]: for k from 1 to 9 do li:=[op(li),seq(diff( subs(a4=a1,e1[k]), var[i] ),i=1..9)] od: M:=matrix(9,9,li); simplify( det(M) ); # -64 # -------- # 26040609 # # M is the matrix of $\Phi|_{0\times E}$ # # Check the value of det(M) in the paper: "Check det(Phi): ", simplify( -64/(3^12 * 49)/det(M) ); # ====================================================================================== # # [X1,Y1] and tau1[k,j] are $f_{u,a^{[1]}(u)}$ and $\tau(a^{[1]}$ respectively # with all coefficients $a_{i,n}$, $b_{i,n}$, $n=0,1$, computed # # S1[k,j] is $S_{k,j}^{[1]}$ # T1[k,j] is $T_{k,j}^{[1]}$ X1 := factor(subs(so1,X1a)); Y1 := factor(subs(so1,Y1a)); T1 := collect([seq([seq(0,j=1..3)], k=1..9 )],u): for k from 1 to 9 do for j from 1 to 3 do tau1[k,j] := simplify(subs(so1,tau1a[k,j])); od; S1[k,1] := S0[k,1] + simplify( subs(so1,s1star[k]) )*u; S1[k,2] := S0[k,2] + simplify( subs(so1,s2star[k]) )*u; S1[k,3] := S0[k,3] + simplify( subs(so1,s1a[k]) )*u; for j from 1 to 3 do i:=ii[k,j]; tau1[k,j] := simplify(subs(so1,tau1a[k,j])); T1[k,j] := collect( expand(simplify( t[i] + tau1[k,j]*u + S1[k,j]*u^2 )), u ); od; od: T1; X1o := subs(om=o,X1): Y1o := subs(om=o,Y1): T1o := subs(om=o,T1): save X1o,Y1o,T1o, "12tp-approx1.sav"; # ====================================================================================== # # Check that [X1,Y1] is a solution up to $O(u^3)$. # for k from 1 to 9 do for j from 1 to 3 do xx[j] := convert(series( subs( T = T1[k,j], X1 ), u,4),polynom); yy[j] := convert(series( subs( T = T1[k,j], Y1 ), u,4),polynom); od; print( factor(xx[1]-xx[2]), factor(xx[2]-xx[3]), factor(yy[1]-yy[2]), factor(yy[2]-yy[3]) ); od: # ====================================================================================== # # The recursion step [m-1] ---> [m] for m = 2 # # ====================================================================================== su1 := seq( a[k]=subs(so1,a[k])+u*a[k], k=1..9 ), seq( b[k]=subs(so1,b[k])+u*b[k], k=1..9 ); # # su1 replaces the indeterminate $a_{i,1}$ (which is represented by a[i] in the maple code) # by $a_{i,1} + a_{i,2} u$ with the computed value of $a_{i,1}$ and # indeterminate $a_{i,2}$, which is still represented by a[i] in the maple code # X2a := factor(subs(su1,X1a)); Y2a := factor(subs(su1,Y1a)); for k from 1 to 9 do for j from 1 to 3 do tau2a[k,j] := simplify(subs(su1,tau1a[k,j])) od: od: # # [X2a,Y2a] and tau2a[k,j] are $f_{u,a^{[2]}(u)}$ and $\tau(a^{[2]}$ respectively # with indeterminate coefficients $a_{k,2}$ and $b_{k,2}$ # # ====================================================================================== # Compute $s_{k,j}$, $s_{k,j}^*$, and Equation (7). for k from 1 to 9 do i1 := ii[k,1]; i2 := ii[k,2]; i3 := ii[k,3]; se := series( subs( T = t[i1] + tau2a[k,1]*u + (S1[k,1] + s1*u^2)*u^2, X2a ),u,5): xx1 := simplify(expand(op(7,se))); # coef of u^3 se := series( subs( T = t[i2] + tau2a[k,2]*u + (S1[k,2] + s2*u^2)*u^2, X2a ),u,5): xx2 := simplify(expand(op(7,se))); se := series( subs( T = t[i3] + tau2a[k,3]*u + (S1[k,3] + s3*u^2)*u^2, X2a ),u,5): xx3 := simplify(expand(op(7,se))); se := series( subs( T = t[i1] + tau2a[k,1]*u + (S1[k,1] + s1*u^2)*u^2, Y2a ),u,5): yy1 := simplify(expand(op(7,se))); se := series( subs( T = t[i2] + tau2a[k,2]*u + (S1[k,2] + s2*u^2)*u^2, Y2a ),u,5): yy2 := simplify(expand(op(7,se))); se := series( subs( T = t[i3] + tau2a[k,3]*u + (S1[k,3] + s3*u^2)*u^2, Y2a ),u,5): yy3 := simplify(expand(op(7,se))); so13 := solve({xx1=xx3,yy1=yy3},{s1,s3}); so23 := solve({xx2=xx3,yy2=yy3},{s2,s3}); s1star[k] := factor(subs(so13,s1)); ### $s_{k,1}^*$ s2star[k] := factor(subs(so23,s2)); ### $s_{k,2}^*$ s1a[k] := factor(subs(so13,s3)); ### $s_{k,1}$ s2a[k] := factor(subs(so23,s3)); ### $s_{k,2}$ e2[k] := factor(s1a[k]-s2a[k]); ### the equation (7) print([k],expand(simplify(e2[k]))); od: # ====================================================================================== # Solve Equation (7): so2 := solve({a4=a1, a7=0, b1=0, seq(b[k]=0,k = 4..9), seq(e2[k]=0,k=1..9)}); # # subs(so2,a[i]) is $a_{i,2}$ # subs(so2,b[i]) is $b_{i,2}$ # # ====================================================================================== # # [X2,Y2] and tau2[k,j] are $f_{u,a^{[2]}(u)}$ and $\tau(a^{[2]}$ respectively # with all coefficients $a_{i,n}$, $b_{i,n}$, $n=0,1,2$, computed # # S2[k,j] is $S_{k,j}^{[2]}$ # T2[k,j] is $T_{k,j}^{[2]}$ X2 := factor(subs(so2,X2a)); Y2 := factor(subs(so2,Y2a)); T2 := [seq([seq(0,j=1..3)], k=1..9 )]: for k from 1 to 9 do for j from 1 to 3 do tau2[k,j] := simplify(subs(so2,tau2a[k,j])); od; S2[k,1] := S1[k,1] + simplify( subs(so2,s1star[k]) )*u^2; S2[k,2] := S1[k,2] + simplify( subs(so2,s2star[k]) )*u^2; S2[k,3] := S1[k,3] + simplify( subs(so2,s2a[k]) )*u^2; for j from 1 to 3 do i:=ii[k,j]; tau2[k,j] := simplify(subs(so2,tau2a[k,j])); T2[k,j] := collect( expand(simplify( t[i] + tau2[k,j]*u + S2[k,j]*u^2 )), u ); od; od: T2; X2o := subs(om=o,X2): Y2o := subs(om=o,Y2): T2o := subs(om=o,T2): save X2o,Y2o,T2o, "12tp-approx2.sav"; # ====================================================================================== # # Check that [X2,Y2] is a solution up to $O(u^4)$. # for k from 1 to 9 do for j from 1 to 3 do xx[j] := convert(series( subs( T = T2[k,j], X2 ), u,5),polynom); yy[j] := convert(series( subs( T = T2[k,j], Y2 ), u,5),polynom); od; print( factor(xx[1]-xx[2]), factor(xx[2]-xx[3]), factor(yy[1]-yy[2]), factor(yy[2]-yy[3]) ); od: # ====================================================================================== # # The recursion step [m-1] ---> [m] for m = 3 # # ====================================================================================== su2 := seq( a[k]=subs(so2,a[k])+u*a[k], k=1..9 ), seq( b[k]=subs(so2,b[k])+u*b[k], k=1..9 ); # # su2 replaces the indeterminate $a_{i,2}$ (which is represented by a[i] in the maple code) # by $a_{i,2} + a_{i,3} u$ with the computed value of $a_{i,2}$ and # indeterminate $a_{i,3}$, which is still represented by a[i] in the maple code # X3a := factor(subs(su2,X2a)); Y3a := factor(subs(su2,Y2a)); for k from 1 to 9 do for j from 1 to 3 do tau3a[k,j] := simplify(subs(su2,tau2a[k,j])) od: od: # # [X3a,Y3a] and tau3a[k,j] are $f_{u,a^{[3]}(u)}$ and $\tau(a^{[3]}$ respectively # with indeterminate coefficients $a_{k,3}$ and $b_{k,3}$ # # ====================================================================================== # Compute $s_{k,j}$, $s_{k,j}^*$, and Equation (7). for k from 1 to 9 do i1 := ii[k,1]; i2 := ii[k,2]; i3 := ii[k,3]; se := series( subs( T = t[i1] + tau3a[k,1]*u + (S2[k,1] + s1*u^3)*u^2, X3a ),u,6): xx1 := simplify(expand(op(9,se))); # coef of u^4 se := series( subs( T = t[i2] + tau3a[k,2]*u + (S2[k,2] + s2*u^3)*u^2, X3a ),u,6): xx2 := simplify(expand(op(9,se))); se := series( subs( T = t[i3] + tau3a[k,3]*u + (S2[k,3] + s3*u^3)*u^2, X3a ),u,6): xx3 := simplify(expand(op(9,se))); se := series( subs( T = t[i1] + tau3a[k,1]*u + (S2[k,1] + s1*u^3)*u^2, Y3a ),u,6): yy1 := simplify(expand(op(9,se))); se := series( subs( T = t[i2] + tau3a[k,2]*u + (S2[k,2] + s2*u^3)*u^2, Y3a ),u,6): yy2 := simplify(expand(op(9,se))); se := series( subs( T = t[i3] + tau3a[k,3]*u + (S2[k,3] + s3*u^3)*u^2, Y3a ),u,6): yy3 := simplify(expand(op(9,se))); so13 := solve({xx1=xx3,yy1=yy3},{s1,s3}); so23 := solve({xx2=xx3,yy2=yy3},{s2,s3}); s1star[k] := factor(subs(so13,s1)); ### $s_{k,1}^*$ s2star[k] := factor(subs(so23,s2)); ### $s_{k,2}^*$ s1a[k] := factor(subs(so13,s3)); ### $s_{k,1}$ s2a[k] := factor(subs(so23,s3)); ### $s_{k,2}$ e3[k] := factor(s1a[k]-s2a[k]); ### the equation (7) print([k],expand(simplify(e3[k]))); od: # ====================================================================================== # Solve Equation (7): so3 := solve({a4=a1, a7=0, b1=0, seq(b[k]=0,k = 4..9), seq(e3[k]=0,k=1..9)}); # # subs(so3,a[i]) is $a_{i,3}$ # subs(so3,b[i]) is $b_{i,3}$ # # ====================================================================================== # # [X3,Y3] and tau3[k,j] are $f_{u,a^{[3]}(u)}$ and $\tau(a^{[3]}$ respectively # with all coefficients $a_{i,n}$, $b_{i,n}$, $n=0,1,2,3$, computed # # S3[k,j] is $S_{k,j}^{[3]}$ # T3[k,j] is $T_{k,j}^{[3]}$ X3 := factor(subs(so3,X3a)); Y3 := factor(subs(so3,Y3a)); T3 := [seq([seq(0,j=1..3)], k=1..9 )]: for k from 1 to 9 do for j from 1 to 3 do tau3[k,j] := simplify(subs(so3,tau3a[k,j])); od; S3[k,1] := S2[k,1] + simplify( subs(so3,s1star[k]) )*u^3; S3[k,2] := S2[k,2] + simplify( subs(so3,s2star[k]) )*u^3; S3[k,3] := S2[k,3] + simplify( subs(so3,s2a[k]) )*u^3; for j from 1 to 3 do i:=ii[k,j]; tau3[k,j] := simplify(subs(so3,tau3a[k,j])); T3[k,j] := collect( expand(simplify( t[i] + tau3[k,j]*u + S3[k,j]*u^2 )), u ); od; od: T3; X3o := subs(om=o,X3): Y3o := subs(om=o,Y3): T3o := subs(om=o,T3): save X3o,Y3o,T3o, "12tp-approx3.sav"; # ====================================================================================== # # Check that [X3,Y3] is a solution up to $O(u^5)$. # for k from 1 to 9 do for j from 1 to 3 do xx[j] := convert(series( subs( T = T3[k,j], X3 ), u,6),polynom); yy[j] := convert(series( subs( T = T3[k,j], Y3 ), u,6),polynom); od; print( factor(xx[1]-xx[2]), factor(xx[2]-xx[3]), factor(yy[1]-yy[2]), factor(yy[2]-yy[3]) ); od: # ====================================================================================== # # The recursion step [m-1] ---> [m] for m = 4 # # ====================================================================================== su3 := seq( a[k]=subs(so3,a[k])+u*a[k], k=1..9 ), seq( b[k]=subs(so3,b[k])+u*b[k], k=1..9 ); # # su3 replaces the indeterminate $a_{i,3}$ (which is represented by a[i] in the maple code) # by $a_{i,3} + a_{i,4} u$ with the computed value of $a_{i,3}$ and # indeterminate $a_{i,4}$, which is still represented by a[i] in the maple code # X4a := factor(subs(su3,X3a)); Y4a := factor(subs(su3,Y3a)); for k from 1 to 9 do for j from 1 to 3 do tau4a[k,j] := simplify(subs(su3,tau3a[k,j])) od: od: # # [X4a,Y4a] and tau4a[k,j] are $f_{u,a^{[4]}(u)}$ and $\tau(a^{[4]}$ respectively # with indeterminate coefficients $a_{k,4}$ and $b_{k,4}$ # # ====================================================================================== # Compute $s_{k,j}$, $s_{k,j}^*$, and Equation (7). for k from 1 to 9 do i1 := ii[k,1]; i2 := ii[k,2]; i3 := ii[k,3]; se := series( subs( T = t[i1] + tau4a[k,1]*u + (S3[k,1] + s1*u^4)*u^2, X4a ),u,7): xx1 := simplify(expand(op(11,se))); # coef of u^5 se := series( subs( T = t[i2] + tau4a[k,2]*u + (S3[k,2] + s2*u^4)*u^2, X4a ),u,7): xx2 := simplify(expand(op(11,se))); se := series( subs( T = t[i3] + tau4a[k,3]*u + (S3[k,3] + s3*u^4)*u^2, X4a ),u,7): xx3 := simplify(expand(op(11,se))); se := series( subs( T = t[i1] + tau4a[k,1]*u + (S3[k,1] + s1*u^4)*u^2, Y4a ),u,7): yy1 := simplify(expand(op(11,se))); se := series( subs( T = t[i2] + tau4a[k,2]*u + (S3[k,2] + s2*u^4)*u^2, Y4a ),u,7): yy2 := simplify(expand(op(11,se))); se := series( subs( T = t[i3] + tau4a[k,3]*u + (S3[k,3] + s3*u^4)*u^2, Y4a ),u,7): yy3 := simplify(expand(op(11,se))); so13 := solve({xx1=xx3,yy1=yy3},{s1,s3}); so23 := solve({xx2=xx3,yy2=yy3},{s2,s3}); s1star[k] := factor(subs(so13,s1)); ### $s_{k,1}^*$ s2star[k] := factor(subs(so23,s2)); ### $s_{k,2}^*$ s1a[k] := factor(subs(so13,s3)); ### $s_{k,1}$ s2a[k] := factor(subs(so23,s3)); ### $s_{k,2}$ e4[k] := factor(s1a[k]-s2a[k]); ### the equation (7) print([k],expand(simplify(e4[k]))); od: # ====================================================================================== # Solve Equation (7): so4 := solve({a4=a1, a7=0, b1=0, seq(b[k]=0,k = 4..9), seq(e4[k]=0,k=1..9)}); # # subs(so4,a[i]) is $a_{i,4}$ # subs(so4,b[i]) is $b_{i,4}$ # # ====================================================================================== # # [X4,Y4] and tau4[k,j] are $f_{u,a^{[4]}(u)}$ and $\tau(a^{[4]}$ respectively # with all coefficients $a_{i,n}$, $b_{i,n}$, $n=0,...,4$, computed # # S4[k,j] is $S_{k,j}^{[4]}$ # T4[k,j] is $T_{k,j}^{[4]}$ X4 := factor(subs(so4,X4a)); Y4 := factor(subs(so4,Y4a)); T4 := [seq([seq(0,j=1..3)], k=1..9 )]: for k from 1 to 9 do for j from 1 to 3 do tau4[k,j] := simplify(subs(so4,tau4a[k,j])); od; S4[k,1] := S3[k,1] + simplify( subs(so4,s1star[k]) )*u^4; S4[k,2] := S3[k,2] + simplify( subs(so4,s2star[k]) )*u^4; S4[k,3] := S3[k,3] + simplify( subs(so4,s2a[k]) )*u^4; for j from 1 to 3 do i:=ii[k,j]; tau4[k,j] := simplify(subs(so4,tau4a[k,j])); T4[k,j] := collect( expand(simplify( t[i] + tau4[k,j]*u + S4[k,j]*u^2 )), u ); od; od: T4; X4o := subs(om=o,X4): Y4o := subs(om=o,Y4): T4o := subs(om=o,T4): save X4o,Y4o,T4o, "12tp-approx4.sav"; # ====================================================================================== # # Check that [X4,Y4] is a solution up to $O(u^6)$. # for k from 1 to 9 do for j from 1 to 3 do xx[j] := convert(series( subs( T = T4[k,j], X4 ), u,7),polynom); yy[j] := convert(series( subs( T = T4[k,j], Y4 ), u,7),polynom); od; print( [[k]], factor(xx[1]-xx[2]), factor(xx[2]-xx[3]), factor(yy[1]-yy[2]), factor(yy[2]-yy[3]) ); od: # ====================================================================================== # # The recursion step [m-1] ---> [m] for m = 5 # # ====================================================================================== su4 := seq( a[k]=subs(so4,a[k])+u*a[k], k=1..9 ), seq( b[k]=subs(so4,b[k])+u*b[k], k=1..9 ); # # su4 replaces the indeterminate $a_{i,4}$ (which is represented by a[i] in the maple code) # by $a_{i,4} + a_{i,5} u$ with the computed value of $a_{i,4}$ and # indeterminate $a_{i,5}$, which is still represented by a[i] in the maple code # X5a := factor(subs(su4,X4a)); Y5a := factor(subs(su4,Y4a)); for k from 1 to 9 do for j from 1 to 3 do tau5a[k,j] := simplify(subs(su4,tau4a[k,j])) od: od: # # [X5a,Y5a] and tau5a[k,j] are $f_{u,a^{[5]}(u)}$ and $\tau(a^{[5]}$ respectively # with indeterminate coefficients $a_{k,5}$ and $b_{k,5}$ # # ====================================================================================== # Compute $s_{k,j}$, $s_{k,j}^*$, and Equation (7). for k from 1 to 9 do i1 := ii[k,1]; i2 := ii[k,2]; i3 := ii[k,3]; se := series( subs( T = t[i1] + tau5a[k,1]*u + (S4[k,1] + s1*u^5)*u^2, X5a ),u,8): xx1 := simplify(expand(op(13,se))); # coef of u^6 se := series( subs( T = t[i2] + tau5a[k,2]*u + (S4[k,2] + s2*u^5)*u^2, X5a ),u,8): xx2 := simplify(expand(op(13,se))); se := series( subs( T = t[i3] + tau5a[k,3]*u + (S4[k,3] + s3*u^5)*u^2, X5a ),u,8): xx3 := simplify(expand(op(13,se))); se := series( subs( T = t[i1] + tau5a[k,1]*u + (S4[k,1] + s1*u^5)*u^2, Y5a ),u,8): yy1 := simplify(expand(op(13,se))); se := series( subs( T = t[i2] + tau5a[k,2]*u + (S4[k,2] + s2*u^5)*u^2, Y5a ),u,8): yy2 := simplify(expand(op(13,se))); se := series( subs( T = t[i3] + tau5a[k,3]*u + (S4[k,3] + s3*u^5)*u^2, Y5a ),u,8): yy3 := simplify(expand(op(13,se))); so13 := solve({xx1=xx3,yy1=yy3},{s1,s3}); so23 := solve({xx2=xx3,yy2=yy3},{s2,s3}); s1star[k] := factor(subs(so13,s1)); ### $s_{k,1}^*$ s2star[k] := factor(subs(so23,s2)); ### $s_{k,2}^*$ s1a[k] := factor(subs(so13,s3)); ### $s_{k,1}$ s2a[k] := factor(subs(so23,s3)); ### $s_{k,2}$ e5[k] := factor(s1a[k]-s2a[k]); ### the equation (7) print([k],expand(simplify(e5[k]))); od: # ====================================================================================== # Solve Equation (7): so5 := solve({a4=a1, a7=0, b1=0, seq(b[k]=0,k = 4..9), seq(e5[k]=0,k=1..9)}); # # subs(so5,a[i]) is $a_{i,5}$ # subs(so5,b[i]) is $b_{i,5}$ # # ====================================================================================== # # [X5,Y5] and tau5[k,j] are $f_{u,a^{[5]}(u)}$ and $\tau(a^{[5]}$ respectively # with all coefficients $a_{i,n}$, $b_{i,n}$, $n=0,...,5$, computed # # S5[k,j] is $S_{k,j}^{[5]}$ # T5[k,j] is $T_{k,j}^{[5]}$ X5 := factor(subs(so5,X5a)); Y5 := factor(subs(so5,Y5a)); T5 := [seq([seq(0,j=1..3)], k=1..9 )]: for k from 1 to 9 do for j from 1 to 3 do tau5[k,j] := simplify(subs(so5,tau5a[k,j])); od; S5[k,1] := S4[k,1] + simplify( subs(so5,s1star[k]) )*u^5; S5[k,2] := S4[k,2] + simplify( subs(so5,s2star[k]) )*u^5; S5[k,3] := S4[k,3] + simplify( subs(so5,s2a[k]) )*u^5; for j from 1 to 3 do i:=ii[k,j]; tau5[k,j] := simplify(subs(so5,tau5a[k,j])); T5[k,j] := collect( expand(simplify( t[i] + tau5[k,j]*u + S5[k,j]*u^2 )), u ); od; od: X5o := subs(om=o,X5): Y5o := subs(om=o,Y5): T5o := subs(om=o,T5): save X5o,Y5o,T5o, "12tp-approx5.sav"; # ====================================================================================== # # Check that [X5,Y5] is a solution up to $O(u^7)$. # for k from 1 to 9 do for j from 1 to 3 do xx[j] := convert(series( subs( T = T5[k,j], X5 ), u,8),polynom); yy[j] := convert(series( subs( T = T5[k,j], Y5 ), u,8),polynom); od; print( [[k]], factor(xx[1]-xx[2]), factor(xx[2]-xx[3]), factor(yy[1]-yy[2]), factor(yy[2]-yy[3]) ); od: # ====================================================================================== # # The recursion step [m-1] ---> [m] for m = 6 # # ====================================================================================== su5 := seq( a[k]=subs(so5,a[k])+u*a[k], k=1..9 ), seq( b[k]=subs(so5,b[k])+u*b[k], k=1..9 ); # # su5 replaces the indeterminate $a_{i,5}$ (which is represented by a[i] in the maple code) # by $a_{i,5} + a_{i,6} u$ with the computed value of $a_{i,5}$ and # indeterminate $a_{i,6}$, which is still represented by a[i] in the maple code # X6a := factor(subs(su5,X5a)); Y6a := factor(subs(su5,Y5a)); for k from 1 to 9 do for j from 1 to 3 do tau6a[k,j] := simplify(subs(su5,tau5a[k,j])) od: od: # # [X6a,Y6a] and tau6a[k,j] are $f_{u,a^{[6]}(u)}$ and $\tau(a^{[6]}$ respectively # with indeterminate coefficients $a_{k,6}$ and $b_{k,6}$ # # ====================================================================================== # Compute $s_{k,j}$, $s_{k,j}^*$, and Equation (7). for k from 1 to 9 do i1 := ii[k,1]; i2 := ii[k,2]; i3 := ii[k,3]; se := series( subs( T = t[i1] + tau6a[k,1]*u + (S5[k,1] + s1*u^6)*u^2, X6a ),u,9): xx1 := simplify(expand(op(15,se))); # coef of u^7 se := series( subs( T = t[i2] + tau6a[k,2]*u + (S5[k,2] + s2*u^6)*u^2, X6a ),u,9): xx2 := simplify(expand(op(15,se))); se := series( subs( T = t[i3] + tau6a[k,3]*u + (S5[k,3] + s3*u^6)*u^2, X6a ),u,9): xx3 := simplify(expand(op(15,se))); se := series( subs( T = t[i1] + tau6a[k,1]*u + (S5[k,1] + s1*u^6)*u^2, Y6a ),u,9): yy1 := simplify(expand(op(15,se))); se := series( subs( T = t[i2] + tau6a[k,2]*u + (S5[k,2] + s2*u^6)*u^2, Y6a ),u,9): yy2 := simplify(expand(op(15,se))); se := series( subs( T = t[i3] + tau6a[k,3]*u + (S5[k,3] + s3*u^6)*u^2, Y6a ),u,9): yy3 := simplify(expand(op(15,se))); so13 := solve({xx1=xx3,yy1=yy3},{s1,s3}); so23 := solve({xx2=xx3,yy2=yy3},{s2,s3}); s1star[k] := factor(subs(so13,s1)); ### $s_{k,1}^*$ s2star[k] := factor(subs(so23,s2)); ### $s_{k,2}^*$ s1a[k] := factor(subs(so13,s3)); ### $s_{k,1}$ s2a[k] := factor(subs(so23,s3)); ### $s_{k,2}$ e6[k] := factor(s1a[k]-s2a[k]); ### the equation (7) print([k],expand(simplify(e6[k]))); od: # ====================================================================================== # Solve Equation (7): so6 := solve({a4=a1, a7=0, b1=0, seq(b[k]=0,k = 4..9), seq(e6[k]=0,k=1..9)}); # # subs(so6,a[i]) is $a_{i,6}$ # subs(so6,b[i]) is $b_{i,6}$ # # ====================================================================================== # # [X6,Y6] and tau6[k,j] are $f_{u,a^{[6]}(u)}$ and $\tau(a^{[6]}$ respectively # with all coefficients $a_{i,n}$, $b_{i,n}$, $n=0,...,6$, computed # # S6[k,j] is $S_{k,j}^{[6]}$ # T6[k,j] is $T_{k,j}^{[6]}$ X6 := factor(subs(so6,X6a)); Y6 := factor(subs(so6,Y6a)); T6 := [seq([seq(0,j=1..3)], k=1..9 )]: for k from 1 to 9 do for j from 1 to 3 do tau6[k,j] := simplify(subs(so6,tau6a[k,j])); od; S6[k,1] := S5[k,1] + simplify( subs(so6,s1star[k]) )*u^6; S6[k,2] := S5[k,2] + simplify( subs(so6,s2star[k]) )*u^6; S6[k,3] := S5[k,3] + simplify( subs(so6,s2a[k]) )*u^6; for j from 1 to 3 do i:=ii[k,j]; tau6[k,j] := simplify(subs(so6,tau6a[k,j])); T6[k,j] := collect( expand(simplify( t[i] + tau6[k,j]*u + S6[k,j]*u^2 )), u ); od; od: X6o := subs(om=o,X6): Y6o := subs(om=o,Y6): T6o := subs(om=o,T6): save X6o,Y6o,T6o, "12tp-approx6.sav"; # ====================================================================================== # # Check that [X6,Y6] is a solution up to $O(u^8)$. # for k from 1 to 9 do for j from 1 to 3 do xx[j] := convert(series( subs( T = T6[k,j], X6 ), u,9),polynom); yy[j] := convert(series( subs( T = T6[k,j], Y6 ), u,9),polynom); od; print( [[k]], factor(xx[1]-xx[2]), factor(xx[2]-xx[3]), factor(yy[1]-yy[2]), factor(yy[2]-yy[3]) ); od: # ====================================================================================== # # The recursion step [m-1] ---> [m] for m = 7 # # ====================================================================================== su6 := seq( a[k]=subs(so6,a[k])+u*a[k], k=1..9 ), seq( b[k]=subs(so6,b[k])+u*b[k], k=1..9 ); # # su6 replaces the indeterminate $a_{i,6}$ (which is represented by a[i] in the maple code) # by $a_{i,6} + a_{i,7} u$ with the computed value of $a_{i,6}$ and # indeterminate $a_{i,7}$, which is still represented by a[i] in the maple code # X7a := factor(subs(su6,X6a)); Y7a := factor(subs(su6,Y6a)); for k from 1 to 9 do for j from 1 to 3 do tau7a[k,j] := simplify(subs(su6,tau6a[k,j])) od: od: # # [X7a,Y7a] and tau7a[k,j] are $f_{u,a^{[7]}(u)}$ and $\tau(a^{[7]}$ respectively # with indeterminate coefficients $a_{k,7}$ and $b_{k,7}$ # # ====================================================================================== # Compute $s_{k,j}$, $s_{k,j}^*$, and Equation (7). for k from 1 to 9 do i1 := ii[k,1]; i2 := ii[k,2]; i3 := ii[k,3]; se := series( subs( T = t[i1] + tau7a[k,1]*u + (S6[k,1] + s1*u^7)*u^2, X7a ),u,10): xx1 := simplify(expand(op(17,se))); # coef of u^8 se := series( subs( T = t[i2] + tau7a[k,2]*u + (S6[k,2] + s2*u^7)*u^2, X7a ),u,10): xx2 := simplify(expand(op(17,se))); se := series( subs( T = t[i3] + tau7a[k,3]*u + (S6[k,3] + s3*u^7)*u^2, X7a ),u,10): xx3 := simplify(expand(op(17,se))); se := series( subs( T = t[i1] + tau7a[k,1]*u + (S6[k,1] + s1*u^7)*u^2, Y7a ),u,10): yy1 := simplify(expand(op(17,se))); se := series( subs( T = t[i2] + tau7a[k,2]*u + (S6[k,2] + s2*u^7)*u^2, Y7a ),u,10): yy2 := simplify(expand(op(17,se))); se := series( subs( T = t[i3] + tau7a[k,3]*u + (S6[k,3] + s3*u^7)*u^2, Y7a ),u,10): yy3 := simplify(expand(op(17,se))); so13 := solve({xx1=xx3,yy1=yy3},{s1,s3}); so23 := solve({xx2=xx3,yy2=yy3},{s2,s3}); s1star[k] := factor(subs(so13,s1)); ### $s_{k,1}^*$ s2star[k] := factor(subs(so23,s2)); ### $s_{k,2}^*$ s1a[k] := factor(subs(so13,s3)); ### $s_{k,1}$ s2a[k] := factor(subs(so23,s3)); ### $s_{k,2}$ e7[k] := factor(s1a[k]-s2a[k]); ### the equation (7) print([k],expand(simplify(e7[k]))); od: # ====================================================================================== # Solve Equation (7): so7 := solve({a4=a1, a7=0, b1=0, seq(b[k]=0,k = 4..9), seq(e7[k]=0,k=1..9)}); # # subs(so7,a[i]) is $a_{i,7}$ # subs(so7,b[i]) is $b_{i,7}$ # # ====================================================================================== # # [X7,Y7] and tau7[k,j] are $f_{u,a^{[7]}(u)}$ and $\tau(a^{[7]}$ respectively # with all coefficients $a_{i,n}$, $b_{i,n}$, $n=0,...,7$, computed # # S7[k,j] is $S_{k,j}^{[7]}$ # T7[k,j] is $T_{k,j}^{[7]}$ X7 := factor(subs(so7,X7a)); Y7 := factor(subs(so7,Y7a)); T7 := [seq([seq(0,j=1..3)], k=1..9 )]: for k from 1 to 9 do for j from 1 to 3 do tau7[k,j] := simplify(subs(so7,tau7a[k,j])); od; S7[k,1] := S6[k,1] + simplify( subs(so7,s1star[k]) )*u^7; S7[k,2] := S6[k,2] + simplify( subs(so7,s2star[k]) )*u^7; S7[k,3] := S6[k,3] + simplify( subs(so7,s2a[k]) )*u^7; for j from 1 to 3 do i:=ii[k,j]; tau7[k,j] := simplify(subs(so7,tau7a[k,j])); T7[k,j] := collect( expand(simplify( t[i] + tau7[k,j]*u + S7[k,j]*u^2 )), u ); od; od: X7o := subs(om=o,X7): Y7o := subs(om=o,Y7): T7o := subs(om=o,T7): save X7o,Y7o,T7o, "12tp-approx7.sav"; # ====================================================================================== # # Check that [X7,Y7] is a solution up to $O(u^9)$. # for k from 1 to 9 do for j from 1 to 3 do xx[j] := convert(series( subs( T = T7[k,j], X7 ), u,10),polynom); yy[j] := convert(series( subs( T = T7[k,j], Y7 ), u,10),polynom); od; print( [[k]], factor(xx[1]-xx[2]), factor(xx[2]-xx[3]), factor(yy[1]-yy[2]), factor(yy[2]-yy[3]) ); od: