在相平面中绘制 ode 解

Question

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

请在下面找到：将 numpy 导入为 np 将 matplotlib.pyplot 导入为 pl 从 numpy 导入 sin, cos 将 numpy 导入为 np 将 matplotlib.pyplot 导入为 plt 导入 scipy.integrate 作为集成导入 matplotlib.animation 作为动画来自数学导入*

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：相平面的空图是什么导致了问题。

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

Answer 1

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

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

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

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()

输出：

在相平面中绘制 ode 解

Plotting ode solution in phase plane

python

numpy

matplotlib

ode