python - 在相平面中绘制 ode 解

Question

我正在尝试使用求解器集成.odeint 为我的 ode 绘制解决方案，但是当我尝试绘制它时我没有获得解决方案。我看不出我的代码的表述有什么问题。

请在下面找到： import numpy as np import matplotlib.pyplot as pl from numpy import sin, cos import numpy as np import matplotlib.pyplot as plt import scipy.integrate as integration import matplotlib.animation as animation from math import*

g = 9.81 
l = 1.6
l_big = 2.0
l_small = 1.6
m = 0.01
alpha = 0.4
k = 100
def sh2(r1,t):

    theta1,omega1 = r1 
    sh2_theta1 = omega1/(l + ((1/2)*alpha*(1+np.tanh(theta1*omega1*k))))**2

    sh2_omega1 = -g*(l + ((1/2)*alpha*(1+np.tanh(theta1*omega1*k))))*sin(theta1)

    #return sh2_theta1, sh2_omega1

    return np.array([sh2_theta1, sh2_omega1],float)

init_state = np.radians([69.0,0])
dt = 1/40
time = np.arange(0,10.0,dt)
timexo = np.arange(0,10.0,dt)

state2 = integrate.odeint(sh2,init_state,time)
#print(state2)
print(len(state2),len(timexo))

plt.plot(timexo[0:2500],state2[0:2500])
plt.show()

#phase plot attempt 
# initial values
x_0 = 0 # intial angular position
v_0 = 1.0 # initial angular momentum
t_0 = 0 # initial time

# initial y-vector from initial position and momentum
y0 = np.array([x_0,v_0]) 


# max value of time and points in time to integrate to
t_max = 10
N_spacing_in_t = 1000

# create vector of time points you want to evaluate
t = np.linspace(t_0,t_max,N_spacing_in_t)

# create vector of positions for those times
y_result = np.zeros((len(t), 2))

# loop through all demanded time points
for it, t_ in enumerate(t):

    # get result of ODE integration up to the demanded time
    y = integrate.odeint(sh2,init_state,t_)

    # write result to result vector
    y_result[it,:] = y

# get angle and angular momentum
angle = y_result[:,0]
angular_momentum = y_result[:,1]

# plot result
pl.plot(angle, angular_momentum,'-',lw=1)
pl.xlabel('angle $x$')
pl.ylabel('angular momentum $v$')

pl.gcf().savefig('pendulum_single_run.png',dpi=300)

pl.show()

此代码的输出：

图 1：随着时间的推移，ode 解的正确图
图 2：导致问题的相平面的空图。

欢迎任何提示。在不太重要的一点上，我的第一个情节给出了两条线，但我只期待蓝线。

标签： pythonnumpymatplotlibode

Answer 1

该图是空的，因为积分器在您的 for 循环中返回零。为什么首先使用 for 循环？如果您随着时间的推移进行集成，就像您在代码的第一部分中所做的那样，一切都会正常工作。 注意：您不需要导入 matplotlib.pyplot 两次。

关于 plot1 中的两条线：对象state2[0:2500]的尺寸为 2x2500。因此出现两条线。如果您只想要一行，请使用np.transpose然后state2[0]或state2[1]访问这些行。

下面的代码将解决您的问题。我添加了一个plt.figure()命令来生成两个图，而无需等待第一个凝块被关闭。

import matplotlib.pyplot as plt
import numpy as np
from numpy import sin
import scipy.integrate as integrate
from math import *

g = 9.81
l = 1.6
l_big = 2.0
l_small = 1.6
m = 0.01
alpha = 0.4
k = 100
def sh2(r1,t):

    theta1,omega1 = r1 
    sh2_theta1 = omega1/(l + ((1/2)*alpha*(1+np.tanh(theta1*omega1*k))))**2

    sh2_omega1 = -g*(l + ((1/2)*alpha*(1+np.tanh(theta1*omega1*k))))*sin(theta1)

    #return sh2_theta1, sh2_omega1

    return np.array([sh2_theta1, sh2_omega1],float)

init_state = np.radians([69.0,0])
dt = 1/40
time = np.arange(0,10.0,dt)
timexo = np.arange(0,10.0,dt)

state2 = integrate.odeint(sh2,init_state,time)

print(len(state2),len(timexo))
state2_plot = np.transpose(state2[0:2500])
plt.plot(timexo[0:2500],state2_plot[1])

#phase plot attempt 
# initial values
x_0 = 0.0 # intial angular position
v_0 = 1.0 # initial angular momentum
t_0 = 0 # initial time

# initial y-vector from initial position and momentum
y0 = np.array([x_0,v_0]) 

# max value of time and points in time to integrate to
t_max = 10
N_spacing_in_t = 1000

# create vector of time points you want to evaluate
t = np.linspace(t_0,t_max,N_spacing_in_t)

# create vector of positions for those times

y_result = integrate.odeint(sh2, init_state, t)

# get angle and angular momentum
angle = y_result[:,0]
angular_momentum = y_result[:,1]

# plot result
fig = plt.figure()
plt.plot(angle, angular_momentum,'-',lw=1)
plt.xlabel('angle $x$')
plt.ylabel('angular momentum $v$')

plt.gcf().savefig('pendulum_single_run.png',dpi=300)
plt.show()

输出：