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:
```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
```python
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
]
```