solve_ivp в SciPy с использованием метода Randau

Question

Я пытаюсь решить ODE первого порядка (один ODE, а не система) в Python, используя solve_ivp в SciPy с методом Randau.ODE очень жесткий, а именно, я использую метод solve_ivp обновленного SciPY.Для сложных задач SciPY рекомендует метод интеграции Randau.Прикрепленный фрагмент кода:

for gamma in gamma_array:
    for q_nondim in q_nondim_array:
        zs = []
        us = []
        r = solve_ivp(model, (0,1),[0], method = 'Radau',jac = jacobian,t_eval= [0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8,0.9, 1.0])
        if r.status <0 :
            print(r.message)
            print(r.success)
        us = np.array(r.y)
        zs = np.array(r.t)
        pressure = pressure_generator(us)
        pressure_array[i,j]= pressure[-1]
        j = j+1
    i = i+1
    j =0

Здесь модель - это вызываемая функция, которая вычисляет du / dz:

def model(t,y):  
    A1 = (1+t_by_a)*(1+t_by_a)  
    A2 = A1 + y*y +2*y  
    A3 = 2*gamma*(np.power((3+1/n)*q_nondim,n))  
    Dr1 = 2*A1*(1+y)/(A2*A2) 
    Dr2 = 2*(1+y)/(A2) 
    Dr3 = 2/(1+y) 
    Dr4 = 2/((1+y)*(1+y)*(1+y)) 
    A = np.power((1+y),(1+3*n)) 
    dudz = -A3/(A*(Dr1+Dr2-Dr3-Dr4)) 
    return dudz

jacobian - это функция, которая возвращает якобиан ODE в видемассив формы (1,1)

def jacobian(t,y):
    C = 2*gamma*pow((3+1/n)*q_nondim,n) 
    a= t_by_a
    b=((C*(a**2 + 2*a + (y + 1)**2))*((y + 1)**(1 - 3*n))* ((a**4)* (3*n*(y**2 + 2*y + 2) - 2) + 4*(a**3)*(3*n*(y**2 + 2*y + 2) - 2) + (a**2)* (3*n*(y**4 + 4*y**3 + 13*y**2 + 18*y + 12) - 2*(2*y**4 + 8*y**3 + 15*y**2 + 14*y + 9)) + 2*a*((y + 1)**2)*(3*n*(y**2 + 2*y + 4) - 2*(2*y**2 + 4*y + 5)) + 2*(3*n - 4)*(y + 1)**4))/(2*a*(a + 2)*((a**2)*(y**2 + 2*y + 2) + 2*a*(y**2 + 2*y + 2) + 2*(y + 1)**2)**2)
    c = np.array([b]) #convert a list into a 1D array
    d = np.reshape(c,(1,1)) 
    return c

Генератор давления - это функция, которая принимает решение ОДУ и преобразует его в требуемую переменную (в данном случае давление).

def pressure_generator(u):

A1 = (1+t_by_a)*(1+t_by_a)
A2 = A1 + u*u +2*u
A3 = A1/A2
F1 = A3  + np.log(A3)
A4 = 1/((1+u)*(1+u))
F2 = A4 + np.log(A4)
p = (1/gamma)*(F1 - F2)
return p

Однако интеграция здесь, похоже, терпит неудачу.Значение r.success равно False, а значение r.message:

Требуемый размер шага меньше интервала между числами.

Я попытался снова, изменив массив t_eval на:

t_eval= [0.0,0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95, 1.0]

, но интеграция не удалась даже тогда с тем же сообщением.Более того, я получил другое сообщение об ошибке:

Value error: setting array element with a sequence.

, относящееся к:

pressure_array[i,j]= pressure[-1]

строка кода.

Пожалуйста, помогите мне всделать эту интеграцию и решить ODE.Поскольку solve_ivp является новым дополнением к SciPy, литература, относящаяся к нему, довольно ограничена, и официальная документация мало помогает.

Между тем, ODE приводится ниже: ODE_For_SciPY