# Can't solve second order differential equation with scipy

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

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
]

• Please clarify your specific problem or provide additional details to highlight exactly what you need. As it's currently written, it's hard to tell exactly what you're asking.
– Community Bot
Commented Feb 14, 2022 at 14:40

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.

• Thanks for your help! I corrected the mistake and added some context to the question. Now I will need some time to think about your answer :)
– Mr P
Commented Feb 16, 2022 at 19:04
• The friction force needs to be given as $-C_R|v|v$, so that it always works to slow down. Changing just that removes the singularity. /// Please check your model logic, using depth typically means that falling increases the depth, reversely with height, where falling decreases the coordinate. Your equations appear to mostly use the "height" concept, at least with the direction of the gravity component. Commented Feb 16, 2022 at 19:29
• That was the clue I needed. I didn't use the same coordinate system for depth and velocity. Now it works perfectly. Thanks a lot!!!
– Mr P
Commented Feb 17, 2022 at 1:38