Solve 1st order ODE in using `scipy`

Question

I've been trying to solve the following equation $$ y(t)=-A\cdot\frac{\mathrm{d} y}{\mathrm{d} t}+B\cdot\left(\frac{\mathrm{d} y}{\mathrm{d} t}\right)^{2}+C \\ y(t=0)=y_{0}\\ $$ where $A$, $B$, and $C$ are positive numbers. The derivative at $t=0$ is also known but not really needed. I tried to solve it using scipy.integrate.solve_ivp but all the tutorials I watched don't have the derivate in the function and I can't figure out how to decouple it.

Answer 1

One preliminary comment is that supplying at the initial point both $y$ and $y'$ is redundant. If $y'(t_0)$ is given then obviously the ODE fully defines $y(t_0)$. On the other hand, if $y(t_0)$ is given then in general there are two different values of $y'(t_0)$ that satisfy the equation since the ODE can be interpreted as a quadratic equation for $y'$.

Extending this argument further, we solve for $y'$ in terms of $y$ and parameters $A$,$B$,$C$:

$$ y' = - \frac{A}{2B} \pm \frac{1}{2B} \sqrt{A^2 - 4 B (C-y)} $$

There are two branches here, and each of them defines an ODE in the canonical form $y_t=F(t,y)$ which can be solved with any standard method. Below is a short Python code solving the problem for illustrative values of the parameters: A=1, B=-1, C=1, y(0)=-100. There are two branches of the solution, and in this situation, one would usually need to invoke some physics arguments to establish if both branches are meaningful.

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

#-parameters
A=1; B=-1; C=1

#-two branches
Fp = lambda t, y: -0.5*A/B + (0.5/B)*np.sqrt(A*A - 4*B*(C-y))
Fm = lambda t, y: -0.5*A/B - (0.5/B)*np.sqrt(A*A - 4*B*(C-y))

t_eval = np.arange(0, 10., 0.1)
solp = solve_ivp(Fp, [0, 10.], [-100.], t_eval=t_eval)
solm = solve_ivp(Fm, [0, 10.], [-100.], t_eval=t_eval)

plt.plot(solp.t, solp.y[0])
plt.plot(solm.t, solm.y[0])
plt.xlabel('t')
plt.ylabel('y(t)')
plt.show()