Python 越级轨迹代码答案

【问题标题】：Python leapfrog trajectory codePython 越级轨迹代码
【发布时间】：2018-02-16 23:06:59
【问题描述】：

我是一名物理专业的学生，但在编程方面是新手。几年前，我已经学会了如何为引力场中的粒子的跳跃式积分器编写C代码，但是记忆模糊，我现在正在做的是用Python编写粒子的跳跃式积分器的代码在一定的磁场中。有人告诉我 Boris 算法更适合这种模拟，但我决定先用我之前学到的东西进行实验，即越级积分器。但是 C 和 Python 的语法差异太大（至少对我而言），所以我不能只将 C 代码转换为 Python，我必须编写一个新代码。所以我不确定算法是否正确。

我的代码是这样的，

# -*- coding: utf-8 -*-
"""
Created on Thu Feb 15 19:00:55 2018

@author: Heptacle
"""

import numpy as np
import matplotlib.pyplot as plt 

q=1.6e-19 # unit charge
m=1.67e-27 # proton mass
xs=1 # x_star
B0=1 # maximum magnetic field
b=B0/xs

initial_position = np.array((0.1, 0,0))         # Initial position vector of Particle
initial_velocity = np.array((0, 0,0.1))           # Initial velocity vector of Particle

num_steps = 4000
time_values = np.linspace(0, 1000, num_steps)
dt = time_values[1] - time_values[0]

positions = np.zeros((num_steps, 3))   
positions[0] = initial_position

velocities = np.zeros((num_steps, 3))    
velocities[0] = initial_velocity


def acc(x,v):
    if np.abs(x[0])<=1:
        B=(0,b*x[0],0)
    elif x[0]>=1:
        B=(0,B0,0)
    else:
        B=(0,-B0,0)
    a=q*np.cross(v,B)/m
    return a

vh=np.zeros((num_steps, 3))
vh[0]=velocities[0]+acc(positions[0],velocities[0])*dt/2

accs = np.zeros((num_steps, 3))
accs[0] = acc(positions[0],velocities[0])

for i in range(num_steps - 1):
    positions[i+1]=positions[i]+dt*vh[i]
    vh[i+1]=vh[i]+dt*acc(positions[i+1],vh[i])
    velocities[i+1]=vh[i+1]-dt/2*acc(positions[i+1],velocities[i])





####################
##### PLOTTING #####
####################
x_vals = positions[:,0]
y_vals = positions[:,1]
z_vals = positions[:,2]

plt.figure()
plt.plot(x_vals, z_vals, color = "blue", label = "Particle trajectory")
plt.legend(loc = "upper right")
plt.title("Orbit Plots")
plt.xlim((min(x_vals), max(x_vals)))
plt.ylim((min(z_vals), max(z_vals)))
plt.xlabel("x position ")
plt.ylabel("z position ")

plt.show()

但它不能正常工作。 result plot

z 值无情地增加。这似乎是一些计算问题，但我找不到确切的问题。有人可以帮帮我吗？

【问题讨论】：

标签： python simulation

【解决方案1】：

我不习惯您使用的跳跃式算法版本。但是，我测试了您的代码，我认为罪魁祸首是 q 和 m 变量。他们采用的极小的值很可能会导致数字问题。事实上，我的 python 解释器甚至对此发出警告：

RuntimeWarning: overflow encountered in divide
  a=q*np.cross(v,B)/m

这类问题在数值模拟中非常常见，可以通过使用所谓的缩减单位轻松解决（参见例如here）。在您的示例中，设置 q = 1 和 m = 1 似乎会产生现实的results。

编辑：我想补充一点，使用缩减单位并不像将所有常量设置为 1 那样简单，因为不能单独选择物理测量单位。例如，在分子模拟（我的领域）中，习惯上将粒子直径、能量尺度和质量设置为 1。然后，所有其他数量都以这些单位表示。在您的情况下，设置 q = 1 和 m = 1 将更改 B0 和 dt 的值。

【讨论】：