一个小球碰撞的模拟

qibbxxt

最近看本老外的书,上面提到一个例子:一个球静止躺在一个方形的桌上,球被撞击一下,然后无限的运动,问在什么情况下,小球可以回到原来的点,老外给出的方法就是把平面区域展开的方法,给了一个例子,但是我模拟了下,感觉不对(按老外的观点NN次就可以回到原点,但是模拟的结果不是,我运行了2*NN次，大致回到那个位置，但是好像不在那个点，不知是不时计算精度的问题),以下是我的程序,我愿意和有兴趣的读者一起讨论,

1. 
   
2. clear; clc; close all
   
3. step = 0.01;
   
4. N=1000;
   
5. xdir=1;
   
6. ydir=1;
   
7. xlow=1;
   
8. xup=2;
   
9. hold on
   
10. grid on
    
11. set(gca,'GridLineStyle','-');
    
12. daspect([1,1,1]);
    
13. % [x,y] = meshgrid((xlow+.1):step:(xup-.1));
    
14. % pos=unidrnd(numel(x));
    
15. % x0=x(pos);
    
16. % y0=y(pos);
    
17. x0=1.2;
    
18. y0=1.5;
    
19. xs=x0;ys=y0;
    
20. x2=x0;
    
21. y2=y0;
    
22. tanfz=3;
    
23. tanfm=7;
    
24. % tanfz=randint(1,1,[1,10]);
    
25. % tanfm=randint(1,1,[1,10]);
    
26. tany=tanfz/tanfm;
    
27. axis([0 16 0 8]);
    
28. plot(x0,y0,'go','markersize',5);
    
29. h=plot(x0,y0,'r.','markersize',20);
    
30. h1=plot(x0,y0,'g.','markersize',20);
    
31. set(gca,'xtick',0:1:16,'ytick',0:1:8);
    
32. gc=gcd(tanfz,tanfm);
    
33. NN=tanfm/gc/step;
    
34. for i=1:2*NN
    
35.     x1=x0;
    
36.     y1=y0;
    
37.     x3=x2;
    
38.     y3=y2;
    
39.     x0=x0+xdir*step;
    
40.     y0=y0+ydir*step*tany;
    
41.     x2=x2+step;
    
42.     y2=y2+step*tany;
    
43.     if x0>=xup || x0<=xlow
    
44.         xdir = -xdir;
    
45.     end
    
46.     if y0>=xup || y0<=xlow
    
47.         ydir = -ydir;
    
48.     end
    
49.     set(h,'xdata',x0,'ydata',y0);
    
50.     set(h1,'xdata',x2,'ydata',y2);
    
51.     plot([x1,x0],[y1,y0],'k-','linewidth',1);
    
52.     plot([x3,x2],[y3,y2],'y--','linewidth',1);
    
53.     pause(0.001);
    
54.     if dist([xs,ys],[x0;y0])<eps
    
55.         break
    
56.     end
    
57. end
    

复制代码