Science 版 (精华区)

发信人: qpcwth (独翅鸟), 信区: Science
标  题: 《分形艺术》64
发信站: 哈工大紫丁香 (2001年11月03日18:29:17 星期六), 站内信件

第九章 微分方程系统
9.4若斯勒混沌
    1976年若斯勒对洛仑兹模型进行再建模，化简出一个新的方程组，仍然能够产生混
沌。洛仑 兹混沌吸引子是对称的，而若斯勒混沌吸引子丧失了对称性，整体结构颇像单
侧曲面——麦 比乌斯(A.F.Mobius,1790-1868)带。
    若斯勒方程的具体形式为
dx/dt=-y-z,
dy/dt=x+ay,
dz/dt=b+z(x-c)，
其中a,b和c都是正的参数，方程中只含有一个非线性项z x，而洛仑兹方程含有两个非线
性项。从形式上看若斯勒方程的确比洛仑兹方程简单。 显然，这也是一个自治的常微分
方程，右端不显含时间。
    先对系统作定性分析，将原系统划分为两个相互关联的子系统，第一个子系统由前
两个方程 代表，第二个子系统由第三个方程代表。当z足够小时，第一个方程中的z可 
以暂时忽略，这时第一个子系统可以暂时简化为
dx/dt=-y,
dy/dt=x+ay,
    它实际上是一个二阶线性振子，将第一个方程代入第二个方程立即得到如下微分方
程：
d^2x/dt^2-a(dx/dt)+x=0.
    学过振动理论的对此方程一定是非常熟悉的，当a为正数时，它代表负阻尼振子的振
 动情况。在相平面(x,y)上考虑问题，原点(0，0)是一个不稳定焦点(focus) ，轨线从
原点附近向外盘旋，圈越转越大。大到一定程度，整个系统的非线性(在第三个方 程中
体现出来)就必须考虑了。非线性起什么作用呢?非线性限制了轨线无限制向外盘旋。如
 果没有非线性，轨线很快奔向无穷远处。非线性项的存在阻止了这一点。我们看第三个
方程 ，当x大于c时，z的系数变正，取b为正值，这时第三个方程 所代表的子系统变得
不稳定。
    非线性项起作用的结果是z值增大，而z增大的结果又会导致x减小， 因为z增大后，
第一个方程中右端变小，以至于成为负值，于是x减小。x 逐步变小，以至于小于c，这
样又会使z变小，相轨道逐渐落入( x,y)平面，并接近原点。这时第三个方程又是次要的
了，第一个子系统 起支配作用，轨道又向外盘旋，再升高，再收缩，再降低，等等。定
性上分析，轨道运动大 概就是这个样子：在(x,y)平面上拉伸，在z轴方向折叠。两个子
系 统奇妙地耦合在一起(通过非线性)，子系统的稳定性交替变化，相互制约。
    若斯勒系统的数值积分程序与洛仑兹系统的完全一样，没必要单独写出来。不过，
现在要考 虑另一件重要的事情：吸引子的任意投影。在此之前我们只能做XOY，XOZ，Y
OZ三个 特殊方向上的投影，而这是很不够的。在计算三维微分方程系统时，得到三维数
据(x ,y,z)，任意角度投影的思想就是将这三维数据压缩成两维，并使这两 维的数据仍
然包含原来三维数据的信息，也就是说想找到一种变换，使(x,y ,z)→(X_N,Y_N)，只要
指定两个投影角度，就能顺利得出变换公式。设投影 角度分别是θ和φ，则投影算法为
：
X_N=xcosθ-ysinθ,
Y_N=xsinθsinφ+ycosθsinφ+ zcosφ,
Z_N=xsinθcosφ+ycosθcosφ- zsinφ.
    如果不考虑透视关系，只用前面的两式就可以了。如果考虑透视效果，还需用第三
式。我们 这里只用前两式。现在给出涉及任意角度投影的计算程序：
{Rossler.PAS, Huajie,1993}
Program RossAttrwithProj;
Uses Graph,Crt;
const
CenX=300;CenY=400; {origin coordinates in my screen}
multip=20;
Procedure Projt(var x,y,z,p,q,xN,yN,zN:real);
begin
xN:=X*cos(p)-Y*sin(p);    {theta=p, phi=q}
yN:=X*sin(p)*sin(q)+Y*cos(p)*sin(q)+Z*cos(q);
zN:=cos(q)*(X*sin(p)+Y*cos(p))-Z*sin(q);
end;
var
gd,gm，lx,ly,lz,n:integer;
x,y,z，a,b,c,delta,u,v,xN,yN,zN:real;
kx1,kx2,kx3,kx4,ky1,ky2,ky3,ky4,kz1,kz2,kz3,kz4:real;
begin
sound(500);delay(200);nosound;
u:=20;v:=45;{Two angles to look (project) the figure:
u=theta presents the angle with YOZ plane,
v=phi presents the angle with XOY plane.    }
u:=pi*u/180.0;v:=pi*v/180.0;
gd:=vga; gm:=vgahi;initgraph(gd,gm,’D:\\PASCAL’);
if graphresult <> grok then halt(1);
x:=15;y:=0;z:=0; {drawing X axis}
Projt(x,y,z,u,v,xN,yN,zN);
lx:=round(xN*multip);ly:=round(yN*multip);
Line(CenX,CenY,CenX+lx,CenY-ly);
OutTextXY(CenX+lx+4,CenY-ly+2,’X’);
x:=0;y:=15;z:=0; {drawing Y axis}
Projt(x,y,z,u,v,xN,yN,zN);
lx:=round(xN*multip);ly:=round(yN*multip);
Line(CenX,CenY,CenX+lx,CenY-ly);
OutTextXY(CenX+lx-10,CenY-ly+2,’Y’);
x:=0;y:=0;z:=15; {drawing Z axis}
Projt(x,y,z,u,v,xN,yN,zN);
lx:=round(xN*multip);
ly:=round(yN*multip);
Line(CenX,CenY,CenX+lx,CenY-ly);
OutTextXY(CenX+lx+4,CenY-ly+2,’Z’);
x:=2;y:=2;z:=2;  {initial point}
a:=0.40;b:=0.3;c:=4.5; {parameters}
delta:=0.002; {integraion step}
repeat
kx1:=delta*(-y-z);
kx2:=kx1;kx3:=kx1;kx4:=kx1;
ky1:=delta*(x+a*y);
ky2:=delta*(x+a*(y+ky1/2));
ky3:=delta*(x+a*(y+ky2/2));
ky4:=delta*(x+a*(y+ky3));
kz1:=delta*(b-c*z+x*z);
kz2:=delta*(b+(x-c)*(z+kz1/2));
kz3:=delta*(b+(x-c)*(z+kz2/2));
kz4:=delta*(b+(x-c)*(z+kz3));
x:=x+(kx1+2*kx2+2*kx3+kx4)/6;
y:=y+(ky1+2*ky2+2*ky3+ky4)/6;
z:=z+(kz1+2*kz2+2*kz3+kz4)/6;
Projt(x,y,z,u,v,xN,yN,zN);
lx:=round(xN*multip);
ly:=round(yN*multip);
 {lz:=round(zN*multip);}
PutPixel(lx+CenX,CenY-ly,15); {PutPixel(lx+CenX,CenY-lz,15);}
until KeyPressed;
sound(500);delay(200);nosound;
readln;CloseGraph;
end.

--
心事浩茫连广宇,于无声处听惊雷

※ 来源:·哈工大紫丁香 bbs.hit.edu.cn·[FROM: 202.118.229.154]