Как избежать ошибок с сингулярными якобианами вsolve_bvp? ⇐ Python

[Anonymous]

Я пытаюсь замедлить ракету в космосе, когда она приближается к Марсу для выхода на последнюю орбиту. Я пытаюсь использоватьsolve_bvp для описания положения и скорости ракеты при ее замедлении. Код следующий.
При запуске этого кода решателю не удается сходиться, и появляется сообщение: «При решении системы коллокации обнаружен сингулярный якобиан». Как я могу это исправить? Спасибо за помощь.

Код: Выделить всё

import numpy as np
from scipy.integrate import solve_ivp, solve_bvp
import matplotlib.pyplot as plt

# Define necessary constants
Cd = .2
A = 11.4  # m^2
G = 6.673 * 10**-11 # Nm**2/kg**2
m1 = 5.97219 * 10**24  # kg
re = 6.371 * 10**6  # m
Isp = 300  # sec
Isp2 = 450 # sec
g0 = 9.81  # m/s^2
p0 = 101325  # Pa
M = .0289652  # kg/mol
R = 8.31446  # J/(mol*K)
T0 = 288.15  # K
L = .0065  # K/m
dry_mass = 21000 # kg
slv_mass = 2300 # kg
payload_mass = 2180 # kg
stage1_fuel_mass = 200000 # kg
stage2_fuel_mass = 50000 # kg
deceleration_fuel_mass = stage2_fuel_mass
total_mass = dry_mass + stage1_fuel_mass + stage2_fuel_mass + payload_mass + slv_mass + deceleration_fuel_mass
total_mass_2 = stage2_fuel_mass + payload_mass + slv_mass + deceleration_fuel_mass
m2dot = 2100 # kg/s
m2dot2 = 400 # kg/s

total_mass_3 = payload_mass + slv_mass + deceleration_fuel_mass

required_velocity = 3000
#ta,tb = sol2.t[-1],sol2.t[-1]+900
ta,tb = 0,900
#va,vb = sol2.y[1][-1], required_velocity
va,vb = 40000, required_velocity

def fun3(t,S):
h,v = S
dhdt = v
dvdt = (m2dot2*(-Isp2*g0+v))/((total_mass_3)-m2dot2*t)
return np.vstack([dhdt,dvdt])

def bc(Sa,Sb):
bc1 = Sa[1] - va
bc2 = Sb[1] - vb
return np.array([bc1,bc2])

t_bvp = np.linspace(ta,tb,1000)
S = np.zeros((2,t_bvp.size))
sol3 = solve_bvp(fun3, bc, t_bvp, S, max_nodes = 10000)

print(sol3)

Я пробовал изменять различные заданные граничные значения. Например, если я изменю первое граничное условие на Sa[0] вместо Sa[1], сделав его ссылкой на первое значение переменной состояния, просто в качестве эксперимента. При этом в коде не возникает проблем с точкой Якобиана, и если я увеличиваю максимальное количество узлов, это фактически решается. Конечно, это без граничных значений, которые я хочу использовать. Я думаю, это означает, что с моими границами что-то не так, и я не совсем уверен, что это может быть.

Подробнее здесь: https://stackoverflow.com/questions/784 ... -solve-bvp