solving ODE $y'' = -10 - abs(y')\times y'$ using python

Question

Need to solve: $y'' = -10 - |y'|\times y'$ which it says to be modelling a fall of a ball.

I wrote a code for it which goes:

import numpy as np
import matplotlib.pyplot as plt

# DE definition
def get_y2dot(ydot):
    return -10 - abs(ydot)*ydot

# solution definition
def y(t, timestep, ydot0, y0):
    ydot = ydot0
    y = y0
    for time in np.arange(0,t,timestep):
        y2dot = get_y2dot(ydot)
        ydot += y2dot * timestep
        y += ydot * timestep
    return y

# ICs
Ydot0 = 0 # initial velocity
Y0 = 10 # initial height
timestep = 0.01

time = np.arange(0,10, timestep)
ys = [y(i, timestep, Ydot0, Y0) for i in time]

plt.plot(time,ys)
plt.xlabel("t")
plt.ylabel("y")

But I get negative $y$ which is wrong but I'm not sure where I'm going wrong in the code.

I guess more importantly, I don't know how to solve the DE. I tried: $$let \ y'=Y, \ y''=Y'$$ I didn't know how to deal with the absolute value so I said: $$ for \ Y>0: Y'+Y^2 = -10 $$ but I'm not sure what to do with the $Y^2$.

This post almost answered my question: Solving second-order nonlinear nonhomogeneous differential equation But I don't know how they solve $\int y''y'dx$. My thought process was: $$ \frac{d^2y}{dx^2} \frac{dy}{dx} + a. y^3 \frac{dy}{dx} = b \frac{dy}{dx} \rightarrow$$ $$ \require{cancel} \int \frac{d^2y}{dx^2} \frac{dy}{\cancel{dx}}\cancel{dx} + a \int y^3 \frac{dy}{\cancel{dx}}\cancel{dx} = b \int \frac{dy}{\cancel{dx}}\cancel{dx} $$

Answer 1

I don't know why you think $y<0$ is wrong. The ball accelerates downward, approaching its terminal velocity of $y_{\infty}^{\prime}=-\sqrt{10}$ so its height will be a little above $y\approx y_0-\sqrt{10}\,t$. Let's solve the equation analytically: $$y^{\prime\prime}=-10+y^{\prime\,2}$$ For this problem I have dispensed with the absolute value beause we know the ball will be falling. Let $y=y^{\prime}$. Then $y^{\prime\prime}=\frac{dv}{dt}=\frac{dv}{dy}\frac{dy}{dt}=v\frac{dv}{dy}$ so $$v\frac{dv}{dy}=v^2-10$$ $$\int\frac{v\,dv}{v^2-10}=\frac12\ln\lvert v^2-10\rvert=\int dy=y+C_1$$ $$v^2-10=\pm e^{2C_1}e^{2y}$$ When $t=0$, $v=0$, $y=y_0$, so $\pm e^{2C_1}=-10e^{-2y_0}$. It is falling so $$\frac{dy}{dt}=v=-\sqrt{v^2}=-\sqrt{10\left(1-e^{2(y-y_0)}\right)}$$ If we make the substitution $u=e^{-(y-y_0)}$ then $du=-e^{-(y-y_0)}dy=-u\,dy$ and $$\int\frac{dy}{\sqrt{1-e^{2(y-y_0)}}}=\int\frac{-\frac{du}u}{\sqrt{1-\frac1{u^2}}}=-\int\frac{du}{\sqrt{u^2-1}}$$ Now let $u=\cosh\theta$, $\sqrt{u^2-1}=\sinh\theta$, $du=\sinh\theta\,d\theta$ so $$\begin{align}\int\frac{dy}{\sqrt{1-e^{2(y-y_0)}}}&=-\int d\theta=-\theta=-\cosh^{-1}u=-\cosh^{-1}\left(e^{-(y-y_0)}\right)\\ &=\int-\sqrt{10}\,dt=-\sqrt{10}\,t+C_2\end{align}$$ At $t=0$, $y=y_0$, $\cosh^{-1}\left(e^0\right)=\cosh^{-1}(1)=0$ so $C_2=0$ and $$\cosh^{-1}\left(e^{-(y-y_0)}\right)=\sqrt{10}\,t$$ Solving for $y$, $$y=y_0-\ln\left(\cosh\left(\sqrt{10}\,t\right)\right)$$ In your python code you seem to be solving the whole differential equation from scratch to get the solution at each time point. This requires $\frac12n(n-1)$ derivative evaluations to find the solutions at each of $n$ time points. Instead I recommend putting each value of $y$ into the array to be used for plotting during a single sweep through all time points thus requiring only $n-1$ derivative evaluations. Here is my Matlab code which attempts to reproduce your solution and also plots the exact solution:

% ball.m

% DE definition
get_y2dot = @(ydot) -10-abs(ydot).*ydot;

% ICs
Ydot0 = 0;
Y0 = 10;
timestep = 0.01;

time = [0:timestep:10];
ys = zeros(size(time));
% Numeriacal solution
y = Y0;
ydot = Ydot0;
ys(1) = Y0;
for k = 2:length(time),
    ydot = ydot+get_y2dot(ydot)*timestep;
    y = y+ydot*timestep;
    ys(k) = y;
end
% Exact solution
y_exact = Y0-log(cosh(sqrt(10)*time));
plot(time,ys,'r-',time,y_exact,'b--');
title('Falling Ball');
xlabel('t');
ylabel('y');
legend('Numerical','Exact','Location','best');

And my plot: