最简单而引人注目的混沌，莫过于三体运动。仅仅三颗星体的运动，就能变得复杂而眩目。这种复杂曾令数学家们在百年间困惑不已。如果只有两个天体，那么一切是多么简单，18世纪的伯努利就已解出了运动的所有可能轨迹，用合适的坐标，就能用简单的曲线描述。但仅仅是多了一个天体，就要等到19世纪的庞加莱，才给出了差强人意的答案：没有漂亮的解（正式术语是三体系统是不可积的）。这并非因为人类的智慧所限，而是从本质上来说，三个天体之间的运动轨迹不可能用简单的式子表达。自然并不像原来期盼的那么简单，它的复杂性令人绝望。小说《三体》中的三体人，大概最明白这一点。
下面是一个三体运动的GIF动画:
这里使用我定义语法的脚本代码生成混沌图像.相关软件参见:YChaos生成混沌图像.如果你对数学生成图形图像感兴趣,欢迎加入QQ交流群: 367752815
想要实现三体,先要了解下万有引力下的行星运动.比如地球如何绕太阳在椭圆型的轨道上运动:
[ScriptLines] k = [static]0.5*q*(u*u + v*v) p = [static]-g*w*q/sqrt(x*x + y*y) e = [static]k+p r = x*x + y*y d = sqrt(r) a = g*w/r i = -a*x/d j = -a*y/d x = x + u*t + 0.5*i*t*t y = y + v*t + 0.5*j*t*t u = u + i*t v = v + j*t p = -g*w*q/sqrt(x*x + y*y) k = e - p s = sqrt(2*k/q)/sqrt(u*u + v*v) u = u*s v = v*s [Variables] g=100.000000 q=0.100000 t=0.001000 u=4.000000 v=0.000000 w=1.000000 x=5.000000 y=8.000000
上述代码可以生成如下图像:
图像中,紫色的区域为近日端,绿色的区域为远日端.代码中有若干个参数,其中uv代表速度,q为质量,xy表示位置,t表示间隔时间。如果将t设置大一些,会得到一个圆形的轨道:
讲完基础，接下来是三体图像的脚本生成代码。在我之前的博客中，已有若干篇关于三体实现的程序介绍。这里实现的三体是最简单的一种情况：假定三体中的两个星体位置不动，只有一个可动星体在万有引力下的运动轨迹。并且既然是生成图像,所以对其的运算处理只在二维平面上。就算有了很大的简化，其脚本代码也十分复杂，这是我写的所有混沌图像生成脚本中最复杂的一个。写它的时候相当纠结，因为我总记不清哪个变量代表哪个意思，调试了好久才搞对。
代码如下:
[ScriptLines] k=[static]0.5*w*(u*u + v*v) p=[static]-g*m*w / sqrt((a - x)^2 + (b - y)^2) q=[static]-g*n*w / sqrt((c - x)^2 + (d - y)^2) e=[static]k + p + q f=(a - x)^2 + (b - y)^2 h=(c - x)^2 + (d - y)^2 r=g*m/f s=g*n/h i=r*(a - x)/sqrt(f) + s*(c - x)/sqrt(h) j=r*(b - y)/sqrt(f) + s*(d - y)/sqrt(h) x=x + u*t + 0.5*i*t*t y=y + v*t + 0.5*j*t*t u=u + i*t v=v + j*t f=(a - x)^2 + (b - y)^2 h=(c - x)^2 + (d - y)^2 p=-g*m*w/sqrt(f) q=-g*n*w/sqrt(h) k=e - p - q k=if(k<0, 0, k) s=sqrt(2*k/w)/sqrt(u*u + v*v) u=u*s v=v*s [Variables] a=-10.000000 b=0.000000 c=10.000000 d=0.000000 g=50.000000 m=1.000000 n=1.000000 t=0.000100 u=0.000000 v=2.000000 w=0.100000 x=5.000000 y=5.000000
代码中有若干个参数,其含义如下:
g表示万有引力系数
a,b表示固定星体1的位置
c,d表示固定星体2的位置
m表示固定星体1的质量
n表示固定星体2的质量w表示运动星体的质量
u,v表示运动星体的速度
x,y表示运动星体的位置
t表示间隔时间
代码有些难懂,还是请大家欣赏图像吧:
这种情况下的三体运动就是个纠结的过程,让我想到成语:朝三暮四,朝秦暮楚,东成西就
如果情况在复杂点呢.比如一个运动质点,在六个固定质点下的运动.这给我感觉像是"六道轮回"
其脚本十分难写:
[ScriptLines] p0x=[static]r*sin(0.0) p1x=[static]r*sin(PI/3) p2x=[static]r*sin(PI*2/3) p3x=[static]r*sin(PI*3/3) p4x=[static]r*sin(PI*4/3) p5x=[static]r*sin(PI*5/3) p0y=[static]r*cos(0.0) p1y=[static]r*cos(PI/3) p2y=[static]r*cos(PI*2/3) p3y=[static]r*cos(PI*3/3) p4y=[static]r*cos(PI*4/3) p5y=[static]r*cos(PI*5/3) m0=[static]1.0 m1=[static]1.0 m2=[static]1.0 m3=[static]1.0 m4=[static]1.0 m5=[static]1.0 k=[static]0.5*w*(u*u + v*v) d0x=[static]p0x - x d0y=[static]p0y - y q0=[static]-g*m0*w/sqrt(d0x*d0x + d0y*d0y) d1x=[static]p1x - x d1y=[static]p1y - y q1=[static]-g*m1*w/sqrt(d1x*d1x + d1y*d1y) d2x=[static]p2x - x d2y=[static]p2y - y q2=[static]-g*m2*w/sqrt(d2x*d2x + d2y*d2y) d3x=[static]p3x - x d3y=[static]p3y - y q3=[static]-g*m3*w/sqrt(d3x*d3x + d3y*d3y) d4x=[static]p4x - x d4y=[static]p4y - y q4=[static]-g*m4*w/sqrt(d4x*d4x + d4y*d4y) d5x=[static]p5x - x d5y=[static]p5y - y q5=[static]-g*m5*w/sqrt(d5x*d5x + d5y*d5y) e=[static]k + q0 + q1 + q2 + q3 + q4 + q5 f0=(p0x - x)*(p0x - x) + (p0y - y)*(p0y - y) r0=g*m0/f0 f1=(p1x - x)*(p1x - x) + (p1y - y)*(p1y - y) r1=g*m1/f1 f2=(p2x - x)*(p2x - x) + (p2y - y)*(p2y - y) r2=g*m2/f2 f3=(p3x - x)*(p3x - x) + (p3y - y)*(p3y - y) r3=g*m3/f3 f4=(p4x - x)*(p4x - x) + (p4y - y)*(p4y - y) r4=g*m4/f4 f5=(p5x - x)*(p5x - x) + (p5y - y)*(p5y - y) r5=g*m5/f5 i=r0*(p0x - x)/sqrt(f0) j=r0*(p0y - y)/sqrt(f0) i=i + r1*(p1x - x)/sqrt(f1) j=j + r1*(p1y - y)/sqrt(f1) i=i + r2*(p2x - x)/sqrt(f2) j=j + r2*(p2y - y)/sqrt(f2) i=i + r3*(p3x - x)/sqrt(f3) j=j + r3*(p3y - y)/sqrt(f3) i=i + r4*(p4x - x)/sqrt(f4) j=j + r4*(p4y - y)/sqrt(f4) i=i + r5*(p5x - x)/sqrt(f5) j=j + r5*(p5y - y)/sqrt(f5) x=x + u*t + 0.5*i*t*t y=y + v*t + 0.5*j*t*t u=u + i*t v=v + j*t d0x=p0x - x d0y=p0y - y q0=-g*m0*w/sqrt(d0x*d0x + d0y*d0y) d1x=p1x - x d1y=p1y - y q1=-g*m1*w/sqrt(d1x*d1x + d1y*d1y) d2x=p2x - x d2y=p2y - y q2=-g*m2*w/sqrt(d2x*d2x + d2y*d2y) d3x=p3x - x d3y=p3y - y q3=-g*m3*w/sqrt(d3x*d3x + d3y*d3y) d4x=p4x - x d4y=p4y - y q4=-g*m4*w/sqrt(d4x*d4x + d4y*d4y) d5x=p5x - x d5y=p5y - y q5=-g*m5*w/sqrt(d5x*d5x + d5y*d5y) k=e - q0 - q1 - q2 - q3 - q4 - q5 k=if(k < 0.0, 0, k) s=sqrt(2*k/w)/sqrt(u*u + v*v) u=u*s v=v*s [Variables] g=50.000000 r=1.000000 t=0.000050 u=0.000000 v=0.000000 w=0.100000 x=0.100000 y=-0.500000
相关文章: