I am trying to solve a problem in python using scipy's solve_ivp. The system of ODEs I am trying to solve is for coupled where I am solving for two time-dependent variables: f and T. The derivative of T is in the following format: $$ \frac{dT}{dt}=g(t)+\frac{T}{f} \times \frac{df}{dt} $$ Where $g(t)$ is some function (irrelevant in this example), where that relation can of course be interpreted through the product rule as: $\frac{d(fT)}{dt}=g(t)\times f$.

However, my initial condition for $f(t)$ is the following: $f(t=0)=0$. Which causes quite a problem for solve_ivp, as $\frac{dT}{dt}$ becomes a NaN.

I am also quite reluctant to slightly alter the initial condition and give it a symbolic value of $f(t=0)=\epsilon$, where $\epsilon \rightarrow 0$. As $\frac{df}{dt}$ is initially quite large due to the nature of the system, and if I were to divide by a very small $f$ then the gradient of $\frac{dT}{dt}$ becomes too large, making my system too stiff. And if I were to increase too much the value of $\epsilon$, then I fear the solution may not be too accurate, particularly in the earlier stages of the simulation.

The problem statement is too complex to solve analytically, so I am looking for a numerical solution, does anyone here know of a nice trick/getaround to solve this issue?


  • $\begingroup$ It's not clear how you have gotten $\frac{dT}{dt}$, but since the RHS of your system of ODE's is undefined (i.e. NaN) there is no well-defined solution. (In general to have a well defined solution RHS of your system of ODE's must be Lipschitz continuous.) So the question is if there is a way to reformulate the model or determine the differential equations in a way such that the RHS is well behaved for the range of integration. It's hard to say how this could be addressed with out more context about the model, variables and how you have come to this differential equation. $\endgroup$
    – user9794
    Commented Aug 23, 2023 at 14:51
  • 1
    $\begingroup$ To be more specific than @user9794: This sort of situation typically happens because you made a mistake in deriving your equations. If the model is wrong, choosing a different numerical technique will not help you. $\endgroup$ Commented Aug 23, 2023 at 19:36
  • $\begingroup$ You mention a coupled system, but all I see is 1 equation and 2 unknowns $f$ and $T$. Where is the second equation? As far as the division by zero goes, your problem may (under sone assumptions about what you're trying to model) be reformulated as $f d_tT = fg + T d_t f$. Then $f(0) = 0$ results in $T d_t f = 0$. As an example you can now use implicit Euler to get $$f^{k+1} \frac{T^{k+1}-T^k}{\tau} = f^{k+1} g^{k+1} + T^{k+1}\frac{f^{k+1}-f^k}{\tau}.$$ $\endgroup$
    – lightxbulb
    Commented Aug 23, 2023 at 21:24
  • $\begingroup$ You did not say how $f$ is computed. If it is a known function, then use an implicit scheme where $f$ remains at the RHS thus avoiding a singularity when $f$ is 0. Otherwise, I would also check that you have a well-defined problem, i.e. does the solution exist and is it unique: en.wikipedia.org/wiki/Picard%E2%80%93Lindel%C3%B6f_theorem. $\endgroup$
    – Laurent90
    Commented Aug 23, 2023 at 22:01


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.