Блог

Я пытаюсь решить простую механическую проблему, которая является проблемой нелинейного простого маятника.

Для этого я должен использовать odeint из scipy.

Для напоминания я должен решить нелинейную ОДУ, которая: theta'' w²*sin(theta) = 0

Вот что у меня есть на данный момент, я должен определить функцию, которая возвращает положение маятника, то есть theta(t) :

 def ODE(omega,theta0,tf):
    """solution au bout tf avec odeint"""
    
    def F(Y,t,omega):
        """second membre de l'EDO"""
        return np.array([Y[1],-omega**2*np.sin(Y[0])])
    
    Y0 = np.array([theta0,0.])
    Y = odeint(F, Y0, tf, args=(omega,))
    return Y
 

Но когда я пытаюсь выполнить проверку с помощью этого, для сравнения с линейным решением:

 # test
theta0 = 0.1 
omega = 2.58
N = 15 
T = np.linspace(0,2*np.pi/omega,N)

    
# solution with odeint
ThetaODE = ODE(omega,theta0,T)


# linear solution
ThetaL = theta0*np.cos(omega*T)

# plot
plt.figure(figsize=(15,6))
plt.subplot(1,2,1)
plt.plot(T,ThetaODE,label="ODEint")
plt.plot(T,ThetaL,label="Lineaire")
plt.xlabel('t')
plt.legend()

plt.subplot(1,2,2)
plt.plot(T,ThetaL-ThetaODE)
plt.title('Ecart entre ODEint et la sol Lineaire')
plt.xlabel('t')
 

У меня такая ошибка:

 ValueError                                Traceback (most recent call last)
<ipython-input-128-d8b2832c3f6c> in <module>
     22 
     23 plt.subplot(1,2,2)
---> 24 plt.plot(T,ThetaL-ThetaODE)
     25 plt.title('Ecart entre ODEint et la sol Lineaire')
     26 plt.xlabel('t')

ValueError: operands could not be broadcast together with shapes (15,) (15,2) 
 

Должно быть что-то, что я неправильно понимаю в отношении odeint, есть идеи, что я делаю не так?
Спасибо!

Вы должны понимать, как работает библиотека odeint. Кстати, вместо этого вы должны использовать библиотеку solve_ivp : https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html#scipy.integrate.solve_ivp

В вашей системе есть две переменные: theta и ее производная по времени. Поэтому метод интеграции ODE odeint или solve_ivp вернет массив ( Y в вашей функции ODE ) с двумя измерениями. Одно измерение — это количество переменных (2), другое — количество временных точек, которые были вычислены между ними (чтобы контролировать это, обратитесь к документации).

Итак, здесь должна быть ваша сюжетная линия plt.plot(T,ThetaL-ThetaODE[:,0]) .

Совет здесь заключается в том, чтобы прочитать документы и поискать в Интернете похожие проблемы. Я думаю, вы бы легко нашли ответ. Существует множество руководств по интеграции ODE, поэтому, пожалуйста, найдите время, чтобы ознакомиться с некоторыми из них. Например, в этом есть несколько интересных примеров: https://scipy-cookbook.readthedocs.io/items/idx_ordinary_differential_equations.html

Здесь показано, как работать с pendulum: https://scientific-python.readthedocs.io/en/latest/notebooks_rst/3_Ordinary_Differential_Equations/04_Exercices/00_Tutorials/02_ODE_tutorial_pendulum.html

Вот графическая часть вашего кода с исправлениями:

 # plot the two solution components
plt.figure(figsize=(15,6))
plt.subplot(1,2,1)
plt.plot(T,ThetaODE[:,0],label=r"$theta$")
plt.plot(T,ThetaODE[:,1],label=r"$dot{theta}$")
plt.xlabel('t')
plt.grid()
plt.legend()

# compare the time evolution of the angle with the inear solution
plt.figure(figsize=(15,6))
plt.subplot(1,2,1)
plt.plot(T,ThetaODE[:,0] ,label="ODEint", linestyle='-', marker=' ')
plt.plot(T,ThetaL, label="Lineaire", linestyle='--', marker=None)
plt.xlabel('t')
plt.legend()

plt.subplot(1,2,2)
plt.plot(T,ThetaL-ThetaODE[:,0])
plt.title('Ecart entre ODEint et la sol Lineaire')
plt.xlabel('t')