solve_ivp on a differential equation containing squared functions

Question

I was trying to solve (x’(t))^2+ax^2(t)-b=0 numerically in python. I know I can use scipy’s solve_ivp for simpler equations by defining the functions like :

def f(t, y):
   return y ** 2

which would represent the equation y’(t)=y**2. How can I write the first equation like this? I’m aware that potentially the equation can be simplified with a substituition (maybe y=x^2 or something) but I’m just interested in solving it directly (and numerically) in this case.

I tried naively solving the first equation for x’(t), meaning the definition of f in python contains math.sqrt(). This gives me a domain error (at some stage the input is negative presumably). This error occurred with many different values of an and b so I’m at a loss for how to continue.

Answer 1

ODE with non-smooth right sides can have defects under numerical solution.

In the given situation an exact solution can have segments where it is constant, x'=0, and other segments where it solves x''+ax=0, and these modes can switch at points where both mode equations are simultaneously satisfied. This gives an infinity of solutions.

So it would be best if you could select the solution behavior from the start and solve an ODE with less singularities.

Answer 2

It seemed to me that you can follow your approach of taking the square root, but making sure the values are complex to begin with to avoid the domain error.

import numpy as np
from scipy.integrate import solve_ivp

def f(t, x, a, b):
    return np.sqrt(b - a*x**2)
    # return -np.sqrt(b - a*x**2)

a = 2 + 0j
b = 1 + 0j
t0, tf = 0, 10
y0 = [1 + 0j]

res = solve_ivp(f, (t0, tf), y0, args=(a, b))  # solves successfully

Based on experiments in Mathematica, this doesn't seem to produce a correct solution, though. But then this becomes more of a math question, which you might ask on the Mathematics SE site.