#r04.maple # # Frank Sottile # Jonathon Hauenstein # # This file will study the different four-bar linkages interpolating nine # given points. It begins with numerical solutions to Problem 4 in the # paper of Morgan, Sommese, and Wampler. # # This studies those 2,2 curves with two full circles # interface(quiet=true): with(plots): with(plottools): plotsetup(x11): read("RealSols.maple"): # # Real solutions stored in file as list of Points # mr:=proc(r) return( trunc(1000000*r)/1000000 ) end proc: # # Interpolating poionts # P:=[[0.25, 0.00], [0.52, 0.10], [0.80, 0.70], [1.20, 1.00], [1.40, 1.30], [1.10, 1.48], [0.70, 1.40], [0.20, 1.00], [0.02, 0.40]]: P:=[seq([mr(P[ii][1]),mr(P[ii][2])],ii=1..9)]: Delta:=[seq([P[ii][1]-P[1][1], P[ii][2]-P[1][2]],ii=1..9)]: Deltac:=[seq(P[ii][1]-P[1][1] + I*P[ii][2]-P[1][2],ii=1..9)]: # 3, #j:=1 # # # 14,15,16,28,30,54,60 Tw:=[22,23,28,29,35,42,43,51,58,60,68,85,88,96,106,114,118,120,126,129,130,131,132,134,138,142,148,153,157,160,162,169,175,180,181,194,195,204,207,209,210,213,219,222,241,266,268,276,278,281,284,287,288,291,301,304,310,318,324,335,338,343,354,364]: i:=Tw[j]; S:=[seq(mr(Points[i][ii]),ii=1..12)]: #Digits:=10: #S:=[seq(Points[i][ii],ii=1..12)]: # The four vectors x:=[S[1],S[2]]: a:=[S[3],S[4]]: n:=[S[5],S[6]]: y:=[S[7],S[8]]: b:=[S[9],S[10]]: m:=[S[11],S[12]]: # Make them into complex numbers xc:=x[1]+I*x[2]: ac:=a[1]+I*a[2]: nc:=n[1]+I*n[2]: yc:=y[1]+I*y[2]: bc:=b[1]+I*b[2]: mc:=m[1]+I*m[2]: # # Check that equation (9) holds. # evalf(simplify(ac*conjugate(xc)-nc)): evalf(simplify(bc*conjugate(yc)-mc)): # # Check that the remaining equations hold # for j from 2 to 9 do solve({(conjugate(ac)-conjugate(Deltac[j]))*xc*gam + (ac-Deltac[j])*conjugate(xc)*cgam + Deltac[j]*(conjugate(ac)-conjugate(xc)) + conjugate(Deltac[j])*(ac-xc) -Deltac[j]*conjugate(Deltac[j]), (conjugate(bc)-conjugate(Deltac[j]))*yc*gam + (bc-Deltac[j])*conjugate(yc)*cgam + Deltac[j]*(conjugate(bc)-conjugate(yc)) + conjugate(Deltac[j])*(bc-yc) -Deltac[j]*conjugate(Deltac[j])}, {gam,cgam}): end do: ########################################################################### # The Anchor points A:=[P[1][1]+a[1], P[1][2]+a[2]]: B:=[P[1][1]+b[1], P[1][2]+b[2]]: # # Two legs of four-bar mechanism, as vectors # u:=[x[1]-a[1], x[2]-a[2]]: v:=[y[1]-b[1], y[2]-b[2]]: # # The w-bar (as a vector), and its length # w:=[y[1]-x[1], y[2]-x[2]]: Pars:=solve({al*w[1]+be*w[2]=x[1], al*w[2]-be*w[1]=x[2]},{al,be}): ##################################################################### # # p and q are the positions of the endpoints of the u-bar and the v-bar # p := (s,t) -> [A[1]+(2*s*t*u[1]+(s^2-t^2)*u[2])/(s^2+t^2), A[2]+(2*s*t*u[2]-(s^2-t^2)*u[1])/(s^2+t^2)]: q := (s,t) -> [B[1]+(2*s*t*v[1]+(s^2-t^2)*v[2])/(s^2+t^2), B[2]+(2*s*t*v[2]-(s^2-t^2)*v[1])/(s^2+t^2)]: # # The one equation is that dist(p,q) = |w|. # W:=[q(t,1)[1]-p(s,1)[1],q(t,1)[2]-p(s,1)[2]]: X:=simplify(subs(op(Pars),[p(s,1)[1]-al*W[1]-be*W[2], p(s,1)[2]-al*W[2]+be*W[1]])): G:=(p(s,1)[1]-q(t,1)[1])^2+(p(s,1)[2]-q(t,1)[2])^2-(w[1]^2+w[2]^2): G:=numer(simplify(G)): #implicitplot(G,s=-8..8, t=-8..8,grid=[100,100],color=red); Gi:=(p(s,1)[1]-q(1,t)[1])^2+(p(s,1)[2]-q(1,t)[2])^2-(w[1]^2+w[2]^2): Gi:=numer(simplify(Gi)): #implicitplot(Gi,s=-8..8, t=-8..8,grid=[100,100],color=blue); At:=coeff(G,t^2): Bt:=coeff(G,t): Ct:=subs(t=0,G): t1:= unapply((-Bt + sqrt(Bt^2-4*At*Ct))/2/At,s): t2:= unapply((-Bt - sqrt(Bt^2-4*At*Ct))/2/At,s): ##################################################################### x:=unapply(X,s,t): s0:=0:s1:=1/2: fb:=[plottools[curve]([A,p(s0,1),q(t1(s0),1),B],color=green), plottools[curve]([p(s0,1),x(s0,t1(s0)),q(t1(s0),1)],color=red), plottools[curve]([A,p(s1,1),q(t2(s1),1),B],color=green), plottools[curve]([p(s1,1),x(s1,t2(s1)),q(t2(s1),1)],color=red)]: Gr:=[pointplot([A,B],symbol=cross,color=black), pointplot([A,B],symbol=circle,color=red)]: Po:=pointplot(P,symbol=circle,color=blue): plotP1:=[plot([p(s,1)[], s=-1..1], color=green), plot([p(1,s)[], s=-1..1], color=green)]: plotQ1:=[plot([q(t1(s),1)[], s=-1..1], color=blue), plot([q(t1(1/s),1)[], s=-1..1], color=blue)]: plotR1:=[plot([x(s,t1(s))[], s=-1..1], color=red), plot([x(1/s,t1(1/s))[], s=-1..1], color=red)]: #display(plotP1,plotQ1,plotR1,Po,Gr,scaling=constrained); plotP2:=[plot([p(s,1)[], s=-1..1], color=green), plot([p(1,s)[], s=-1..1], color=green)]: plotQ2:=[plot([q(t2(s),1)[], s=-1..1], color=blue), plot([q(t2(1/s),1)[], s=-1..1], color=blue)]: plotR2:=[plot([x(s,t2(s))[], s=-1..1], color=brown), plot([x(1/s,t2(1/s))[], s=-1..1], color=brown)]: #display(plotP1,plotQ1,plotR1,plotP2,plotQ2,plotR2,Po,Gr,scaling=constrained); display(plotR1,plotR2,fb,Po,Gr,scaling=constrained);