[Python] Scipy solve_ivp issue: divergence time step problem

I’m currently trying to solve a stiff ODE problem in Python using the LSODA method from solve_ivp.

The governing equation is u' = Lu

matrix information of u

matrix information of L

here, u is defined as complex 2M-vector, and L is complex 2M x 2M matrix and L varies by time. The definition of matrices is complex, but it is well calculated externally, and I only solve the problem with given datas. The main code implementation is as follows.

def compute_matrices(psi):

Q = M - N * q(psi)

A = ((2 * np.pi) ** 2) * (
N @ (G_22(psi) @ N + G_23(psi) @ M)
+ M @ (G_23(psi) @ N + G_33(psi) @ M)
)

B = -2 * np.pi * 1j * chi_r(psi) * (
N @ (G_22(psi) + q(psi) * G_23(psi))
+ M @ (G_23(psi) + q(psi) * G_33(psi))
)

C = (
(-2 * np.pi * 1j) * (
chi_rr(psi) * (N @ G_22(psi) + M @ G_23(psi))
+ qchi_r_derivative(psi) * (N @ G_23(psi) + M @ G_33(psi))
)
- ((2 * np.pi) ** 2) * chi_r(psi) * (
(N @ G_12(psi) + M @ G_31(psi)) @ Q
)
- np.pi * 1j * (
Q @ J_theta(psi) + J_theta(psi) @ Q
+ (P_r(psi) / chi_r(psi)) * (N @ J(psi) + J(psi) @ N)
)
)

D = (chi_r(psi) ** 2) * (
G_22(psi) + q(psi) * G_23(psi)
+ q(psi) * (G_23(psi) + q(psi) * G_33(psi))
)

E = (
chi_r(psi) * (
chi_rr(psi) * (G_22(psi) + q(psi) * G_23(psi))
+ qchi_r_derivative(psi) * (G_23(psi) + q(psi) * G_33(psi))
)
- 2 * np.pi * 1j * (chi_r(psi) ** 2) * (
(G_12(psi) + q(psi) * G_31(psi)) @ Q
)
+ P_r(psi) * J(psi)
)

H = (
chi_rr(psi) * (
chi_rr(psi) * G_22(psi) + qchi_r_derivative(psi) * G_23(psi)
)
+ qchi_r_derivative(psi) * (chi_rr(psi) * G_23(psi) + qchi_r_derivative(psi) * G_33(psi))

+ 2 * np.pi * 1j * chi_r(psi) * (
chi_rr(psi) * (Q @ G_12(psi) - G_12(psi) @ Q)
+ qchi_r_derivative(psi) * (Q @ G_31(psi) - G_31(psi) @ Q)
)
+ ((2 * np.pi * chi_r(psi)) ** 2) * (Q @ G_11(psi) @ Q)
+ (
P_r(psi) * ((J(psi) * chi_rr(psi) / chi_r(psi)) + J_r(psi))
+ q_r(psi) * chi_r(psi) * J_theta(psi)
)
)

A_inv = np.linalg.inv(A)
F = D - B.conj().T @ A_inv @ B
K = E - B.conj().T @ A_inv @ C
G = H - C.conj().T @ A_inv @ C
F_inv = np.linalg.inv(F)

L_upper = np.hstack([-F_inv @ K, F_inv])
L_lower = np.hstack([G - K.conj().T @ F_inv @ K, K.conj().T @ F_inv])
L = np.vstack([L_upper, L_lower])

return L, F, K


The details of the matrices is too long to show here, but I've checked this is not the problem of interpolation or inversion of matrices.

def ode_system(psi, U_flat):

U = np.ascontiguousarray(U_flat).reshape(2 * size, -1).view(dtype=np.complex128)
L, _, _ = compute_matrices(psi)
U_dot = L @ U # (2M, N)

return np.ascontiguousarray(U_dot).view(dtype=np.float64).ravel()

psi_range = np.linspace(1e-1, 0.2, 20)
_, F0, K0 = compute_matrices(psi_range[0])

u0_complex = np.vstack([ # (2 * size, size)
np.zeros((size, size), dtype=complex),
F0 @ np.eye(size, dtype=complex)
])

u0_real = u0_complex.view(dtype=np.float64).ravel()
ucrit = 1e4
current_t = psi_range[0]
final_t = psi_range[-1]
sing_flag = False
test_flag = False
solutions_t = []
solutions_y = []
current_u_real = u0_real.copy()

sol = solve_ivp(
ode_system,
[current_t, final_t],
current_u_real,
method='LSODA',
atol=1e-3,
rtol=1e-5,
max_step=1e-3,
min_step=1e-5,
dense_output=True
)

solutions_t.append(sol.t)
solutions_y.append(sol.y)


I’ve encountered various issues during implementation, and the main ones are as follows.

1. When I compute the Euclidean norm of the solution, at some point it suddenly grows and then abruptly drops to zero. I've found the complex part of the solution matrix becomes suddenly zero. And this phenomenon appears at different position by integration rangel. It’s difficult to pinpoint the cause of this phenomenon. Even if I switch the solver to RK45, Radau, or BDF, the same blow-up and drop-to-zero issue still occurs. This even makes the different solution of IVP problem.

t = [0.1, 0.2] integration result

t = [0.1, 0.15] integration result

1. LSODA is supposed to advance quickly through non-stiff regions and reduce its step size in stiff regions, but compared to RK45, it runs nearly 100 times slower. Since RK45 seems to yield a proper solution, I’m wondering whether LSODA is really moving fast enough in the non-stiff portions.

I'm trying to solve a complex problem, which is n-coupled 1st order ODE. But, I think this is a numerical issue.

Continue reading...