在阶段计划中绘制ode解决方案

Question

我试着用解算器绘制出我的颂歌的解决方案集成.odeint然而，当我不想得到一个解决方案的时候。 我看不出我的代码在哪里是错误的。在

请查找以下内容： 将numpy作为np导入 进口matplotlib.pyplot作为pl 从纽比进口罪恶，科斯 将numpy作为np导入 进口matplotlib.pyplot作为plt 进口整合作为整体 进口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循环中返回零。为什么首先使用for循环？ 如果随着时间的推移进行集成，就像在代码的第一部分中所做的那样，一切都会正常工作。 注意：您不需要导入matplotlib.pyplot两次。在

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

下面的代码将解决您的问题。我添加了一个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()

输出：