利用Python研究平方反比斥力场中粒子的运动,以 alpha 粒子在重核场中的运动为例。
Astro_Lee
编辑于 2021年02月04日 22:10

彭芳麟《计算物理基础》第五章练习第2题

研究平方反比斥力场中粒子的运动。以%5Calpha粒子在重核场中的运动为例,设重核位于力心且固定不动,%5Calpha粒子的质量为m,它到重核的距离为r,所受到库仑斥力F%3Dk%2Fr%5E2k为由库仑定律确定的常量。要求:

本题微分方程

%5Cbegin%7Beqnarray%7D%0A%5Cleft%5C%7B%5Cbegin%7Bmatrix%7D%20%0A%5Cfrac%7Bd%5E2r%7D%7Bdt%5E2%7D%20%26%20%3D%20%26%20r%20%5Cleft(%5Cfrac%7Bd%5Ctheta%7D%7Bdt%7D%5Cright)%5E2%2B%5Cfrac%7Bk%7D%7Br%5E2%7D%20%5C%5C%20%20%0A%5Cfrac%7Bd%5E2%5Ctheta%7D%7Bdt%5E2%7D%20%26%20%3D%20%26%20-%5Cfrac%7B2%7D%7Br%7D%5Cfrac%7Bdr%7D%7Bdt%7D%5Cfrac%7Bd%5Ctheta%7D%7Bdt%7D%0A%5Cend%7Bmatrix%7D%5Cright.%0A%5Cend%7Beqnarray%7D


(1)画出%5Calpha粒子在不同初始条件下的轨道,通过改变初始条件研究影响散射角的因素。

r_1%3Dr%2C%5C%20r_2%3Dr_1%5E%7B%26%2339%3B%7D%3D%5Cfrac%7Bdr_1%7D%7Bdt%7D%2C%5C%20r_2%5E%7B%26%2339%3B%7D%3Dr_1%5E%7B%26%2339%3B%26%2339%3B%7D%3D%5Cfrac%7Bd%5E2r_1%7D%7Bdt%5E2%7D

%5Ctheta_1%3D%5Ctheta%2C%5C%20%5Ctheta_2%3D%5Ctheta_1%5E%7B%26%2339%3B%7D%3D%5Cfrac%7Bd%5Ctheta_1%7D%7Bdt%7D%2C%5C%20%5Ctheta_2%5E%7B%26%2339%3B%7D%3D%5Ctheta_1%5E%7B%26%2339%3B%26%2339%3B%7D%3D%5Cfrac%7Bd%5E2%5Ctheta_1%7D%7Bdt%5E2%7D

r_1%5E%7B%26%2339%3B%7D%3Dr_2%2C%5C%20r_2%5E%7B%26%2339%3B%7D%3Dr_1(%5Ctheta_2)%5E2%2B%5Cfrac%7Bk%7D%7Br_1%5E2%7D

%5Ctheta_1%5E%7B%26%2339%3B%7D%3D%5Ctheta_2%2C%5C%20%5Ctheta_2%5E%7B%26%2339%3B%7D%3D-%5Cfrac%7B2%7D%7Br_1%7Dr_2%5Ctheta_2


代码实现:

代码块
Python
自动换行
复制代码
# 导入所需模块
import numpy as np
import matplotlib.pyplot as plt
plt.rcParams['font.sans-serif']=['SimHei'] # 用来正常显示中文标签
plt.rcParams['axes.unicode_minus']=False # 用来正常显示负号
plt.rcParams['axes.linewidth'] = 1.5
plt.rcParams['lines.linewidth'] = 3
plt.rcParams['xtick.labelsize'] = 20
plt.rcParams['ytick.labelsize'] = 20
from scipy.integrate import odeint
复制成功
代码块
Python
自动换行
复制代码
#直角坐标系下的微分方程
def fun(y,t):
    k = 3
    return np.array([y[1],#x一阶导
                     k*y[0]/np.sqrt(y[0]**2+y[2]**2)**3,#x二阶导
                    y[3],#y一阶导
                     k*y[2]/np.sqrt(y[0]**2+y[2]**2)**3])#y二阶导

#初始条件
y0=np.array([[-10,1,10,0],
    [-10,1,-10,0],
    [-10,1,2,0],
    [-10,1,-2,0],
    [-10,2,2,0],
    [-10,2,-2,0],
    [-10,1,4,0],
    [-10,1,-4,0],
    [-10,1,15,0],
    [-10,1,-15,0],
    [-10,1,7,0],
    [-10,2,-7,0],
    [-10,2,7,0],
    [-10,1,-7,0],
    [-10,1,0,0.4],
    [-10,1,0,-0.4]])

t = np.arange(0,35,0.1)

plt.figure(figsize=(15,8),dpi=300)
for i in range(len(y0)):
    res = odeint(fun, y0[i], t)
    x,y = res[:,0],res[:,2]
    plt.plot(x,y)
    
plt.plot(0,0,'ko',ms=10)
plt.text(0.5,-0.4,'靶粒子',fontsize=20)
plt.axis([-10,20,-20,20])
复制成功


图像:

散射结果

(2)根据解决问题的需要来选择坐标系,本题就是选择直角坐标系而不是极坐标系。

如有错误,欢迎指正!