> restart; with(LinearAlgebra); assume(omega, real, omega > 0); > G := 9; > z := (xi^2+xi/(1+xi^2))/(1+xi^2); `output redirected...`> print(); # input placeholder > C := `<,>`(1-z, seq(sin((n-1)*Pi*z), n = 2 .. G)); `output redirected...`> print(); # input placeholder > g := Transpose(C); `output redirected...`> print(); # input placeholder > A := Multiply(C, g); `output redirected...`> print(); # input placeholder > B := xi^2*A; `output redirected...`> print(); # input placeholder > for nn to G do for hh to G do A[nn, hh] := evalf(int(A[nn, hh], xi = 0 .. infinity)); B[nn, hh] := evalf(int(B[nn, hh], xi = 0 .. infinity)) end do end do; > for nn to G do C(nn, 1) := evalf(int(C(nn, 1), xi = 0 .. infinity)) end do; > f := Transpose(C); > H := MatrixInverse(A); `output redirected...`> print(); # input placeholder > M := Multiply(H, B); `output redirected...`> print(); # input placeholder > N := Multiply(H, C); > gh := `<|>`(N, Vector(1 .. G, 0), -M); `output redirected...`> print(); # input placeholder > fgh2 := `<|>`(h1, K1-1, (7*K2*(1/11))*f); > cv := `<|>`(1, Vector[row](1 .. G+1, 0)); `output redirected...`> print(); # input placeholder > capB := `<,>`(fgh2, cv, gh); `output redirected...`> print(); # input placeholder > Sb := CharacteristicMatrix(capB, s); `output redirected...`> print(); # input placeholder > carpol2 := LinearAlgebra[Determinant](Sb); `output redirected...`> print(); # input placeholder > eq2 := subs(s = I*omega, carpol2); `output redirected...`> print(); # input placeholder > eq21 := evalc(Re(eq2)); `output redirected...`> print(); # input placeholder > eq22 := evalc(Im(eq2)); `output redirected...`> print(); # input placeholder > ee2 := solve({eq21, eq22}, {K1, K2}); `output redirected...`> print(); # input placeholder > hio := subs(ee2[1], ee2[2], fgh2); `output redirected...`> print(); # input placeholder > capA := `<,>`(hio, cv, gh); > kv := `<,>`(1, seq(v[nn], nn = 2 .. G+2)); `output redirected...`> print(); # input placeholder > Iden := IdentityMatrix(G+2); `output redirected...`> print(); # input placeholder > junk := Multiply(capA-lambda*Iden, kv); > %; > righteigenvector := solve(subs(lambda = I*omega, {seq(junk(nn, 1), nn = 2 .. G+2)}), {seq(v[nn], nn = 2 .. G+2)}); `output redirected...`> print(); # input placeholder > mo := `<,>`(1, seq(m[nn], nn = 2 .. G+2)); `output redirected...`> print(); # input placeholder > junke := Multiply(capA-lambda*Iden, mo); > righteigenvectorminus := solve(subs(lambda = -I*omega, {seq(junke(nn, 1), nn = 2 .. G+2)}), {seq(m[nn], nn = 2 .. G+2)}); `output redirected...`> print(); # input placeholder > kw := `<|>`(1, Vector[row]([seq(w[nn], nn = 2 .. G+2)])); `output redirected...`> print(); # input placeholder > junk1 := Multiply(kw, capA-lambda*Iden); > lefteigenvector := solve(subs(lambda = I*omega, {seq(junk1(1, nn), nn = 1 .. G+1)}), {seq(w[nn], nn = 2 .. G+2)}); `output redirected...`> print(); # input placeholder > kwo := `<|>`(1, Vector[row]([seq(wo[nn], nn = 2 .. G+2)])); `output redirected...`> print(); # input placeholder > junke1 := Multiply(kwo, capA-lambda*Iden); > lefteigenvectorminus := solve(subs(lambda = -I*omega, {seq(junke1(1, nn), nn = 1 .. G+1)}), {seq(wo[nn], nn = 2 .. G+2)}); `output redirected...`> print(); # input placeholder > alpha := `<,>`(seq(diff(y[n](t), t), n = 1 .. G+2)); `output redirected...`> print(); # input placeholder > beta := `<,>`(seq(y[n](t), n = 1 .. G+2)); `output redirected...`> print(); # input placeholder > ee := alpha-Multiply(capB, beta); > de[1] := ee(1, 1)-epsilon*h2*y[1](t)^2-epsilon^2*(h3*y[1](t)^3+y[2](t)*Delta); `output redirected...`> print(); # input placeholder > de[2] := ee(2, 1); `output redirected...`> print(); # input placeholder > for j to G do de[j+2] := ee(j+2, 1) end do; `output redirected...`> print(); # input placeholder > vm11 := Vector(G+2, [1, seq(v[nn], nn = 2 .. G+2)]); `output redirected...`> print(); # input placeholder > vm12 := Vector(G+2, [1, seq(m[nn], nn = 2 .. G+2)]); `output redirected...`> print(); # input placeholder > wm11 := Vector(G+2, [1, seq(w[nn], nn = 2 .. G+2)]); # # > wm12 := Vector(G+2, [1, seq(wo[nn], nn = 2 .. G+2)]); > for j to G+2 do for i to G+2 do expand(subs(y[i](t) = y[i](seq(T[k](t), k = 0 .. 2)), de[j])); expand(subs({seq(diff(T[k](t), t) = epsilon^k, k = 0 .. 2)}, %)); convert(taylor(expand(subs({seq(T[k](t) = T[k], k = 0 .. 2)}, %)), epsilon = 0, 3), polynom); de[j] := subs((D[1](y[j]))(seq(T[mmm], mmm = 0 .. 2)) = diff(y[j](seq(T[s], s = 0 .. 2)), T[0]), (D[2](y[j]))(seq(T[mmm], mmm = 0 .. 2)) = diff(y[j](seq(T[s], s = 0 .. 2)), T[1]), (D[3](y[j]))(seq(T[mmm], mmm = 0 .. 2)) = diff(y[j](seq(T[s], s = 0 .. 2)), T[2]), %) end do end do; > for ff to G+2 do for j to G+2 do de[ff] := convert(taylor(expand(subs(y[j](seq(T[ss], ss = 0 .. 2)) = add(y[j, tt](T[0], T[1], T[2])*epsilon^tt, tt = 0 .. 2), de[ff])), epsilon = 0, 3), polynom) end do end do; `output redirected...`> print(); # input placeholder > for ff to G+2 do for j to G+2 do de[ff] := convert(taylor(combine(expand(subs(y[j, 0](seq(T[ss], ss = 0 .. 2)) = B1(T[1], T[2])*vm11(j, 1)*exp(I*omega*T[0])+B2(T[1], T[2])*vm12(j, 1)*exp(-I*omega*T[0]), de[ff]))), epsilon, 3), polynom) end do end do; > for ff to G+2 do h1[ff] := coeff(coeff(de[ff], epsilon, 1), exp((1.*I)*omega*T[0])) end do; > for ff to G+2 do h[ff] := coeff(coeff(de[ff], epsilon, 1), exp(-(1.*I)*omega*T[0])) end do; > junk1 := sum('wm11[kk]*h1[kk]', 'kk' = 1 .. G+2); > junk2 := sum('wm12[kk]*h[kk]', 'kk' = 1 .. G+2); > eqs := solve({junk1, junk2}, {diff(B1(T[1], T[2]), T[1]), diff(B2(T[1], T[2]), T[1])}); ??> print(); # input placeholder > for ff to G+2 do t1[ff] := coeff(de[ff], epsilon, 1) end do; > for j to G+2 do z1sol[j] := {y[j, 1](T[0], T[1], T[2]) = B1(T[1], T[2])^2*s[3*j-2]*exp((2*I)*omega*T[0])+B2(T[1], T[2])^2*s[3*j-1]*exp(-(2*I)*omega*T[0])+s[3*j]*B1(T[1], T[2])*B2(T[1], T[2])} end do; > for ff to G+2 do t1[ff] := subs(eqs, convert(taylor(combine(expand(subs(seq(z1sol[j], j = 1 .. G+2), t1[ff]))), epsilon, 3), polynom)) end do; > for ff to G+2 do t1[ff] := subs(exp((2*I)*omega*T[0]) = ev2, exp((2.*I)*omega*T[0]) = ev2, exp(-(2*I)*omega*T[0]) = en2, exp(-(2.*I)*omega*T[0]) = en2, exp((1.*I)*omega*T[0]) = ev1, exp(I*omega*T[0]) = ev1, exp(-(1.*I)*omega*T[0]) = en1, exp(-I*omega*T[0]) = en1, t1[ff]) end do; > for rr to G+2 do eqe(3*rr-2) := coeff(t1[rr], ev2); eqe(3*rr-1) := coeff(t1[rr], en2); eqe(3*rr) := subs(ev2 = 0, en2 = 0, ev1 = 0, en1 = 0, t1[rr]) end do; > sol := solve({seq(eqe(mm), mm = 1 .. 3*(G+2))}, {seq(s[nn], nn = 1 .. 3*(G+2))}); ??> print(); # input placeholder > for ff to G+2 do t2[ff] := simplify(subs(seq(z1sol[j], j = 1 .. G+2), coeff(de[ff], epsilon, 2))) end do; for ff to G+2 do t2[ff] := subs(exp((2*I)*omega*T[0]) = ev2, exp((2.*I)*omega*T[0]) = ev2, exp(-(2*I)*omega*T[0]) = en2, exp(-(2.*I)*omega*T[0]) = en2, exp((1.*I)*omega*T[0]) = ev1, exp(I*omega*T[0]) = ev1, exp(-(1.*I)*omega*T[0]) = en1, exp(-I*omega*T[0]) = en1, t2[ff]) end do; > for ff to G+2 do h1[ff] := coeff(t2[ff], ev1) end do; > for ff to G+2 do h[ff] := coeff(t2[ff], en1) end do; > W1 := subs(sol, ee2, righteigenvector, righteigenvectorminus, lefteigenvector, lefteigenvectorminus, sum('wm11[kk]*h1[kk]', 'kk' = 1 .. G+2)); > W2 := subs(sol, ee2, righteigenvector, righteigenvectorminus, lefteigenvector, lefteigenvectorminus, sum('wm12[kk]*h[kk]', 'kk' = 1 .. G+2)); > junk11 := subs(phi(T[1], T[2]) = phi, R(T[1], T[2]) = R, subs(diff(R(T[1], T[2]), T[2]) = Rdot, diff(phi(T[1], T[2]), T[2]) = phidot, eval(subs(B1(T[1], T[2]) = (1/2)*R(T[1], T[2])*exp((1.*I)*phi(T[1], T[2])), B2(T[1], T[2]) = (1/2)*R(T[1], T[2])*exp(-(1.*I)*phi(T[1], T[2])), W1)))); > junk22 := subs(phi(T[1], T[2]) = phi, R(T[1], T[2]) = R, subs(diff(R(T[1], T[2]), T[2]) = Rdot, diff(phi(T[1], T[2]), T[2]) = phidot, eval(subs(B1(T[1], T[2]) = (1/2)*R(T[1], T[2])*exp((1.*I)*phi(T[1], T[2])), B2(T[1], T[2]) = (1/2)*R(T[1], T[2])*exp(-(1.*I)*phi(T[1], T[2])), W2)))); > eqsubs1 := {CR11 = coeff(subs(Rdot = 0, phidot = 0, junk11), R), CR13 = coeff(subs(Rdot = 0, phidot = 0, junk11), R^3), Cpd1 = simplify(coeff(junk11, phidot)), Crd1 = simplify(coeff(junk11, Rdot))}; > eqsubs2 := {CR21 = coeff(subs(Rdot = 0, phidot = 0, junk22), R), CR23 = coeff(subs(Rdot = 0, phidot = 0, junk22), R^3), Cpd2 = simplify(coeff(junk22, phidot)), Crd2 = simplify(coeff(junk22, Rdot))}; > eqsf := {Crd1*Rdot+Cpd1*phidot+CR11*R+CR13*R^3, Crd2*Rdot+Cpd2*phidot+CR21*R+CR23*R^3}; > junksol := solve(eqsf, {Rdot, phidot}); > Rsol := subs(junksol, Rdot); > phisol := subs(junksol, phidot); > Rsol := subs(eqsubs1, eqsubs2, Rsol); > %; > phisol := subs(eqsubs1, eqsubs2, phisol); > tyu := expand(simplify(subs(omega = 1, h1 = 0.646e-1, h2 = 0.293e-1, h3 = -0.98e-1, Rsol))); ??> print(); # input placeholder > limcye := evalc(Re(tyu)); ??> print(); # input placeholder > limcy := solve(limcye, R); ??> print(); # input placeholder >