轮胎型相图求助
曾经在一篇论文上看到过一种相图，作者将它叫做轮胎型相图，目前自己做的图 也是这样，但是要追究原因，但是又一时找不到那篇论文了，请教各位有没有人看过类似的论文，或者帮我分析一下这类相图的原因 如下图形式，系统存在阻尼，但是速度和位移随时间不会为零 :@L :@L 没人理我，这里面是不是存在内部激励呢？可是没有外加激励 我遇到过，但是是我的参数错了之后得到的，这个有可能是概周期响应，所以你试一下Poincare图看一下是什么样子的？ 呵呵，自治系统的，不知道怎么画它的Poincare映射呢目前不过我肯定见过这类图象，即便是我自己的算错了
回复 #5 咕噜噜 的帖子
小咕，在maple里面有个函数，好像就是poincare，可以直接画Hamilton系统的Poincare截面，你可以把你的系统写成Hamilton函数，然后用这个函数试试！ ^_^，maple不熟悉啊，用惯了matlab，有时间看看回复 #7 咕噜噜 的帖子
maple里面确实有，但是好像首先要化为Hamiton系统的形式，然后才能够求解出来 为方便起见，把这个贴出来，小咕可以直接试试，应该不难的！有本《混沌与Melnikov方法》有一章专讲Hamilton系统的，可以参考参考
DEtools - 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 , [...],..., where is a list of numbers representing the initial values for .
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(..).
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= variables constituting the 2-D plane of
the phase space where the 2PS is plotted
scene= variables constituting the 2-D plane and
related ranges for the 2PS plot
scene= 2-D plane for plotting the 2PS and a third
variable, xk, the cross-variable or the
third axis in 3PS
scene= 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= or scene=, with or without extra ranges (for displaying a plot of a particular region), as in scene=. 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= 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, {},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, {},
stepsize=.1,iterations=3,scene=);
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, {},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+ q1q2 - - 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 do
ics := 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 do
F4 := poincare(H,t=-150..150,ics,
stepsize=.1,iterations=2,scene=);
end do:
FF4 := array([,F4,F4],,F4,F4]]):
plots(FF4);
回复 #9 octopussheng 的帖子
好的，我一定仔细看看，琢磨一下，多谢，但是这个化简一般都能做到吗，会不会有一些系统根本就不能化简为Hamiton系统的形式回复 #10 咕噜噜 的帖子
那肯定有一些是非ｈａｍｉｔｉｏｎ系统就不能够化了关于Hamliton系统
对平面系统：dx/dt=f(x), x=', f=' ————————(1)设以原点位简单中心
定义：平面某区域D上所定义的实值函数P(x1,x2),为r+1阶连续可导函数，称为上述系统(1)
若沿着系统(1)解曲线，P(x1,x2)等于常数，又，如果P(x1,x2)存在，则称系统(1)可积。
如果域D内不含充满螺线的区域，以及没有极限的分界线，因此(1)在D内存在全局积分，从而具有全局的积分因之m(x1,x2)，使得方程
f1(x1,x2)dx2-f2(x1,x2)dx1=0
在乘以m(x1,x2)以后，化为恰当的微分方程，由于(mf1)x1+(mf2)x2（这里的x1和x2都是下标）等于0,因此（1）右端若乘以m(x1,x2)，即可化为全微分方程，又称Hamilton系统
这个就是Hamilton系统的定义！
如果还有别的定义、说法，请补充 真实的，耗散可忽略不计的物理过程才可以表示成Hamilton体系，
系统存在阻尼是不能表示为标准Hamilton系统的吧？
我觉得这种相图要从你选用的数值求解微分方程算法本身的特性上去考虑 请问给出了一个Hamiltonian 系统，怎么写出它的Hamiltonian function？ 一个简单的例子：Duffing振子：x''+x-x^3=-e*d*x'
式中，e为小参数。
未扰系统为x''+x-x^3=0
———> x'=y
y'=-x+x^3
则其Hamliton函数可写成：H(x,y)=1/2*y^2+1/2*x^2-1/4*x^4=h