python - 求解非线性二阶 ODE

Question

我被困在这里更正我的代码。我有一个用于解决相对论运动的二阶非线性微分方程，但我找不到解决方案。这是我迄今为止尝试过的：

import matplotlib.pyplot as plt
from scipy.integrate import odeint


def f(y, t, params):
    v, a = y      # unpack current values of y
    e, E, tau, m = params  # unpack parameters
    gamma = (1-(v/c)**2)**(1/2)
    derivs = [a,      # list of dy/dt=f functions
             -((((gamma**2)*theta+1)/(tau*(gamma**5)))*a)+((e*E)/(tau*(gamma**6)*m))-(6*(gamma**2)*v*(a**2))]
    return derivs

# Parameters
tau = 6*10**-24         # Given values
E = 10**6         # electric field
e = 1.602*10**-19     # electron charge
m = 9.109*10**-31 #electron mass
c = 1
#gamma = (1-(v/c)**2)**(1/2)

# Initial values
v0 = 0.0     # initial angular displacement
a0 = 0.0     # initial angular velocity

# Bundle parameters for ODE solver
params = [e, E, tau, m]

# Bundle initial conditions for ODE solver
y0 = [v0, a0]

# Make time array for solution
tStop = 0.2
tInc = 0.05
t = np.arange(0., tStop, tInc)

t = 0
dt = 0.0001
while v0 < 0.99 :
    t+= dt
    v0+= t
    v=v0
    print(t,v)

# Call the ODE solver
v = odeint(f, y0, t, args=(params,))

while v.all() < 0.99 :
  tStop = t + dt
  print(t, tStop, v)

# Plot results
fig = plt.figure(1, figsize=(8,8))

# Plot theta as a function of time
ax1 = fig.add_subplot(311)
ax1.plot(t, v[:,0])
ax1.set_xlabel('time')
ax1.set_ylabel('v')

# Plot omega as a function of time
ax2 = fig.add_subplot(312)
ax2.plot(t, v[:,1])
ax2.set_xlabel('time')
ax2.set_ylabel('a')

# Plot a vs v
ax3 = fig.add_subplot(313)
twopi = 2.0*np.pi
ax3.plot(v[:,0]%twopi, v[:,1], '.', ms=1)
ax3.set_xlabel('v')
ax3.set_ylabel('a')
ax3.set_xlim(0., twopi)

plt.tight_layout()
plt.show()

最初的问题也是在给定这个 ODE 的情况下找到 v = 0.99c 的时间：

4. 考虑一个在大小为 E = 106 V/m 的均匀电场中最初处于静止状态的电子的运动。运动的相对论方程是：

(/m) − τ /dt (6 /dt) = /dt ()

其中=02 =6×10−24 6

a) 电子达到 v = 0.99c 需要多长时间？b) 在线性和对数图中将速度作为时间的函数绘制。

谢谢你的时间 ！

标签： pythonmatplotlibscipy