下面是一个三体运动的GIF动画:
这里使用我定义语法的脚本代码生成混沌图像.相关软件参见:.如果你对数学生成图形图像感兴趣,欢迎加入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+pr = x*x + y*yd = sqrt(r)a = g*w/ri = -a*x/dj = -a*y/dx = x + u*t + 0.5*i*t*ty = y + v*t + 0.5*j*t*tu = u + i*tv = v + j*tp = -g*w*q/sqrt(x*x + y*y)k = e - ps = sqrt(2*k/q)/sqrt(u*u + v*v)u = u*sv = v*s[Variables]g=100.000000q=0.100000t=0.001000u=4.000000v=0.000000w=1.000000x=5.000000y=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 + qf=(a - x)^2 + (b - y)^2h=(c - x)^2 + (d - y)^2r=g*m/fs=g*n/hi=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*ty=y + v*t + 0.5*j*t*tu=u + i*tv=v + j*tf=(a - x)^2 + (b - y)^2h=(c - x)^2 + (d - y)^2p=-g*m*w/sqrt(f)q=-g*n*w/sqrt(h)k=e - p - qk=if(k<0, 0, k)s=sqrt(2*k/w)/sqrt(u*u + v*v)u=u*sv=v*s[Variables]a=-10.000000b=0.000000c=10.000000d=0.000000g=50.000000m=1.000000n=1.000000t=0.000100u=0.000000v=2.000000w=0.100000x=5.000000y=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.0m1=[static]1.0m2=[static]1.0m3=[static]1.0m4=[static]1.0m5=[static]1.0k=[static]0.5*w*(u*u + v*v)d0x=[static]p0x - xd0y=[static]p0y - yq0=[static]-g*m0*w/sqrt(d0x*d0x + d0y*d0y)d1x=[static]p1x - xd1y=[static]p1y - yq1=[static]-g*m1*w/sqrt(d1x*d1x + d1y*d1y)d2x=[static]p2x - xd2y=[static]p2y - yq2=[static]-g*m2*w/sqrt(d2x*d2x + d2y*d2y)d3x=[static]p3x - xd3y=[static]p3y - yq3=[static]-g*m3*w/sqrt(d3x*d3x + d3y*d3y)d4x=[static]p4x - xd4y=[static]p4y - yq4=[static]-g*m4*w/sqrt(d4x*d4x + d4y*d4y)d5x=[static]p5x - xd5y=[static]p5y - yq5=[static]-g*m5*w/sqrt(d5x*d5x + d5y*d5y)e=[static]k + q0 + q1 + q2 + q3 + q4 + q5f0=(p0x - x)*(p0x - x) + (p0y - y)*(p0y - y)r0=g*m0/f0f1=(p1x - x)*(p1x - x) + (p1y - y)*(p1y - y)r1=g*m1/f1f2=(p2x - x)*(p2x - x) + (p2y - y)*(p2y - y)r2=g*m2/f2f3=(p3x - x)*(p3x - x) + (p3y - y)*(p3y - y)r3=g*m3/f3f4=(p4x - x)*(p4x - x) + (p4y - y)*(p4y - y)r4=g*m4/f4f5=(p5x - x)*(p5x - x) + (p5y - y)*(p5y - y)r5=g*m5/f5i=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*ty=y + v*t + 0.5*j*t*tu=u + i*tv=v + j*td0x=p0x - xd0y=p0y - yq0=-g*m0*w/sqrt(d0x*d0x + d0y*d0y)d1x=p1x - xd1y=p1y - yq1=-g*m1*w/sqrt(d1x*d1x + d1y*d1y)d2x=p2x - xd2y=p2y - yq2=-g*m2*w/sqrt(d2x*d2x + d2y*d2y)d3x=p3x - xd3y=p3y - yq3=-g*m3*w/sqrt(d3x*d3x + d3y*d3y)d4x=p4x - xd4y=p4y - yq4=-g*m4*w/sqrt(d4x*d4x + d4y*d4y)d5x=p5x - xd5y=p5y - yq5=-g*m5*w/sqrt(d5x*d5x + d5y*d5y)k=e - q0 - q1 - q2 - q3 - q4 - q5k=if(k < 0.0, 0, k)s=sqrt(2*k/w)/sqrt(u*u + v*v)u=u*sv=v*s[Variables]g=50.000000r=1.000000t=0.000050u=0.000000v=0.000000w=0.100000x=0.100000y=-0.500000
相关文章: