洛伦兹公式的数值实验，用python的jupyterbook写

python
洛伦兹公式中。

（1）只有均匀磁场情况下，粒子运动轨迹。

（2）均匀电场和磁场（二者相互垂直），粒子运动轨迹。

电场，电荷分布，初始条件等等自己定。

比如每隔0.01秒，算一次加速度，一次速度，一次坐标，计算10秒的吧
（**能明显看出粒子运动轨迹的形状**）

看不懂问题可以问问，希望帮忙写写代码，我还不会自己写代码。

import  numpy as np
from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d import Axes3D


p0=np.array([0,0,0]) # 初始位置 m 
v0=np.array([1,0,1])  #初速度    m/s
m=1   #  粒子质量   kg
e=np.array([0,1,0]) # 电场强度 N/c
b=np.array([0,0,3])  # 磁感应强度 T
q=1                 # 电荷量  C
ts=0.01         # 步长 s
t=10       #      总的模拟时间  s


#  匀强电场和磁场  所以加速度和位置无关
def step(p0,v0,e,b,q):
    global ts

    f=e*q+q*(np.cross(v0,b))
    a=f/m
    p=p0+v0*ts+a*ts**2/2
    v=v0+a*ts

    return  p,v


res1=[]
res1.append(p0.tolist())
res2=[]
res2.append(p0.tolist())

p,v=p0,v0
for i in  range(int(t/ts)):
    p,v=step(p,v,0,b,q)
    res1.append(p.tolist())

res1=np.array(res1) 

fig=plt.figure()
ax1 = fig.add_subplot(121,projection='3d')
ax1.scatter3D(res1[:,0],res1[:,1],res1[:,2])
ax1.set_xlabel('x')
ax1.set_ylabel('y')
ax1.set_zlabel('z')
ax1.set_title('1')

p,v=p0,v0
for i in  range(int(t/ts)):
    p,v=step(p,v,e,b,q)
    res2.append(p.tolist())
res2=np.array(res2)
ax2 = fig.add_subplot(122,projection='3d')
ax2.scatter3D(res2[:,0],res2[:,1],res2[:,2])
ax2.set_xlabel('x')
ax2.set_ylabel('y')
ax2.set_zlabel('z')
ax2.set_title('2')
plt.show()

因为近视将每步长之间看成匀变速运动，所以模拟结果和理论结果又一定的偏差，但这是使用步长这种方法不可避免的，希望代码能帮助到你。

import numpy as np
from scipy import integrate
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

def lorenz(p,t,s,r,b):
    x,y,z = p.tolist()          #无质量点的当前位置(x,y,z)
    print("x,y,z,t:",x,y,z,t)   #帮助理解odeint的执行过程
    return s*(y-x),x*(r-z)-y,x*y-b*z #返回dx/dt,dy/dt,dz/dt

t = np.arange(0,30,0.01)
track1 = integrate.odeint(lorenz,(0.0,1.00,0.0),t,args=(10.0,28.0,2.6))
track2 = integrate.odeint(lorenz,(0.0,1.01,0.0),t,args=(10.0,28.0,2.6))
print("type(track1):",type(track1),"track1.shape:",track1.shape)

fig = plt.figure(figsize=(12,6))
ax = fig.gca(projection='3d')   #获取当前子图，指定三维模式
ax.plot(track1[:,0],track1[:,1],track1[:,2],lw=1.0,color='r')   #画轨迹1         
ax.plot(track2[:,0],track2[:,1],track2[:,2],lw=1.0,color='g')   #画轨迹2         
...
plt.show()

https://blog.csdn.net/SeaBiscuitUncle/article/details/103944319