Нужна помощь в решении нелинейного ODE второго порядка в python

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

У меня есть ODE 2-го порядка:

m = 1220 k = 35600 g = 17.5 a = 450000 

и b составляет от 1000 до 10000 с шагом 500.

 x(0)= 0 x'(0)= 5 m*x''(t) + b*x'(t) + k*x(t)+a*(x(t))^3 = -m*g 

Мне нужно найти наименьшее значение b, чтобы решение никогда не было положительным. Я знаю, как должен выглядеть граф, но я просто не знаю, как использовать odeint для решения дифференциального уравнения. Это код, который у меня есть до сих пор:

 from numpy import * from matplotlib.pylab import * from scipy.integrate import odeint m = 1220.0 k = 35600.0 g = 17.5 a = 450000.0 x0= [0.0,5.0] b = 1000 tmax = 10 dt = 0.01 def fun(x, t): return (b*x[1]-k*x[0]-a*(x[0]**3)-m*g)*(1.0/m) t_rk = arange(0,tmax,dt) sol = odeint(fun, x0, t_rk) plot(t_rk,sol) show() 

Который на самом деле ничего не производит.

Есть предположения? благодаря

2 Solutions collect form web for “Нужна помощь в решении нелинейного ODE второго порядка в python”

Чтобы решить ODE второго порядка с помощью scipy.integrate.odeint , вы должны записать его как систему ОДУ первого порядка:

Я определю z = [x', x] , тогда z' = [x'', x'] , и это ваша система! Конечно, вы должны подключить свои реальные отношения:

 x'' = -(b*x'(t) + k*x(t) + a*(x(t))^3 + m*g) / m 

будет выглядеть так:

z[0]' = -1/m * (b*z[0] + k*z[1] + a*z[1]**3 + m*g)
z[1]' = z[0]

Или просто назовите его d(z) :

 def d(z, t): return np.array(( -1/m * (b*z[0] + k*z[1] + a*z[1]**3 + m*g), # this is z[0]' z[0] # this is z[1]' )) 

Теперь вы можете odeint его в odeint как таковой:

 _, x = odeint(d, x0, t).T 

( _ – пустой заполнитель для переменной x' мы сделали)

Чтобы свести к минимуму b подлежащее ограничению, что максимум x всегда отрицательный, вы можете использовать scipy.optimize.minimize . Я реализую его, фактически максимизируя максимум x , при условии, что ограничение остается отрицательным, потому что я не могу придумать, как свести к минимуму параметр без возможности инвертировать функцию.

 from scipy.optimize import minimize from scipy.integrate import odeint m = 1220 k = 35600 g = 17.5 a = 450000 z0 = np.array([-.5, 0]) def d(z, t, m, k, g, a, b): return np.array([-1/m * (b*z[0] + k*z[1] + a*z[1]**3 + m*g), z[0]]) def func(b, z0, *args): _, x = odeint(d, z0, t, args=args+(b,)).T return -x.max() # minimize negative max cons = [{'type': 'ineq', 'fun': lambda b: b - 1000, 'jac': lambda b: 1}, # b > 1000 {'type': 'ineq', 'fun': lambda b: 10000 - b, 'jac': lambda b: -1}, # b < 10000 {'type': 'ineq', 'fun': lambda b: func(b, z0, m, k, g, a)}] # func(b) > 0 means x < 0 b0 = 10000 b_min = minimize(func, b0, args=(z0, m, k, g, a), constraints=cons) 

Я не думаю, что вы можете решить свою проблему, как указано: ваши начальные условия с x = 0 и x' > 0 означают, что решение будет положительным для некоторых значений, очень близких к начальной точке. Таким образом, нет b для которого решение никогда не будет положительным …

Оставив это в стороне, чтобы решить дифференциальное уравнение второго порядка, вам сначала нужно переписать его как систему двух дифференциальных уравнений первого порядка. Определяя y = x' мы можем переписать ваше единственное уравнение следующим образом:

 x' = y y' = -b/m*y - k/m*x - a/m*x**3 - g x[0] = 0, y[0] = 5 

Поэтому ваша функция должна выглядеть примерно так:

 def fun(z, t, m, k, g, a, b): x, y = z return np.array([y, -(b*y + (k + a*x*x)*x) / m - g]) 

И вы можете решать и строить свои уравнения:

 m, k, g, a = 1220, 35600, 17.5, 450000 tmax, dt = 10, 0.01 t = np.linspace(0, tmax, num=np.round(tmax/dt)+1) for b in xrange(1000, 10500, 500): print 'Solving for b = {}'.format(b) sol = odeint(fun, [0, 5], t, args=(m, k, g, a, b))[..., 0] plt.plot(t, sol, label='b = {}'.format(b)) plt.legend() 

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

  • определить численную оценку производной симплексной функции
  • Заполнение пробелов в массиве numpy
  • Эффективно умножьте плотную матрицу на разреженный вектор
  • Как читать часть двоичного файла с numpy?
  • «ImportError: без модуля scipy» после установки пакета scipy
  • Многопоточные вызовы целевой функции scipy.optimize.leastsq
  • Переключение между различными программами, использующими scipy ode
  • Должен ли я использовать scipy.pi, numpy.pi или math.pi?
  • Python - лучший язык программирования в мире.