最簡單而引人注目的混沌，莫過於三體運動。僅僅三顆星體的運動，就能變得復雜而眩目。這種復雜曾令數學家們在百年間困惑不已。如果只有兩個天體，那么一切是多么簡單，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
相關文章: