Как решить систему из 9 уравнений нелинейной DE с питоном?

Question

Я отчаянно пытаюсь решить (и отобразить график) систему, состоящую из девяти нелинейных дифференциальных уравнений, которые моделируют путь бумеранга. Система следующая:

Все буквы слева являются переменными, остальные являются либо константами, либо известными функциями в зависимости от v_G и w_z

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

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

РЕДАКТИРОВАТЬ: Извините, если я не был достаточно ясно. Поскольку он моделирует путь бумеранга, моя цель не состоит в том, чтобы аналитически решить эту систему (то есть меня не волнует математическое выражение каждой функции), а скорее получить значения каждой функции для определенного временного диапазона (скажем, от t1 = 0 с до t2 = 15 с с интервалом 0,01 с между каждым значением) для отображения графика каждой функции и графика центра масс бумеранга (X, Y, Z - его координаты).

Вот код, который я пробовал:

import scipy.integrate as spi
import numpy as np

#Constants

I3 = 10**-3
lamb = 1
L = 5*10**-1
mu = I3
m = 0.1
Cz = 0.5
rho = 1.2
S = 0.03*0.4
Kz = 1/2*rho*S*Cz
g = 9.81

#Initial conditions

omega0 = 20*np.pi
V0 = 25
Psi0 = 0
theta0 = np.pi/2
phi0 = 0
psi0 = -np.pi/9
X0 = 0
Y0 = 0
Z0 = 1.8

INPUT = (omega0, V0, Psi0, theta0, phi0, psi0, X0, Y0, Z0) #initial conditions 


def diff_eqs(t, INP):  
    '''The main set of equations'''

    Y=np.zeros((9))

    Y[0] = (1/I3) * (Kz*L*(INP[1]**2+(L*INP[0])**2))
    Y[1] = -(lamb/m)*INP[1]
    Y[2] = -(1/(m * INP[1])) * ( Kz*L*(INP[1]**2+(L*INP[0])**2) + m*g) + (mu/I3)/INP[0]
    Y[3] = (1/(I3*INP[0]))*(-mu*INP[0]*np.sin(INP[6]))
    Y[4] = (1/(I3*INP[0]*np.sin(INP[3]))) * (mu*INP[0]*np.cos(INP[5]))
    Y[5] = -np.cos(INP[3])*Y[4]
    Y[6] = INP[1]*(-np.cos(INP[5])*np.cos(INP[4]) + np.sin(INP[5])*np.sin(INP[4])*np.cos(INP[3]))
    Y[7] = INP[1]*(-np.cos(INP[5])*np.sin(INP[4]) - np.sin(INP[5])*np.cos(INP[4])*np.cos(INP[3]))
    Y[8] = INP[1]*(-np.sin(INP[5])*np.sin(INP[3]))


    return Y   # For odeint



t_start = 0.0
t_end = 20
t_step = 0.01
t_range = np.arange(t_start, t_end, t_step)

RES = spi.odeint(diff_eqs, INPUT, t_range)

Тем не менее, я продолжаю получать ту же проблему, как показано здесь и особенно сообщение об ошибке:

Excess work done on this call (perhaps wrong Dfun type)

Я не совсем уверен, что это значит, но похоже, что у решателя проблемы с решением системы. В любом случае, когда я пытаюсь отобразить трехмерную траекторию благодаря координатам XYZ, я просто получаю 3 или 4 точки, где должно быть что-то вроде 2000.

Итак, мои вопросы: - Что-то не так в моем коде? - Если нет, есть ли другой, возможно, более сложный инструмент для решения этой системы? - Если нет, то возможно ли вообще получить то, что я хочу, от этой системы ОДУ?

Заранее спасибо

Answer 1

Есть несколько вопросов:

• если я скопирую код , он не запустится
• указанный вами обходной путь не работает с odeint, данный решение использует оду
• Справочник по scipy для odeint гласит: «Для нового кода используйте scipy.integrate.solve_ivp для решения дифференциального уравнения. "
• вызов RES = spi.odeint (diff_eqs, INPUT, t_range) должен быть соответствует функции head def diff_eqs (t, INP) . Следите за заказ: RES = spi.odeint (diff_eqs, t_range, INPUT)

Есть некоторые проблемы с математическими формулами:

• взгляните на 3-ю формулу на вашей картинке. У него нет термина тенденции, оно начинается с нуля - что это значит?
• Трудно проверить, правильно ли вы перевели формулу в код, поскольку код строго не следует формулам.

Ниже я попробовал решение с scipy solve_ivp . В случае A я могу запустить маятник, но в случае B не найдено никакого значимого решения для бумеранга. Так что проверяйте математику, я думаю, в математических выражениях есть какая-то ошибка.

Для графики используйте pandas, чтобы построить все переменные вместе (см. Код ниже).

import scipy.integrate as spi
import numpy as np
import pandas as pd

def diff_eqs_boomerang(t,Y):   
    INP = Y
    dY = np.zeros((9))
    dY[0] = (1/I3) * (Kz*L*(INP[1]**2+(L*INP[0])**2))
    dY[1] = -(lamb/m)*INP[1]
    dY[2] = -(1/(m * INP[1])) * ( Kz*L*(INP[1]**2+(L*INP[0])**2) + m*g) + (mu/I3)/INP[0]
    dY[3] = (1/(I3*INP[0]))*(-mu*INP[0]*np.sin(INP[6]))
    dY[4] = (1/(I3*INP[0]*np.sin(INP[3]))) * (mu*INP[0]*np.cos(INP[5]))
    dY[5] = -np.cos(INP[3])*INP[4]
    dY[6] = INP[1]*(-np.cos(INP[5])*np.cos(INP[4]) + np.sin(INP[5])*np.sin(INP[4])*np.cos(INP[3]))
    dY[7] = INP[1]*(-np.cos(INP[5])*np.sin(INP[4]) - np.sin(INP[5])*np.cos(INP[4])*np.cos(INP[3]))
    dY[8] = INP[1]*(-np.sin(INP[5])*np.sin(INP[3]))
    return dY   

def diff_eqs_pendulum(t,Y): 
    dY = np.zeros((3))
    dY[0] =  Y[1]
    dY[1] = -Y[0]
    dY[2] =  Y[0]*Y[1]
    return dY

t_start, t_end = 0.0, 12.0

case = 'A'

if case == 'A':         # pendulum
    Y = np.array([0.1, 1.0, 0.0]); 
    Yres = spi.solve_ivp(diff_eqs_pendulum, [t_start, t_end], Y, method='RK45', max_step=0.01)

if case == 'B':          # boomerang
    Y = np.array([omega0, V0, Psi0, theta0, phi0, psi0, X0, Y0, Z0])
    print('Y initial:'); print(Y); print()
    Yres = spi.solve_ivp(diff_eqs_boomerang, [t_start, t_end], Y, method='RK45', max_step=0.01)

#---- graphics ---------------------
yy = pd.DataFrame(Yres.y).T
tt = np.linspace(t_start,t_end,yy.shape[0])
with plt.style.context('fivethirtyeight'): 
    plt.figure(1, figsize=(20,5))
    plt.plot(tt,yy,lw=8, alpha=0.5);
    plt.grid(axis='y')
    for j in range(3):
        plt.fill_between(tt,yy[j],0, alpha=0.2, label='y['+str(j)+']')
    plt.legend(prop={'size':20})