|
|
为方便起见，把这个贴出来，小咕可以直接试试，应该不难的！
有本《混沌与Melnikov方法》有一章专讲Hamilton系统的，可以参考参考
DEtools[poincare] - plot 2-D or 3-D projections of the Poincare surface-of-section (2PS/3PS) of a given Hamiltonian system
Calling Sequence
poincare(H, t=a..b, ics)
poincare(H, t=a..b, ics, options, 3)
Parameters
H - any algebraic expression representing the Hamiltonian
t=a..b - t represents the time and a..b is a numerical range
ics - set of initial conditions for the p's and q's
options - (optional) equations of the form keyword=value
3 - (optional) obtain a 3-D plot (3PS)
Description
The poincare command builds either fast but not so accurate, or slower, as accurate as desired, 2-D or 3-D projection plots of Poincare sections (see optional arguments below). Instead of analytically enforcing the Hamiltonian constraint, the value of the energy of each plotted point is checked against the corresponding H_0. As a complement, another routine, generate_ic, was programmed in order to speed up the preparation of initial conditions for the numerical experiments.
The returned plots can be manipulated using standard Maple facilities available on the icon tool bar, and using the DEtools,zoom command, also part of this package. These facilities usually permit a detailed visual distinction between the KAM surfaces and the layers of stochasticity for, say, typical weakly perturbed systems.
In addition to the returned plot, the following relevant information related to each solution curve is displayed on the screen during the calculation:
- the initial value of the Hamiltonian, p's, and q's for the curve
- the number of intersection points found in the given time interval (for 2-D plots)
- the maximum percentile "energy deviation" of the intersection points
It is useful to know the number of intersection points, since this can indicate the appropriateness of the indicated time interval. Concerning the percentile energy deviation, it is calculated as follows:
(H_0 - H(point))/H_0*100
When H_0=0, only the maximum absolute deviation is displayed. Note that all numerical algorithms lead to values of H different from the initial H_0, especially in the case of optional fast plottings. Percentile deviations below 10^(-8) are displayed as 0. Also, the deviation is calculated for each intersection point, but only the greatest value is displayed. This information gives an idea of the accuracy of the plot, and you can either re-enter the instruction looking for a more accurate but slower plot, or use your own numerical integration algorithm (see below). In typical situations, the smooth curves can be recognized even with energy deviations of the order of H_0/10.
THE ARGUMENTS
The first argument of poincare is the Hamiltonian. Some useful conventions were adopted to represent the p's and q's. All p's and q's must appear as pn or qn where n is a positive integer, as in p1, p2, and the time dependence need not be explicit, as in pn or qn instead of pn(t) or qn(t). The Hamiltonian is assumed to be time independent and the number of degrees of freedom is expected to be less than 10 (20-dimensional phase space), but this restriction can easily be removed.
The solution curves are calculated within a given time interval specified in the second argument as t=a..b. When this time range leads to no intersection points, you can use the 3 option to see how far the intersection plane is from the trajectories, and re-enter the instruction shifting the intersection plane using the shift option (see below).
The third argument, a set, may have any number of lists of initial conditions for the time, p's, and q's, corresponding or not to the same initial value of H; they can be generated by the user or by the generate_ic command. The initial conditions must be inside a set structure as in [n1,n2,n3,...], [...],..., where [n1,n2,n3,...] is a list of numbers representing the initial values for [t, p1,...pk, q1,...qk].
This function is part of the DEtools package, and so it can be used in the form poincare(..) only after executing the command with(DEtools). However, it can always be accessed through the long form of the command by using DEtools[poincare](..).
THE OPTIONAL ARGUMENTS
The optional equations, options, consist of the following.
'iterations'
'method'
'scene'
'shift'
'stepsize'
The optional arguments can be given alone or in conjunction and in any order.
iterations=positive integer
iterations is the number of iterations used in the integration scheme.
When requesting a plot, the numerical algorithm can be iterated as many times as desired by giving the extra argument, iterations = N. By default, iterations = 1.
method=proc
method is a user procedure to be used as the integration method.
The default numerical algorithm used in the integration scheme is basically the fourth order Runge-Kutta of the DEtools package, but this can be changed by using the method = proc option. Here, proc should be a numerical integration algorithm (see DEtools).
scene=list
There are four options for scene:
scene=[xi,xj] variables constituting the 2-D plane of
the phase space where the 2PS is plotted
scene=[xi=a..b,xj=c..d] variables constituting the 2-D plane and
related ranges for the 2PS plot
scene=[xi,xj,xk] 2-D plane for plotting the 2PS and a third
variable, xk, the cross-variable or the
third axis in 3PS
scene=[xi=a..b,xj=c..d,xk=e..f] same as above but including ranges for
the 2PS/3PS plot
The default scene for the plots is the (p1,q1) plane, at q2=0 for the 2PS, or the (p1,q1,q2) 3-D submanifold, for the plot of a 2PS embedded in a 3PS, when the 3 option is indicated. The intersection points constituting the 2PS are obtained by looking for the sign change of a third coordinate, denoted here by the cross variable, by default q2. The default ranges for the display of a plot will include all the calculated intersection points in the case of 2-D plots; or all calculated pieces of projections of solution curves plus the intersection 2-D plane, in the case of 3-D plots.
The 2-D or 3-D submanifolds or the cross variable can be changed by giving the extra argument scene=[x1,x2] or scene=[x1,x2,x3], with or without extra ranges (for displaying a plot of a particular region), as in scene=[x1=a..b,...]. The 2PS is then plotted over the plane formed by the first two variables appearing in the right hand side; and x3, when given, will be the cross variable, or the third axis when the 3 option is also given as an argument.
It is possible to indicate the time, t, as the third variable, in which case a convenient mouse manipulation of the 3-D plot can display the projection of the curves over each (qi,qj) plane. This may be useful in studying the bounded/unbounded properties of a given potential. Furthermore, in the case of a system with three degrees of freedom, the use of the 3 option with scene=[q1,q2,q3] will render the 3-D plot of the physical trajectory. When ranges are indicated, it is still possible to zoom in or out the resulting plot, up to the default ranges mentioned above, by using the DEtools,zoom command.
shift=a real number
shift is a positive/negative shift for the intersection plane in the plots of 2PS/3PS.
The intersection plane over which the 2PS is plotted can be shifted in a positive or negative direction, by indicating shift=s as an additional argument.
stepsize=positive number
stepsize is the size of the step used in the integration scheme.
By default, the step interval is (b-a)/20, where a..b is the range for t. As the stepsize is decreased, the accuracy and the smoothness of the integral curves (as well as the time consumed in the calculations) will increase.
Examples
1. The Toda lattice
with(DEtools):
H := 1/2*(p1^2 + p2^2) + 1/24*(exp(2*q2+2*sqrt(3)*q1)
+ exp(2*q2-2*sqrt(3)*q1) + exp(-4*q2))-1/8;
1 2 1 2 1 / (1/2) \ 1 / (1/2) \
H := - p1 + - p2 + -- exp\2 q2 + 2 3 q1/ + -- exp\2 q2 - 2 3 q1/
2 2 24 24
1 1
+ -- exp(-4 q2) - -
24 8
Create a 2PS over the q2=0 plane with 86 intersection points lying on smooth curves.
poincare(H,t=-150..150, {[0,.1,1.4,.1,0]},stepsize=.05,iterations=5);
Create a 2PS over the q1=0 plane with 98 intersection points. The smoothness of the curves in both (p,q) planes is related to the integrability of the system.
poincare(H,t=-100..100, {[0,.1,1.4,.1,0]},
stepsize=.1,iterations=3,scene=[p2,q2]);
Show the KAM surfaces of solution curves of regular motion. Interesting angles: theta=70, phi=130, Theta=90, Phi=180 and Theta=-20, Phi=75.
poincare(H,t=-100..100, {[0,.1,1.4,.1,0]},stepsize=.1,iterations=4,3);
2. The Henon-Heiles Hamiltonian
H := 1/2*(p1^2+p2^2+q1^2+q2^2)+q1^2*q2-q2^3/3;
1 2 1 2 1 2 1 2 2 1 3
H := - p1 + - p2 + - q1 + - q2 + q1 q2 - - q2
2 2 2 2 3
Generate the initial conditions for the numerical experiments by using generate_ic. (Here, we obtain six sets, related to each value of H_0 respectively, with three different initial conditions each.)
for h in [1/24,1/18,1/12,1/8,1/7,1/6] do
ics[h] := generate_ic(H,{t=0,p2=0.1,q2=-0.2..0.2,q1=-0.2..-0.1,energy=h},3)
end do:
Create the surfaces-of-section with around 300 points, and with percentile H-deviations around 10^(-4). This displays the progressive disintegration of KAM surfaces.
for h in [1/24,1/18,1/12,1/8,1/7,1/6] do
F4[h] := poincare(H,t=-150..150,ics[h],
stepsize=.1,iterations=2,scene=[p2=-.5..0.5,q2=-.5..0.5]);
end do:
FF4 := array([[F4[1/24],F4[1/18],F4[1/12]],[F4[1/8],F4[1/7],F4[1/6]]]):
plots[display](FF4); |
评分
-
1
查看全部评分
-
|
发表于 2007-7-29 16:00
提升卡
变色卡
显身卡
楼主
发表于 2007-7-30 12:06
