Approximating the solution of a non-linear ODE using Python

Question

This is my first time asking a question here, so please tell me if I have made a mistake or if anything is unclear.

I am working on my high school research project on the motion of a ball falling through flowing water and derived eq. 1.
$b$ is the height of the ball in the situation. This equation was derived by adding the effects of gravity, buoyancy, and drag forces.
I want to use the equation to find or approximate the value of $b$ in terms of $t$, given the values of other constants.

Eq. 1: $$\frac{db}{dt}=\int_{0}^{t}{-g}\,dt+\int_{t_1}^{t}\left[\frac{3g\rho_w}{4r^3\rho_i}\left[\left(h_2-b+r\right)^2r-\frac{1}{3}\left(h_2-b+r\right)\right]+\frac{3r\rho_w}{16\rho_i}\left(\frac{db}{dt}-\frac{1}{2}k^2t-\frac{1}{2}k\sqrt{h_0}\right)^2\right]dt$$

$b, h_2$ are function of $t$, time. Everything else $(g, \rho_w, \rho_i, r, k, h_0, t_1)$ is constant.

However, $h_2$ in the formula is a function of $t$ and is given the following formula which includes $b$ :

Eq.2 $$h_2=\frac{A_c+\frac{1}{3}+2\pi r b-2\pi r^2\pm\sqrt{\left(A_c+\frac{1}{3}\right)^2-4\pi A_cr\left(\frac{1}{4}k^2t^2+\frac{1}{2}k\sqrt{h_0}t+h_0-b+r\right)}}{2\pi r}$$

Everything in this formula except $t$ and $b$, so $\pi, A_c, r, k, h_0$, are constants.

The issue I have here is that the eq. 1 has $b^2, b, \frac{db}{dt}, \left(\frac{db}{dt}\right)^2$ inside an integral and is already quite chaotic, plus the formula for $h_2$ in it is described in $b$. To make it even worse, that formula is in the form of a solution of a quadratic equation and is not in a shape which can be easily integrated.

I would like to find the approximated formula of $b\,(t)$ using Python codes.
Any help or suggestions are appreciated! Thank you!

Answer 1

The canonical form of Newton's law for a particle in the classical mechanics is

$ \ddot{x}= f(t,x,\dot{x}) $

That is, the second time derivative of the coordinate x is a function of time, space, and the first time derivative of x. One exception is the "radiation friction" discussed in the literature, but it is not relevant here, so let's leave that special case aside. Note that the RHS function $f$ is just an algebraic or transcendental expression involving local (in time and space) values; in other words, the law of motion is a differential (not integro-differential) equation.

If we take the time derivative, our ODE for $b$ in this problem is turned in the canonical form, i.e.,

$ \ddot{b}= f(t,b,\dot{b}), $

To deal with it, let's write it as a system of 1st order ODEs,

$ \dot{v}= f(t,b,v)\\ \dot{b}= v $

At this point, we can just plug our system of two 1st order ODEs in a standard integrator, and Python has a bunch of those available, e.g., the LSODA integrator from SciPy.

Answer 2

please find my suggestions below.

For Computing b

Ignoring the $\frac{db}{dt}$ on the RHS, This is an ODE with integral terms on the RHS.

$$ \frac{db}{dt} = \int_t^{t1} f(b,h,t) dt $$

First thing is to solve the Integral on RHS, this can be done using numerical quadrature from Scipy.

Approximation of $\frac{du}{dt}$

For the Approximation of $du/dt$ within the integral , use a basic backward difference forumulation $u^{n} - u^{n-1}/\delta t$ and substitute the values for an better approximation of the derivative term. ( As long as your time step is sufficiently low, this appromimation should hold better, However you need to experiment with lot of time steps to see if the solution is diverging due to the numerical errors induced by this approximation)

Numerical Differentiation

The Pseudocode for numerical differentiation will look like

def RHS():
   ## RHS Term here
# initial condition
b0 = 1
  
# values of time
t = np.linspace(0,5)

# solving ODE
y = odeint(RHS, b0, t)

Numerical Integration

Further, you can perform the numerical integration within the RHS function using integrate method from scipy

from scipy import integrate

def terms_in_integral():
  ## all the terms goes here

def RHS():
   integrate.quadrature ( terms_in_integral , start_time ,  end_time )

This should address all the queries that you had mentioned in the question