I’m currently trying to solve a stiff ODE problem in Python using the LSODA method from solve_ivp. I’ve encountered various issues during implementation, and the main ones are as follows. When I compute the Euclidean norm of the solution, at some point it suddenly grows and then abruptly drops to zero. 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. 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. Euclidean norm figure 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...
