Неточный фазовый портрет с матплотлибом - PullRequest
Новая
       16

Неточный фазовый портрет с матплотлибом

0 голосов
/ 26 февраля 2019

Я пытаюсь построить фазовый портрет для уравнения, как определено в моей функции sh2 в коде ниже.Я знаю, что ожидаемый фазовый график должен быть [[ожидаемый фазовый график] [1]] [1] [1]: https://i.stack.imgur.com/y1T9Y.png.

Однако это мой результат: [[результат] [1]: https://i.stack.imgur.com/EuSOm.png.

Я использую integrate.odeint, может кто-нибудь подсказать, что я мог бы изменить в приведенном ниже списке, или если будет лучше использовать другой алгоритм, который даст мне результат, более близкий к ожидаемому,

Пожалуйста, найдите мой код: 

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 = l_big-l_small             
k = 10*(10**40)             

def sh2(r1,t):

    theta1,omega1 = r1 
    sh2_theta1 = omega1

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

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



init_state = np.radians([30.0,0])
dt = 1/10.0
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], '--m', label = r'$\theta = \frac{\pi}{6}$')
plt.xlabel('Time t (s) ')
plt.ylabel('Angular Velocity' '  '  r'$\dot{\theta}$')
plt.show()


#code for phase plot 

# 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 = 10000

# 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_velocity = y_result[:,1]

# plot result
fig = plt.figure()
plt.plot(angle, angular_velocity,'--k',lw=1)
plt.xlabel('Angle' '  ' r'$\theta$')
plt.ylabel(r'Angular Velocity' r' $\dot{\theta}$')
plt.gcf().savefig('pumping.png',dpi=300)
plt.show()

Спасибо за время тура

...