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.
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.