I am trying to decelerate a rocket in space as it approaches Mars to enter its final orbit. I am trying to use solve_bvp to describe the position and the velocity of the rocket as it decelerates. The code is as follows.

When running this code, the solver fails to converge, with a message stating “A singular Jacobian encountered when solving the collocation system.” How can I fix this? Thank you for helping.

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 # Nm2/kg2
m1 = 5.97219 * 1024 # 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*(-Isp2g0+v))/((total_mass_3)-m2dot2t)
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)

I’ve tried changing various boundary values given. For example, if I change the first boundary condition to be Sa[0] instead of Sa[1], making it refer to the first value of the state variable, just as an experiment. When doing this, the code does not have an issue with Jacobian point, and if I increase the maximum number of nodes, it actually solves. Of course, this is without the boundary values I want to use. I think this means there is something wrong with my boundaries, and I am not totally sure what it could be.