Can't solve second order differential equation with scipy

Question

Most of my knowledge about numerically solving differential equations is long forgotten. Unfortunately I stumbled upon a physics problem where I need to do exactly that.

I'm trying to describe the movement of a freediver during freefall. This is the part of the dive when the diver has negative buoyancy and falls towards the targeted depth without any active movement.

The Forces I'm considering are gravitation, buoyancy and drag. I've worked out this this formula for the sum of those forces

$ F_{total} = F_B + F_G + F_D = C_B \left( V_{diver} + V_{tlc} \frac{10 [m]}{d + 10[m]} \right) - mg + C_R v²$

From this I derived the following system of first order ODEs

$ d' = v $
$ v' = \frac{ C_R v² + C_B V_{tlc} \frac{10 [m]}{d + 10[m]} + C_B V_{diver} -m g }{m}$

with
$d_0 = \text{initial depth}$
$v_0 = \text{initial velocity}$

From this I would like to be able to plot how the depth and velocity change over time depending of different initial conditions.

I tried to solve this in python:

import numpy as np
import math
import scipy as sp
from scipy.integrate import odeint
from scipy.integrate import solve_ivp

# Physical constants
rho = 1023.6               # kg/m³ density of saline water
g = 9.807                  # m/s²  gravitational acceleration on earth

# Assumptions about diver
V_diver = 0.062            # m³    volume of diver
V_tlc = 0.006              # m³    total lung capacity
m = 66                     # Kg    weight of diver
A = 0.07                   # m²    crossectional area of diver in diving direction              
C_D = 0.3                  # -     Drag coefficient

# Derived
C_B = rho * g              # Buoyency coefficient
C_R = 0.5 * rho * C_D * A  # Resistive coefficient

# equation
def dSdd(d, S):
    d, v = S
    return [
        v,
        (C_B * V_diver + 10 * C_B * V_tlc / ( d + 10 ) - m * g + C_R * v**2) / m   
        ]
# initial conditions
d_0 = 20
v_0 = 1
S_0 = [d_0, v_0]

# time interval
t = np.linspace(0, 60, 1000)
    
# solution odeint
odeint(dSdd, y0 = S_0, t=t, tfirst=True, full_output = 1) 

# solution solve_ivp
solve_ivp(dSdd, t_span=(0, max(t)), y0=S_0, t_eval=t)

from odeint I get the following output:

/home/marc/.cache/pypoetry/virtualenvs/data-science-6CF2GDM8-py3.9/lib/python3.9/site-packages/scipy/integrate/odepack.py:247: ODEintWarning: Excess work done on this call (perhaps wrong Dfun type). Run with full_output = 1 to get quantitative information.
    warnings.warn(warning_msg, ODEintWarning)

and solve_ivp returns:

message: 'Required step size is less than spacing between numbers.'
     nfev: 614
     njev: 0
      nlu: 0
      sol: None
   status: -1
  success: False
        t: array([0.        , 0.06006006, 0.12012012, 0.18018018, 0.24024024,
       0.3003003 , 0.36036036, 0.42042042, 0.48048048, 0.54054054,
       0.6006006 , 0.66066066, 0.72072072, 0.78078078, 0.84084084,
       0.9009009 , 0.96096096, 1.02102102, 1.08108108, 1.14114114,
       1.2012012 , 1.26126126, 1.32132132, 1.38138138, 1.44144144,
       1.5015015 , 1.56156156, 1.62162162, 1.68168168, 1.74174174,
       1.8018018 , 1.86186186, 1.92192192, 1.98198198, 2.04204204,
       2.1021021 , 2.16216216, 2.22222222, 2.28228228, 2.34234234,
       2.4024024 , 2.46246246, 2.52252252, 2.58258258, 2.64264264,
       2.7027027 , 2.76276276, 2.82282282, 2.88288288, 2.94294294,
       3.003003  , 3.06306306, 3.12312312, 3.18318318, 3.24324324,
       3.3033033 , 3.36336336, 3.42342342, 3.48348348, 3.54354354,
       3.6036036 , 3.66366366, 3.72372372, 3.78378378, 3.84384384,
       3.9039039 , 3.96396396, 4.02402402, 4.08408408, 4.14414414,
       4.2042042 , 4.26426426, 4.32432432, 4.38438438, 4.44444444,
       4.5045045 , 4.56456456, 4.62462462, 4.68468468, 4.74474474,
       4.8048048 , 4.86486486, 4.92492492, 4.98498498, 5.04504505,
       5.10510511, 5.16516517, 5.22522523, 5.28528529, 5.34534535,
       5.40540541, 5.46546547, 5.52552553, 5.58558559, 5.64564565,
       5.70570571, 5.76576577, 5.82582583, 5.88588589, 5.94594595,
       6.00600601, 6.06606607, 6.12612613, 6.18618619, 6.24624625,
       6.30630631, 6.36636637, 6.42642643, 6.48648649, 6.54654655,
       6.60660661, 6.66666667, 6.72672673, 6.78678679, 6.84684685,
       6.90690691, 6.96696697, 7.02702703, 7.08708709, 7.14714715,
       7.20720721, 7.26726727, 7.32732733, 7.38738739, 7.44744745,
       7.50750751, 7.56756757, 7.62762763, 7.68768769, 7.74774775,
       7.80780781, 7.86786787, 7.92792793, 7.98798799])
 t_events: None
        y: array([[ 20.        ,  20.06022325,  20.12077573,  20.18166178,
         20.24288651,  20.3044545 ,  20.36637012,  20.4286379 ,
         20.4912626 ,  20.55424917,  20.61760275,  20.6813287 ,
         20.74543256,  20.80992008,  20.87479721,  20.94007008,
         21.00574504,  21.07182864,  21.13832763,  21.20524893,
         21.27259971,  21.34038729,  21.40861921,  21.47730323,
         21.54644727,  21.61605948,  21.6861482 ,  21.75672196,
         21.82778949,  21.89935974,  21.97144185,  22.04404523,
         22.11719859,  22.19092128,  22.26521681,  22.34008972,
         22.41554563,  22.4915912 ,  22.56823418,  22.64548333,
         22.72334851,  22.80184062,  22.88097161,  22.9607545 ,
         23.04120335,  23.12233331,  23.20416057,  23.28670235,
         23.36997697,  23.4540038 ,  23.53880323,  23.62439676,
         23.71080692,  23.79805729,  23.88617252,  23.97517832,
         24.06510146,  24.15596974,  24.24781206,  24.34065834,
         24.43453958,  24.52948782,  24.62553619,  24.72271883,
         24.82107098,  24.92062892,  25.02142997,  25.12351255,
         25.2269161 ,  25.33168113,  25.43784921,  25.54546297,
         25.65456608,  25.7652033 ,  25.87742041,  25.99126427,
         26.1067828 ,  26.22402496,  26.34304079,  26.46388446,
         26.58669772,  26.71157337,  26.8385711 ,  26.96776102,
         27.09922374,  27.23305032,  27.36934229,  27.50821168,
         27.64978094,  27.79418303,  27.94156136,  28.09206981,
         28.24587274,  28.40314496,  28.56407178,  28.72884895,
         28.8976827 ,  29.07078974,  29.24839722,  29.43074279,
         29.61807456,  29.81083624,  30.0095311 ,  30.21431322,
         30.4254371 ,  30.64326674,  30.86827566,  31.1010469 ,
         31.34227302,  31.59275608,  31.85340768,  32.12524893,
         32.40941044,  32.70713236,  33.01976433,  33.34876555,
         33.69570468,  34.06225995,  34.4510292 ,  34.86534333,
         35.30887002,  35.78633731,  36.30353361,  36.86798076,
         37.48862955,  38.17688489,  38.95181674,  39.84035888,
         40.87817713,  42.12560177,  43.6927668 ,  45.79877743,
         49.03175348,  56.27649463],
       [  1.        ,   1.00544608,   1.01096335,   1.01655377,
          1.02221991,   1.02796387,   1.03378755,   1.03969295,
          1.04568219,   1.05175752,   1.05792131,   1.06417603,
          1.07052428,   1.07696878,   1.08351237,   1.09015801,
          1.09690877,   1.10376784,   1.11073853,   1.11782428,
          1.12502863,   1.13235525,   1.13980793,   1.14739057,
          1.15510719,   1.16296194,   1.17095908,   1.17910299,
          1.18739816,   1.19584922,   1.20446089,   1.21323815,
          1.22221116,   1.2313979 ,   1.2407953 ,   1.2504013 ,
          1.26021483,   1.27023586,   1.28046533,   1.2909052 ,
          1.30155845,   1.31242904,   1.32352196,   1.3348432 ,
          1.34639975,   1.35819961,   1.37025178,   1.38256629,
          1.39515414,   1.40802737,   1.421199  ,   1.43468308,
          1.44849466,   1.46264977,   1.47716548,   1.49205986,
          1.50735197,   1.52306189,   1.53921071,   1.55582052,
          1.5729144 ,   1.59051647,   1.60865184,   1.62734661,
          1.64662791,   1.66652388,   1.68706363,   1.70827733,
          1.7301961 ,   1.75285212,   1.77627853,   1.8005095 ,
          1.82558021,   1.85152684,   1.87838657,   1.9061976 ,
          1.93499912,   1.96483134,   1.99573547,   2.02775866,
          2.06108008,   2.09578494,   2.13190374,   2.16948422,
          2.20859136,   2.24930734,   2.29173157,   2.33598069,
          2.38218857,   2.4305063 ,   2.48110217,   2.53416174,
          2.58988776,   2.64850021,   2.7102363 ,   2.77535047,
          2.84411438,   2.9168169 ,   2.99376414,   3.07527943,
          3.16170334,   3.25408583,   3.3535738 ,   3.45987045,
          3.57297988,   3.69324192,   3.82133214,   3.95826186,
          4.10537813,   4.26436373,   4.4372372 ,   4.6263528 ,
          4.83440054,   5.06440617,   5.31973118,   5.60407278,
          5.92146395,   6.27627341,   6.67633158,   7.13167175,
          7.65433238,   8.26176668,   8.97684267,   9.83515678,
         10.87519331,  12.14316414,  13.75241444,  15.88570982,
         18.80474186,  23.02438119,  29.73016354,  41.82752868,
         70.82032029, 230.53580339]])
 y_events: None

Can anybody point me in the right direction to solve this?

Thank you

Edit
The problem was the use of different coordinate systems for velocity and forces versus depth as pointed out by Lutz Lehmann in the comments.

Changing the equations in the following way made things work perfectly:

$ F_{total} = F_B + F_G + F_D = m g - C_B \left( V_{diver} + V_{tlc} \frac{10 [m]}{d + 10[m]} \right) - C_R v²$

and

$ d' = v $
$ v' = d'' = \frac{ - C_R v² - C_B V_{tlc} \frac{10 [m]}{d + 10[m]} - C_B V_{diver} + m g }{m}$

with
$d_0 = \text{initial depth}$
$v_0 = \text{initial velocity}$

and

def dSdd(d, S):
    d, v = S
    return [
        v,
        - C_R * v ** 2  -  C_B * V_tlc * 10 / ( d + 10)  -  C_B * V_diver  +  m * g / m
    ]

Answer 1

Why did you omit the division by $m$ in the $v$ derivative?

Apparently the dynamic of $$\dot v = ...+C_Rv^2$$ explodes, and this singularity leads to local errors that can not be compensated by reducing the internal step size above the minimal admissible step size. The second error message of solve_ivp literally means that $t+dt$ is not greater than $t$ in floating-point arithmetic.

In the first error message of odeint something similar happens in that the internal step size is so small that the number of necessary steps to reach the next output time becomes larger than the bound in max_step. You can raise that number, but if the fundamental problem is the singularity, then the error is just shifted to some slightly larger time. You might, with a sufficiently large step number, get a different error in that floating-point overflow occurs.

---

Correcting the first omission

def dSdd(d, S):
    d, v = S
    F = C_B * V_diver + 10 * C_B * V_tlc / ( d + 10 ) - m * g + C_R * v**2
    return  v, F/m

still gives raise to the same error messages. One can get the internal steps of solve_ivp by not giving the t_eval argument. Then the last reported time is 8.014644648576457 with a value of v of 1.32492524e+14, clearly demonstrating that the quadratic feed-back dominates the other terms of the dynamic.

---

The distance to the singularity can be estimated. For large $v$ the dynamic $\dot v=av^2$ dominates. This is exactly solvable as $$ v(t)=\frac{v(t_1)}{1-av(t_1)(t-t_1)} $$ so that the singularity can be found at or around $$ t_\ast=t_1+\frac1{av(t_1)}=t_1+\frac{m}{C_Rv(t_1)} $$ For the computed values this give a time offset $t_*-t_1$ of 4.634821e-14.