[关闭]
@Xc-liu 2016-05-03T12:38:40.000000Z 字数 5325 阅读 965

homework_10程序

  • 程序一
  1. x_0=0
  2. y_0=0
  3. v_x_0=0.1
  4. v_y_0=0.1
  5. dt=0.01
  6. end_t=100
  7. x=[]
  8. y=[]
  9. v_x=[]
  10. v_y=[]
  11. t=[]
  12. x.append(x_0)
  13. y.append(y_0)
  14. v_x.append(v_x_0)
  15. v_y.append(v_y_0)
  16. t.append(0)
  18. for i in range(int(end_t/dt)):
  19.     x.append(x[i]+v_x[i]*dt)
  20.     y.append(y[i]+v_y[i]*dt)
  21.     v_x.append(v_x[i])
  22.     v_y.append(v_y[i])
  23.     t.append(t[i]+dt)
  24.     #print abs(x[i+1])+abs(y[i+1])
  25.     if  abs(x[i+1])+abs(y[i+1])>1:
  26.         while abs(abs(x[i+1])+abs(y[i+1])-1)>0.00001:
  27.               X=(x[i]+x[i+1])/2
  28.               Y=(y[i]+y[i+1])/2
  29.               a=(2**0.5/2.0)*abs(X)/X
  30.               b=(2**0.5/2.0)*abs(Y)/Y 
  31.               v_x[i+1]=(1-2*a**2)*v_x[i]-2*a*b*v_y[i]
  32.               v_y[i+1]=(1-2*b**2)*v_y[i]-2*a*b*v_x[i]
  33.               #v_y[i+1]=-((2.0*a**3-2.0*a)/b)*v_x[i+1]-(2*a**2-1)*v_y[i+1]
  34.               if abs(X)+abs(Y)>1:
  35.                  x[i+1]=X
  36.                  y[i+1]=Y
  37.                  continue
  38.               else:
  39.                  x[i]=X
  40.                  y[i]=Y
  41.                  continue
  42.         #if abs(abs(X)+abs(Y)-1)<0.1:
  43.           # break
  46. x_1=[-1]
  47. y_1=[0]
  48. x_2=[-1]
  49. y_2=[0]
  50. dx=0.1
  51. for j in range(20):
  52.     x_1.append(x_1[j]+dx)
  53.     y_1.append(1-abs(x_1[j+1]))
  54.     x_2.append(x_1[j]+dx)
  55.     y_2.append(-(1-abs(x_1[j+1])))
  58. import matplotlib.pyplot as plt
  59. import numpy as np
  60. import math
  61. plt.plot(x,y,'--' ,label='aa')
  62. plt.plot(x_1,y_1,'--',label='bianary')
  63. plt.plot(x_2,y_2,'--',label='bianary')
  64. plt.title(u'',fontsize=14)
  65. plt.xlabel(u'x',fontsize=14)
  66. plt.ylabel(u'y',fontsize=14)
  67. plt.legend(fontsize=12,loc='best')
  68. plt.show()
  • 程序二
  1. import math
  2. x_0=0
  3. y_0=0
  4. v_x_0=1
  5. v_y_0=0.002
  6. dt=0.01
  7. end_t=7000
  8. x=[]
  9. y=[]
  10. v_x=[]
  11. v_y=[]
  12. t=[]
  13. x.append(x_0)
  14. y.append(y_0)
  15. v_x.append(v_x_0)
  16. v_y.append(v_y_0)
  17. t.append(0)
  18. alpha=0.001
  19. v_x_1=[]
  20. x_1_1=[]
  23. for i in range(int(end_t/dt)):
  24.     x.append(x[i]+v_x[i]*dt)
  25.     y.append(y[i]+v_y[i]*dt)
  26.     v_x.append(v_x[i])
  27.     v_y.append(v_y[i])
  28.     t.append(t[i]+dt)
  29.     if abs(x[i+1]-1)<0.0001:
  30.        print x[i+1]
  31.        a=-x[i+1]
  32.        b=0
  33.        v_x[i+1]=(1-2*a**2)*v_x[i]-2*a*b*v_y[i]
  34.        v_y[i+1]=(1-2*b**2)*v_y[i]-2*a*b*v_x[i]
  35.     else:
  36.        if y[i+1]<0:
  37.           if  x[i+1]**2.0+(y[i+1]+alpha/2.0)**2.0>1:
  38.            while abs(x[i+1]**2+(y[i+1]+alpha/2.0)**2-1)>0.0001:
  39.               #print abs(x[i+1]**2+y[i+1]**2-1)              
  40.               X=(x[i]+x[i+1])/2
  41.               Y=(y[i]+y[i+1])/2
  42.               a=X/((X**2+(Y+alpha/2.0)**2)**0.5)
  43.               b=(Y+alpha/2.0)/((X**2+(Y+alpha/2.0)**2)**0.5)
  44.               v_x[i+1]=(1-2*a**2)*v_x[i]-2*a*b*v_y[i]
  45.               v_y[i+1]=(1-2*b**2)*v_y[i]-2*a*b*v_x[i]
  46.               #v_y[i+1]=-((2.0*a**3-2.0*a)/b)*v_x[i+1]-(2*a**2-1)*v_y[i+1]
  47.               if X**2+(Y+alpha/2.0)**2>1:
  48.                  x[i+1]=X
  49.                  y[i+1]=Y
  50.                  continue
  51.               else:
  52.                  x[i]=X
  53.                  y[i]=Y
  54.                  continue
  55.        else:
  56.           if  x[i+1]**2.0+(y[i+1]-alpha/2.0)**2.0>1:
  57.            while abs(x[i+1]**2+(y[i+1]-alpha/2.0)**2-1)>0.0001:
  58.               #print abs(x[i+1]**2+y[i+1]**2-1)
  59.               X=(x[i]+x[i+1])/2
  60.               Y=(y[i]+y[i+1])/2
  61.               a=X/((X**2+(Y-alpha/2.0)**2)**0.5)
  62.               b=(Y-alpha/2.0)/((X**2+(Y-alpha/2.0)**2)**0.5)
  63.               v_x[i+1]=(1-2*a**2)*v_x[i]-2*a*b*v_y[i]
  64.               v_y[i+1]=(1-2*b**2)*v_y[i]-2*a*b*v_x[i]
  65.               #v_y[i+1]=-((2.0*a**3-2.0*a)/b)*v_x[i+1]-(2*a**2-1)*v_y[i+1]
  66.               if X**2+(Y-alpha/2.0)**2>1:
  67.                  x[i+1]=X
  68.                  y[i+1]=Y
  69.                  continue
  70.               else:
  71.                  x[i]=X
  72.                  y[i]=Y
  73.                  continue
  75.     if abs(y[i]-0)<0.0001:
  76.        p=x[i]
  77.        q=v_x[i]
  78.        v_x_1.append(q)
  79.        x_1_1.append(p)
  81. print v_x[30]
  82. x_1=[-1]
  83. y_1=[0]
  84. x_2=[-1]
  85. y_2=[0]
  86. dx=0.01
  87. for j in range(199):
  88.     x_1.append(x_1[j]+dx)
  89.     y_1.append(alpha/2.0+(1-x_1[j+1]**2)**0.5)
  90.     x_2.append(x_2[j]+dx)
  91.     y_2.append(-alpha/2.0-(1-x_2[j+1]**2)**0.5)
  94. import matplotlib.pyplot as plt
  95. import numpy as np
  96. import math
  97. fig= plt.figure(figsize=(7,6))
  98. #fig= plt.figure(figsize=(7,6))
  99. plt.scatter(x_1_1,v_x_1,marker='+',color='b',label='phase space plot')
  100. #plt.plot(x,y,'--' ,label='orbit')
  101. #plt.plot(x_1,y_1,'--',label='bianary')
  102. #plt.plot(x_2,y_2,'--',label='bianary')
  103. plt.title(u'$v_x/v_y=1/0.002,y_0=0,x_0=0,a=0.001,scan-point-y=0$',fontsize=16) 
  104. plt.xlabel(u'x',fontsize=14)
  105. plt.ylabel(u'y',fontsize=14)
  106. plt.legend(fontsize=12,loc='best')
  107. plt.show(fig)
  • 程序三
  1. import math
  2. x_0=1.2
  3. y_0=0
  4. v_x_0=1
  5. v_y_0=2
  6. dt=0.01
  7. end_t=25
  8. x=[]
  9. y=[]
  10. v_x=[]
  11. v_y=[]
  12. t=[]
  13. x.append(x_0)
  14. y.append(y_0)
  15. v_x.append(v_x_0)
  16. v_y.append(v_y_0)
  17. t.append(0)
  18. #v_x_1=[]
  19. #x_1_1=[]
  20. q=1.6*10**-19
  21. m=1*10**-19
  22. b=0.9
  26. for i in range(int(end_t/dt)):
  27.     x.append(x[i]+v_x[i]*dt)
  28.     y.append(y[i]+v_y[i]*dt)
  29.     #v_x.append(v_x[i])
  30.     #v_y.append(v_y[i])
  31.     v_x.append(v_x[i]+q*b*v_y[i]*dt/m)
  32.     v_y.append(v_y[i]-q*b*v_x[i]*dt/m)
  33.     t.append(t[i]+dt)
  34.     #print x[i]**2+y[i]**2
  35.     if  x[i+1]**2+y[i+1]**2>4:
  36.         print  x[i+1]**2+y[i+1]**2
  37.         while abs(x[i+1]**2+y[i+1]**2-4)>0.001:
  38.               #print abs(x[i+1]**2+y[i+1]**2-4) 
  39.               X=(x[i]+x[i+1])/2
  40.               Y=(y[i]+y[i+1])/2
  41.               a=X/(X**2+Y**2)**0.5
  42.               b=Y/(X**2+Y**2)**0.5
  43.               v_x[i+1]=(1-2*a**2)*v_x[i]-2*a*b*v_y[i]
  44.               v_y[i+1]=(1-2*b**2)*v_y[i]-2*a*b*v_x[i]
  45.               #v_y[i+1]=-((2.0*a**3-2.0*a)/b)*v_x[i+1]-(2*a**2-1)*v_y[i+1]
  46.               if X**2+Y**2>4:
  47.                  x[i+1]=X
  48.                  y[i+1]=Y
  49.                  continue
  50.               else:
  51.                  x[i]=X
  52.                  y[i]=Y
  53.                  continue
  54.         #if abs(abs(X)+abs(Y)-1)<0.1:
  55.           # break
  57.     if  x[i+1]**2+y[i+1]**2<1:
  58.         while abs(x[i+1]**2+y[i+1]**2-1)>0.001:
  59.               X=(x[i]+x[i+1])/2
  60.               Y=(y[i]+y[i+1])/2
  61.               a=X/(X**2+Y**2)**0.5
  62.               b=Y/(X**2+Y**2)**0.5
  63.               v_x[i+1]=(1-2*a**2)*v_x[i]-2*a*b*v_y[i]
  64.               v_y[i+1]=(1-2*b**2)*v_y[i]-2*a*b*v_x[i]
  65.               #v_y[i+1]=-((2.0*a**3-2.0*a)/b)*v_x[i+1]-(2*a**2-1)*v_y[i+1]
  66.               if X**2+Y**2<1:
  67.                  x[i+1]=X
  68.                  y[i+1]=Y
  69.                  continue
  70.               else:
  71.                  x[i]=X
  72.                  y[i]=Y
  73.                  continue
  74.     #if abs(y[i])<0.0001:
  75.      #  p=x[i]
  76.      #  q=v_x[i]
  77.      #  v_x_1.append(q)
  78.      #  x_1_1.append(p)
  81. x_1=[-1]
  82. y_1=[0]
  83. x_2=[-1]
  84. y_2=[0]
  85. x_3=[-2]
  86. y_3=[0]
  87. x_4=[-2]
  88. y_4=[0]
  89. dx=0.0001
  90. for j in range(20000):
  91.     x_1.append(x_1[j]+dx)
  92.     y_1.append((1-x_1[j+1]**2)**0.5)
  93.     x_2.append(x_1[j]+dx)
  94.     y_2.append(-(1-x_1[j+1]**2)**0.5)
  96. for k in range(40000):
  97.     x_3.append(x_3[k]+dx)
  98.     y_3.append((4-x_3[k+1]**2)**0.5)
  99.     x_4.append(x_4[k]+dx)
  100.     y_4.append(-(4-x_4[k+1]**2)**0.5)
  102. #x_4=x_4[1:39998]
  103. #y_4=y_4[1:39998]
  105. import matplotlib.pyplot as plt
  106. import numpy as np
  107. import math
  108. fig= plt.figure(figsize=(6,6))
  109. #plt.scatter(x_1_1,v_x_1,marker='+',color='b',label='phase space plot')
  110. plt.plot(x,y,'--',label='orbit')
  111. plt.plot(x_1,y_1,'--',label='binary')
  112. plt.plot(x_2,y_2,'--',label='binary')
  113. plt.plot(x_3,y_3,'--',label='binary')
  114. plt.plot(x_4,y_4,'--',label='binary')
  115. plt.title(u'$v_x/v_y=1/2,x_0=1.2,y_0=0,magnitude=0.9T$',fontsize=14) 
  116. plt.xlabel(u'x',fontsize=14)
  117. plt.ylabel(u'y',fontsize=14)
  118. plt.legend(fontsize=12,loc='best')
  119. plt.show(fig)
添加新批注
在作者公开此批注前，只有你和作者可见。
回复批注