Нахождение нулевых значений с помощью odeint - PullRequest
2 голосов
/ 11 марта 2011

Как я могу найти точку, где первая производная моего уравнения равна 0, используя scipy.integrate.ode?

Я настроил эту функцию, которая получает ответ, но я не уверен в точностии это не может быть самым эффективным способом сделать это.

В основном я использую эту функцию, чтобы найти время, когда снаряд с начальной скоростью перестает двигаться.В системах ODE есть ли лучший способ решить этот ответ?

import numpy as np
from scipy import integrate

def find_nearest(array,value):
    idx=(np.abs(array-value)).argmin()
    return array[idx], idx

def deriv(x,t):
    # This function sets up the following relations
    # dx/dt = v , dv/dt = -(Cp/m)*(4+v^2)
    return np.array([ x[1], -(0.005/0.1) * (4+ (x[1]**2)) ])

def findzero(start, stop, v0):
    time = np.linspace(start, stop, 100000)
    #xinit are initial vaules of equation
    xinit = np.array([0.0,v0]) 
    x = integrate.odeint(deriv,xinit,time)
    # find nearest velocity value nearest to 0
    value, num = find_nearest(x[:,1],0.0001)
    print 'closest value ',
    print value
    print 'reaches zero at time ',
    print time[num]
    return time[num]
# from 0 to 20 seconds with initial velocity of 100 m/s
b = findzero(0.0,20.0,100.0)

Ответы [ 2 ]

4 голосов
/ 11 марта 2011

В общем, хороший подход к решению этой проблемы - переписать ваши уравнения так, чтобы скорость была независимой переменной, а время и расстояние - зависимыми переменными.Тогда вам просто нужно интегрировать уравнения от v = v0 до v = 0.

Однако в приведенном вами примере даже не нужно прибегать к scipy.integrate.Уравнения могут быть легко решены карандашом и бумагой (разделение переменных сопровождается стандартным интегралом).Результат равен

t = (м / (2 Cp)) arctan (v0 / 2)

, где v0 - начальная скорость, и результат arctan должен бытьвзяты в радианах.

Для начальной скорости 100 м / с ответ составляет 15.5079899282 секунды.

0 голосов
/ 11 марта 2011

Я бы использовал что-то вроде scipy.optimize.fsolve(), чтобы найти корни производной.Используя это, можно работать назад, чтобы найти время, необходимое для достижения корня.

Добро пожаловать на сайт PullRequest, где вы можете задавать вопросы и получать ответы от других членов сообщества.
...