Skip to main content
Computational Science

Return to Question

Added solution
Source Link
Mr P
Mr P
  • 31
  • 3

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

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)

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)

/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)

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

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

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
    ]

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

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)

/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)

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

$ 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²$

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)

/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)

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

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
    ]
added 84 characters in body
Source Link
Mr P
Mr P
  • 31
  • 3

I'm trying to describe the movement of a freediver during freefall (falling without moving when buoyancy. 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

I'm trying to describe the movement of a freediver during freefall (falling without moving when buoyancy is negative). The Forces I'm considering are gravitation, buoyancy and drag. I've worked out this this formula for the sum of those forces

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

1. I added context and goal. 2. cleaned up the equations. 3. corrected mistake in implementation as pointed out in the answers
Source Link
Mr P
Mr P
  • 31
  • 3

This is the initial second order equation which I would likeI'm trying to solve for ddescribe the movement of a freediver during freefall (tfalling without moving when buoyancy is negative):. The Forces I'm considering are gravitation, buoyancy and drag. I've worked out this this formula for the sum of those forces

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

following a number of tutorials I rewroteFrom this as aI derived the following system of two first order equations asODEs

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

thenwith
$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)

      message: 'Required step size is less than spacing between numbers.'
     nfev: 626614
     njev: 0
      nlu: 0
      sol: None
   status: -1
  success: False
        t: array([0.       , 0.0060006, 0.0120012, 0.0180018, 0.0240024, 0.030003 ,
       0.0360036, 0.0420042, 0.0480048, 0.0540054, 0.060006 , 0.0660066,
       0.0720072, 0.0780078, 0.0840084, 0.090009 , 0.0960096, 0.1020102,
       0.1080108])
 t_events: None
        y: array([[20.        , 20.00611389, 20.01247168, 20.0191072 , 20.02606478,
        20.03339925, 20.04117593, 20.04947141, 20.05838901, 20.06805067,
        20.07863631, 20.09039598, 20.10366405, 20.11892066, 20.1369306 ,
        20.15915028, 20.18809947, 20.23005987, 20.30820872],
       [ 1.        ,  1.03820542,  1.08182954,  1.13179425,  1.18950952,
         1.25687343,  1.33627217,  1.43064426,  1.54448509,  1.6825812 ,
         1.85388098,  2.07263638,  2.36066427,  2.75061513,  3.2992731 ,
         4.1588621 ,  5.64191939,  8.821076  , 20.38529206]])
 y_events: 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

This is the initial second order equation which I would like to solve for d(t):

$ d = \frac{ 10 C_B V_{TLC} }{m d'' - C_R (d')² - C_B V_{diver} + mg} - 10 $

following a number of tutorials I rewrote this as a system of two first order equations as

$ d' = v $
$ v' = d'' = \frac{ C_R v² + \frac{ 10 C_B V_{TLC} }{(d + 10)} - mg + C_B V_{diver} }{m} $

then 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
    ]

# 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)

      message: 'Required step size is less than spacing between numbers.'
     nfev: 626
     njev: 0
      nlu: 0
      sol: None
   status: -1
  success: False
        t: array([0.       , 0.0060006, 0.0120012, 0.0180018, 0.0240024, 0.030003 ,
       0.0360036, 0.0420042, 0.0480048, 0.0540054, 0.060006 , 0.0660066,
       0.0720072, 0.0780078, 0.0840084, 0.090009 , 0.0960096, 0.1020102,
       0.1080108])
 t_events: None
        y: array([[20.        , 20.00611389, 20.01247168, 20.0191072 , 20.02606478,
        20.03339925, 20.04117593, 20.04947141, 20.05838901, 20.06805067,
        20.07863631, 20.09039598, 20.10366405, 20.11892066, 20.1369306 ,
        20.15915028, 20.18809947, 20.23005987, 20.30820872],
       [ 1.        ,  1.03820542,  1.08182954,  1.13179425,  1.18950952,
         1.25687343,  1.33627217,  1.43064426,  1.54448509,  1.6825812 ,
         1.85388098,  2.07263638,  2.36066427,  2.75061513,  3.2992731 ,
         4.1588621 ,  5.64191939,  8.821076  , 20.38529206]])
 y_events: None

I'm trying to describe the movement of a freediver during freefall (falling without moving when buoyancy is negative). 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{d + 10[m]}{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)

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

Source Link
Mr P
Mr P
  • 31
  • 3