1
$\begingroup$

I got the following function:

f[q1p1_, q2_, p2_, q3_, p3_] = 
 0.9816 q1p1 + 1.1822 I q2 p2 + 1.2514 I q3 p3 - 0.0889594 I p2^3 - 
  0.0970894 I p2 p3^2 - 0.219332 p1 p2 q1 + 0.256743 I p2^2 q2 - 
  0.0970894 I p3^2 q2 + 0.0256743 I p2 q2^2 - 0.0889594 I q2^3 - 
  0.194179 p2 p3 q3 - 0.194179 p3 q2 q3 + 0.0970894 I p2 q3^2 + 
  0.0970894 I q2 q3^2 

and I want to plot the Poincaré section of $f$ where $q1p1=q3=0$ such that $f=0.4$

$\endgroup$
4
  • $\begingroup$ Whenever I've heard the phrase "Poincaré section", it has referred to a slice of the phase space of a dynamical system. Does your problem involve differential equations? $\endgroup$
    – Chris K
    Commented Oct 14, 2016 at 12:12
  • $\begingroup$ Yes it refers to a slice of a certain dynamical system. This is a canonical transformed hamiltonian system $\endgroup$
    – Tigran
    Commented Oct 14, 2016 at 12:17
  • 1
    $\begingroup$ q1 and p1 show up as separate parameters in this equation (see the 6th term), but they are not independent arguments to the function. Is there an error? Also, three variables remain independent; which slice do you want? $\endgroup$ Commented Oct 14, 2016 at 12:54
  • $\begingroup$ Yes there is an error, the sixth term should be q1p1 instead of q1 p1, $\endgroup$
    – Tigran
    Commented Oct 14, 2016 at 16:08

1 Answer 1

1
$\begingroup$

Here is how I would proceed. First, multiply f[ ] by I in order to get real valued function. Next calculate Groebner base and eliminate vars p1 and q1 in terms of q3:

gb = GroebnerBasis[
  Rationalize[{I (0.9816 q1 p1 + 1.1822 I q2 p2 + 1.2514 I q3 p3 - 
        0.0889594 I p2^3 - 0.0970894 I p2 p3^2 - 0.219332 p1 p2 q1 + 
        0.256743 I p2^2 q2 - 0.0970894 I p3^2 q2 + 
        0.0256743 I p2 q2^2 - 0.0889594 I q2^3 - 0.194179 p2 p3 q3 - 
        0.194179 p3 q2 q3 + 0.0970894 I p2 q3^2 + 
        0.0970894 I q2 q3^2) == rez, q1 p1 == q3}, 0], {p3, p2, q2, 
   q3}, {q1, p1}, MonomialOrder -> EliminationOrder]

{889594 p2^3 + 970894 p2 p3^2 - 11822000 p2 q2 - 2567430 p2^2 q2 + 970894 p3^2 q2 - 256743 p2 q2^2 + 889594 q2^3 + 9816000 I q3 - 2193320 I p2 q3 - 12514000 p3 q3 - 1941790 I p2 p3 q3 - 1941790 I p3 q2 q3 - 970894 p2 q3^2 - 970894 q2 q3^2 - 10000000 rez}

Now substitute q3 and f values and solve for auxiliary rez.

new = (rez /. Solve[(gb[[1]] /. q3 -> 0) == 4/10, rez])

{(1/50000000)(-2 + 4447970 p2^3 + 4854470 p2 p3^2 - 59110000 p2 q2 - 12837150 p2^2 q2 + 4854470 p3^2 q2 - 1283715 p2 q2^2 + 4447970 q2^3)}

This is expression you can plot. For example, for slice p3==1

Plot3D[neweq[[1]] /. {p3 -> 1}, {p2, -10, 10}, {q2, -20, 20}]

PlotOutput

$\endgroup$
4
  • $\begingroup$ Thanks, although I made an error. The sixth term in my function $f$ is written as q1 q2 p2 while it should be q1p1 q2. How can I get a 2d plot out of it? If I change the code to Plot2d[...] it gives an error? $\endgroup$
    – Tigran
    Commented Oct 14, 2016 at 16:11
  • $\begingroup$ The same way just change first part to: gb = GroebnerBasis[ Rationalize[{I (0.9816 q1p1 + 1.1822 I q2 p2 + 1.2514 I q3 p3 - 0.0889594 I p2^3 - 0.0970894 I p2 p3^2 - 0.219332 q1p1*q2 + 0.256743 I p2^2 q2 - 0.0970894 I p3^2 q2 + 0.0256743 I p2 q2^2 - 0.0889594 I q2^3 - 0.194179 p2 p3 q3 - 0.194179 p3 q2 q3 + 0.0970894 I p2 q3^2 + 0.0970894 I q2 q3^2) == rez, q1p1 == q3}, 0], {p3, p2, q2, q3}, {q1p1}, MonomialOrder -> EliminationOrder] $\endgroup$
    – Acus
    Commented Oct 14, 2016 at 16:51
  • $\begingroup$ When I use your script all I get is an empty cube [img]s15.postimg.org/dny4jwmiz/Untitled.jpg[/img] If I do it correct i should get something like this [img]scielo.br/img/revistas/jbsms/v24n3/16753f5.gif[/img] $\endgroup$
    – Tigran
    Commented Oct 17, 2016 at 7:55
  • $\begingroup$ I changed neweq to new and I got the 3dplot, Is there a way to get this in a 2d plot like this scielo.br/img/revistas/jbsms/v24n3/16753f5.gif $\endgroup$
    – Tigran
    Commented Oct 17, 2016 at 8:24

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.