Построение одного решения на фазовой плоскости

Я пытаюсь построить решение своей оды с помощью solver integration.odeint, однако я не получаю решения, когда пытаюсь построить его. Я не вижу, где неправильная формулировка моего кода.

Пожалуйста, найдите ниже: 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: правильный график решения оды с течением времени
График 2: пустой график для фазовой плоскости, который вызывает проблемы.

Любые подсказки приветствуются. На менее важном примечании мой первый график дает две линии, хотя я ожидал только синюю линию.


python numpy matplotlib ode
person lara_ca    schedule 23.02.2019    source источник

Ответы (1)


arrow_upward
1
arrow_downward

График пуст, потому что интегратор возвращает ноль в вашем цикле 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()

Выход:

введите здесь описание изображения

person f.wue    schedule 23.02.2019