Моделирование маятника в Python с использованием odeint () не совсем работает как маятник

Я построил моделирование маятника, используя дифференцирование Рунге-Кутты четвертого порядка, где все делается шаг за шагом:

from scipy import *
from matplotlib.pyplot import *
##A pendulum simulation using fourth order 
##Runge-Kutta differentiation

ts=.05 #time step size
td=20 #trial duration
te=int(td/ts) #no of timesteps

mu=0.1 #friction factor
m=1 #mass
g=9.81 #grav. acceleration
l=1 #length

th=[((rand()*2)-1)*pi] #initial angle
om=[0] #initial angular velocity
u=0 #torque

for j in range(te):
    #Euler approximation
    th.append(th[j] + ts*om[j])
    f1 = (-mu*om[j] + m*g*l*sin(th[j]) + u)/(m*(l^2))
    om.append(om[j] + ts*f1)

    #approximation 1 at mid-interval
    th2 = th[j+1] + (ts/2)*om[j+1]
    f2 = (-mu*om[j+1] + m*g*l*sin(th[j+1]) + u)/(m*(l^2))
    om2 = om[j+1] + (ts/2)*f2

    #approximation 2 at mid-interval
    th3 = th2 + (ts/2)*om2
    f3 = (-mu*om2 + m*g*l*sin(th2) + u)/( m*(l^2))
    om3 = om2 + (ts/2)*f3

    #approximation at next time step
    th4 = th3 + (ts)*om3
    f4 = (-mu*om3 + m*g*l*sin(th3) + u)/( m*(l^2))
    om4 = om3 + (ts)*f4

    dth=(om[j] + 2*om[j+1] + 2*om2 + om3)/6
    dom=(f1 + 2*f2 + 2*f3 + f4)/6
    th[j+1] = th[j] + ts*dth
    om[j+1] = om[j] + ts*dom

subplot(211),plot(th),xlabel('Angle'),ylabel('')
subplot(212),plot(om,'r'),xlabel('Angular velocity'),ylabel('')
show()

Это прекрасно работает, но мне не кажется элегантным, поэтому я решил переписать его с помощью функции odeint ():

from scipy import integrate
from scipy import *
from matplotlib.pyplot import *
##A pendulum simulation

th=((rand()*2)-1)*pi #initial angle
om=0 #initial angular velocity
u=0 #torque

y0 = [om, th] #initial values
t = linspace(0, 20, 400) #

def f(y, t):
    mu=0.1 #friction factor
    m=1 #mass
    g=9.81 #grav. acceleration
    l=1 #length

    return (y[0],
        (-mu*y[0] + m*g*l*sin(y[1]) + u)/(m*(l^2)))


r = integrate.odeint(f, y0, t)

subplot(211),plot(t,r[:,1]),xlabel('Angle'),ylabel('')
subplot(212),plot(t,r[:,0],'r'),xlabel('Angular velocity'),ylabel('')
show()

И эта новая версия совсем не похожа на маятник. Кто-нибудь может сказать мне, что я делаю не так?


python differential-equations differentiation
person Fanario Jaihi    schedule 17.06.2012    source источник

Ответы (1)


arrow_upward
1
arrow_downward

У вас есть ошибки в производной функции движения маятника.

В python ^ - это оператор исключающего ИЛИ, а не оператор мощности, которым является **. Так что вы будете использовать l ** 2. Кроме того, AFAICT, первый элемент вектора состояния должен быть y [1]. Из аналогичных работ в прошлом это сработало (предполагалось, что m = 1 и l = 1).

def simple_pendulum_deriv(x, t, g = 9.81, mu = 0.5): 
    nx = np.zeros(2)
    nx[0] = x[1]
    nx[1] = -(g * np.sin(x[0])) - mu*x[1]
    return nx
person Andrew Walker    schedule 17.06.2012
comment
Спасибо за ответ. Я уже нашел ошибку с показателем степени. Это было моей проблемой из-за того, что я перенес код с другого языка и не помнил, что Python делает другие вещи с ^. Что касается вектора состояния, я не уверен, что понимаю, о чем вы. Вы имеете в виду, что я должен переключать угол и угловую скорость? edit: Я переключил их, и это сработало, хотя я не знаю, почему, вектор должен быть спереди назад или что-то в этом роде? Еще раз спасибо за ответ, вы легенда. - person Fanario Jaihi; 17.06.2012