{VERSION 5 0 "IBM INTEL NT" "5.0" } {USTYLETAB {CSTYLE "Maple Input" -1 0 "Courier" 0 1 255 0 0 1 0 1 0 0 1 0 0 0 0 1 }{PSTYLE "Normal" -1 0 1 {CSTYLE "" -1 -1 "Times" 1 12 0 0 0 1 2 2 2 2 2 2 1 1 1 1 }1 1 0 0 0 0 1 0 1 0 2 2 0 1 }{PSTYLE "Maple Plot" -1 13 1 {CSTYLE "" -1 -1 "Times" 1 12 0 0 0 1 2 2 2 2 2 2 1 1 1 1 }3 1 0 0 0 0 1 0 1 0 2 2 0 1 }} {SECT 0 {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 26 "with(linalg): with(p lots):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 121 "# Computes the vector of second derivati ves in (2); x are the k+2 variables,\n# xp are the coresponding first \+ derivatives:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 1047 "ddot:=pro c(x,xp,n)local k,M,r,i,j;\nk:=nops(convert(x,list))-2; M:=matrix(k+2,k +2); r:=vector(k+2);\nfor i to k do for j to k do M[i,j]:=min(i,j)*cos (x[i]-x[j]) end do end do; for i to k do M[k+1,i]:=-i*sin(x[i]);M[i,k+ 1]:=M[k+1,i];\nM[i,k+2]:=i*cos(x[i]);M[k+2,i]:=M[i,k+2];\nr[i]:=-i*sin (x[i])-sum('xp[j]^2*min(i,j)*sin(x[i]-x[j])','j'=1..k) end do;\nr[k+1] :=k+1+sum('j*xp[j]^2*cos(x[j])','j'=1..k);\nr[k+2]:=sum('j*xp[j]^2*sin (x[j])','j'=1..k);\nif x[k+1]=0 and x[k+2]=0 then M[k+1,k+1]:=(n-k)*xp [k+1]^2/(xp[k+1]^2+xp[k+2]^2)+k+1;\nM[k+2,k+2]:=(n-k)*xp[k+2]^2/(xp[k+ 1]^2+xp[k+2]^2)+k+1;\nM[k+1,k+2]:=(n-k)*xp[k+1]*xp[k+2]/(xp[k+1]^2+xp[ k+2]^2);M[k+2,k+1]:=M[k+1,k+2];\nelse M[k+1,k+1]:=(n-k)*x[k+1]^2/(x[k+ 1]^2+x[k+2]^2)+k+1;\nM[k+2,k+2]:=(n-k)*x[k+2]^2/(x[k+1]^2+x[k+2]^2)+k+ 1;\nM[k+1,k+2]:=(n-k)*x[k+1]*x[k+2]/(x[k+1]^2+x[k+2]^2);M[k+2,k+1]:=M[ k+1,k+2];\nr[k+1]:=r[k+1]-(n-k)*x[k+1]*(xp[k+2]*x[k+1]-xp[k+1]*x[k+2]) ^2/(x[k+1]^2+x[k+2]^2)^2;\nr[k+2]:=r[k+2]+(k-n)*x[k+2]*(xp[k+2]*x[k+1] -xp[k+1]*x[k+2])^2/(x[k+1]^2+x[k+2]^2)^2 end if;\nlinsolve(M,r) end pr oc:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 " > " 0 "" {MPLTEXT 1 0 35 "# Advances x and xp by time-step h:" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 288 "kutta:=proc(x,xp,n,h) local k1,k2,k3,k4; k1:=h*ddot(x,xp,n);\nk2:=h*ddot(evalm(x+h/2*xp+h/8*k1),e valm(xp+k1/2),n);\nk3:=h*ddot(evalm(x+h/2*xp+h/8*k1),evalm(xp+k2/2),n) ;\nk4:=h*ddot(evalm(x+h*xp+h/2*k3),evalm(xp+k3),n);\n[evalm(x+h*(xp+k1 /6+k2/6+k3/6)),evalm(xp+k1/6+k2/3+k3/3+k4/6)] end proc:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 147 "# Computes the vector of tau derivatives in (3); x i s now extended by d,\n# xp is extended by Omega, and by the 'infinites imal' values of Xn and Yn:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 738 "impu:=proc(x,xp,n)local k,M,r,i,j;\nk:=nops(convert(x,list))-3; M :=matrix(k+3,k+3); r:=vector(k+3);\nfor i to k do for j to k do M[i,j] :=min(i,j)*cos(x[i]-x[j]) end do end do;\nfor i to k do M[k+1,i]:=-i*s in(x[i]); M[i,k+1]:=M[k+1,i];\nM[i,k+2]:=i*cos(x[i]);M[k+2,i]:=M[i,k+2 ]; M[i,k+3]:=-i*sin(x[i]-x[k]);\nM[k+3,i]:=M[i,k+3]; r[i]:=0 end do;\n r[k+1]:=0;r[k+2]:=0;r[k+3]:=-1;\nM[k+1,k+1]:=(n-k)*xp[k+4]^2/(xp[k+4]^ 2+xp[k+5]^2)+k+1;\nM[k+2,k+2]:=(n-k)*xp[k+5]^2/(xp[k+4]^2+xp[k+5]^2)+k +1;\nM[k+1,k+2]:=(n-k)*xp[k+4]*xp[k+5]/(xp[k+4]^2+xp[k+5]^2);M[k+2,k+1 ]:=M[k+1,k+2];\nM[k+1,k+3]:=k*cos(x[k]); M[k+3,k+1]:=M[k+1,k+3];\nM[k+ 2,k+3]:=k*sin(x[k]); M[k+3,k+2]:=M[k+2,k+3]; M[k+3,k+3]:=k; M:=linsolv e(M,r);\n[seq(M[i],i=1..k+3),xp[k+1],xp[k+2]] end proc:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 70 "# Advances the momenta (plus;infinitesimal' Xn and Yn ) by time step h:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 194 "runge :=proc(x,xp,n,h) local k1,k2,k3,k4; k1:=h*impu(x,xp,n);\nk2:=h*impu(x, evalm(xp+k1/2),n);\nk3:=h*impu(x,evalm(xp+k2/2),n);\nk4:=h*impu(x,eval m(xp+k3),n);\nevalm(xp+k1/6+k2/3+k3/3+k4/6) end proc:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 75 "# Advances x and xp, in about 10 time steps, until sq rt(X^2+Y^2) reaches 1:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 383 " step:=proc(xn,n)local h,k,x; k:=nops(convert(xn[1],list))-2; x:=xn; \nwhile abs(1-sqrt(x[1][k+1]^2+x[1][k+2]^2))>0.0001 do h:=min(rhs(solv e(\{sqrt(x[2][k+1]^2+x[2][k+2]^2)*h+(k+1)/(n+1)*h^2-.1,h>0\},h)[1]),\n (1-sqrt(x[1][k+1]^2+x[1][k+2]^2))*sqrt(x[1][k+1]^2+x[1][k+2]^2+.000000 1)/\n(x[1][k+1]*x[2][k+1]+x[1][k+2]*x[2][k+2]+.0000001));\nx:=kutta(x[ 1],x[2],n,h); h:='h' end do; x end proc:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 68 "# Advan ces the momenta (in about 10 steps) until Omega becomes zero:" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 293 "step2:=proc(xn,n)local h,k, x,a,i; k:=nops(convert(xn[1],list))-3; x:=xn; i:=10;\nwhile abs(x[2] [k+3])>0.0001 do a:=impu(x[1],x[2],n)[k+3];\nh:=min(-x[2][k+3]/a/(i+.0 1),-x[2][k+3]/a); i:=i-1;\nx[2]:=runge(x[1],x[2],n,h); h:='h' end do; \n[[seq(x[1][i],i=1..k+2)],[seq(x[2][i],i=1..k+2)]] end proc:" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 155 "# When sqrt(X^2+Y^2) reaches 1, a new angle is int roduced, Xn and Yn are set to 0\n# and d is set to 1, with the corresp onding momenta properly initialized:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 302 "fix:=proc(x,xp) local k,dd,phi,arcd,xn,xpn,i; k:=nop s(convert(x,list))-2;\ndd:=(x[k+1]*xp[k+1]+x[k+2]*xp[k+2]);\narcd:=(x[ k+1]*xp[k+2]-x[k+2]*xp[k+1]);\nphi:=arctan(x[k+2]/x[k+1]); xn:=[seq(x[ i],i=1..k),phi,0,0,1];\nxpn:=[seq(xp[i],i=1..k),arcd-dd*x[k+1],0,dd,dd *(1-x[k+2]),0,0.0000001];\n[xn,xpn] end proc:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 65 "# \+ Plots the points (assuming that we have the Xn=Yn=0 situation)." }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 137 "points:=proc(x) local k,i,j ; k:=nops(convert(x,list))-2;\n[seq([sum('sin(x[i])','i'=j..k),sum('-c os(x[i])','i'=j..k)],j=1..k+1)] end proc:" }}}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 100 "# Mo ves the chain of n links, from beginning till end (100 links will take several hours,\n# try 25):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 33 "n:=100: xxp:=[[0.00001,0],[0,0]]:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 82 "for j to n do xxp:=step(xxp,n); xxp:=fix(xxp[1],xxp[2 ]); xxp:=step2(xxp,n) end do:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 29 "# Displayes the fi nal answer:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 45 "listplot(poi nts(xxp[1]),scaling=CONSTRAINED);" }}{PARA 13 "" 1 "" {GLPLOT2D 857 694 694 {PLOTDATA 2 "6$-%'CURVESG6#7aq7$$!+oSo$o#!#5$!+gq&4\\)!\")7$$! +o!))))G\"F*$!+kX$>R)F-7$$!*oA) F-7$$\"*\\&*)yPF*$!+SSe#4)F-7$$\"+\\<$)oCF*$!+!\\%z%*zF-7$$\"+\\T\"*yC F*$!+SXz%*yF-7$$\"+\\$*)=m\"F*$!+i)G^z(F-7$$\"+\\b8#z\"F*$!+Yt8&p(F-7$ $\"+\\2JeFF*$!+'=0cf(F-7$$\"+\\R3pPF*$!+It6'\\(F-7$$\"+\\$p')G%F*$!+1C D'R(F-7$$\"+\\Db1WF*$!+a$fiH(F-7$$\"+\\J0uWF*$!+K;E'>(F-7$$\"+\\d?RZF* $!+!z'H'4(F-7$$\"+\\&R$o_F*$!+!)oV'*pF-7$$\"+\\V&G&fF*$!+O9n'*oF-7$$\" +\\p$zi'F*$!+k&**oz'F-7$$\"+\\\\!Q=(F*$!+yT0(p'F-7$$\"+\\>*))f(F*$!+k. 9(f'F-7$$\"+\\BXczF*$!+5V?(\\'F-7$$\"+\\Ta4%)F*$!+3qI(R'F-7$$\"+\\RU$* *)F*$!+7wZ(H'F-7$$\"+zPoe'*F*$!+U\"*p(>'F-7$$\"+Qn/R5!\"*$!+as'z4'F-7$ $\"+)>ho6\"F^s$!+o/F)*fF-7$$\"+yB'*)>\"F^s$!+o!3')*eF-7$$\"+)*3Q(G\"F^ s$!+G(***)z&F-7$$\"+eP#QQ\"F^s$!+ueY*p&F-7$$\"+=A%G[\"F^s$!+9t&**f&F-7 $$\"+)>uzd\"F^s$!+[3T+bF-7$$\"+yl>r;F^s$!+;j%3S&F-7$$\"+y`agF^s$!+!Q'y,^F-7$$\"+\"ob]+#F ^s$!+!>kA+&F-7$$\"+h'[a6#F^s$!+!RvG!\\F-7$$\"+@_*zA#F^s$!+W2^.[F-7$$\" +,p!fO#F^s$!+#HmWq%F-7$$\"+@:c>DF^s$!+IQl0YF-7$$\"+J[#*eEF^s$!+)pHm]%F -7$$\"+h>&**y#F^s$!+?=\\2WF-7$$\"+.+G?HF^s$!+EZM3VF-7$$\"+`=KQIF^s$!+k Q/4UF-7$$\"+$*yr^JF^s$!+y))o4TF-7$$\"+tH-\"G$F^s$!+#RG0,%F-7$$\"+$efgU $F^s$!+kde6RF-7$$\"+RiInNF^s$!+?$)e7QF-7$$\"+R+')*p$F^s$!+Q2Z8PF-7$$\" +*)RlSQF^s$!+YoY9OF-7$$\"+8B\\-SF^s$!+=^y:NF-7$$\"+`Qd%>%F^s$!+?skLF-7$$\"+hW?'e%F^s$!+u'G:A$F-7$$\"+VKPdZF^s$!+ -X+BJF-7$$\"+nC*y#\\F^s$!+m!pW-$F-7$$\"+;]Z)3&F^s$!+EowDHF-7$$\"+_^X@_ F^s$!+c\\lEGF-7$$\"+!QO$[`F^s$!+iJYFFF-7$$\"+)Rl[b&F^s$!+;\">'HEF-7$$ \"+#Rem!fF^s$!+K8,ODF-7$$\"+2LYVjF^s$!+1d0YCF-7$$\"+h;2(o'F^s$!+-W9_BF -7$$\"+\"R(ornF^s$!+QI]_AF-7$$\"+DT$Hl'F^s$!+e1@`@F-7$$\"+:uQ-mF^s$!+) [QL0#F-7$$\"+&zm@(oF^s$!+si/d>F-7$$\"+$pHiX(F^s$!+Ka(e(=F-7$$\"+p7J9\" )F^s$!+o2e+=F-7$$\"+ZzHV&)F^s$!+'o\\-r\"F-7$$\"+>n1i&)F^s$!+-tE5;F-7$$ \"+J'R\\Q)F^s$!+/&[=^\"F-7$$\"+j@R\"R)F^s$!+'e]=T\"F-7$$\"+N+/r()F^s$! +zuL>8F-7$$\"+\">e2S*F^s$!+4^lT7F-7$$\"+u)*4-5F-$!+UR@j6F-7$$\"+l(*zP5 F-$!+iM!)p5F-7$$\"+.C5S5F-$!+V(*H)p*F^s7$$\"+%[:?.\"F-$!+>]d,()F^s7$$ \"+z-CU5F-$!+vg\"oq(F^s7$$\"+dwH$3\"F-$!+(p))\\z'F^s7$$\"+p=.\\6F-$!+8 kRTgF^s7$$\"+R#)GE7F-$!+T\"fkS&F^s7$$\"+R-r-8F-$!+(\\.:w%F^s7$$\"+XSKj 8F-$!+VS9mRF^s7$$\"+Ik(fQ\"F-$!+D%Q@*HF^s7$$\"+:\\Rb8F-$!+>v/S?F^s7$$ \"+ap^$G\"F-$!+B&4[M\"F^s7$$\"+#ePI>\"F-$!+@Tf*=*F*7$$\"+#32l4\"F-$!+& QW$ylF*7$$\"+$QQ1)**F^s$!+\"**42#[F*7$$\"+Qb!*))*)F^s$!+!H7v`$F*7$$\"+ Nym$*zF^s$!+5YtiDF*7$$\"+s,['*pF^s$!+_QG8=F*7$$\"+uI8)*fF^s$!+)4g&Q7F* 7$$\"+G*z!**\\F^s$!+\\t&\\.)!#67$$\"+6ff**RF^s$!+I8&H#[Fdil7$$\"+jN&)* *HF^s$!+=Q1`DFdil7$$\"+qL'***>F^s$!+$))o62\"Fdil7$$\"+7y'*****F*$!+t5D PD!#77$$\"\"!F\\[mF[[m-%(SCALINGG6#%,CONSTRAINEDG" 1 2 0 1 10 0 2 9 1 4 1 1.000000 45.000000 45.000000 0 0 "Curve 1" }}}}}{MARK "26 0 0" 89 }{VIEWOPTS 1 1 0 1 1 1803 1 1 1 1 }{PAGENUMBERS 0 1 2 33 1 1 }