Scipy fsolve doesn't converge to correct value

I'm trying to solve for the initial value of a simple differential equation (dh/dt, trying to find h at t=0). However, using fsolve, the value for the iterations are:

1
1.0
1.0
1.0000000149011612
101.0
nan
1.0000000149011612
101.0
nan
nan
nan
nan
nan
nan
nan

The code iterations follow the same pattern no matter what values I give P,T,kf or kp. Running the code also gives the error warning. So my question is is the issue is in the code itself?

# -*- coding: utf-8 -*-
"""
Created on Mon Apr  3 16:46:15 2023

@author: houle
"""

import numpy as np
import math
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
from scipy.optimize import fsolve

class parametre_antoine():
    A = 8.13484
    B = 1662.48
    C = 238.131
    mmhg_atm = 760
    
prm_antoine = parametre_antoine() 

rho = 0.8 #kg/L
Tin = 110 #C   
R = 1
kp = -200
kf = 300
M = 46.068/1000 #kg/mol
L = 5
M = 46.068 #g/mol
Vtot = np.pi*R**2*L

theta = [rho,R,kp,kf,M,L,prm_antoine,Vtot]

Fin = 20000
u = [Fin]

P = 3.3
T = 300

def Aire(R,L,h):
    return 2*L*np.sqrt(R**2-(R-h)**2)

def dv_dh(R,L,h):
    dVdh = R/np.sqrt(1-(h/R-1)**2)+L*np.sqrt((2*R-h)*h)-(L*(2*R-2*h)*(R-h))/(2*np.sqrt((2*R-h)*h))
    return dVdh

def dh_dt(h,theta,u,T,P):
    [rho,R,kp,kf,M,L,prm_antoine,Vtot] = theta
    [Fin] = u
    dhdt =(Fin- kp*(psat(prm_antoine,T)-P)*Aire(R,L,h)-kf*h**0.5)/rho/ dv_dh(R,L,h)
    print(h)
    return dhdt

def psat(prm,T):
    #T en Celcius
    #Valide entre -114 et 243C
    #Retourne la tension de vapeur en atm
    p_mmhg = 10**(prm.A-prm.B/(T+prm.C))
    p_atm = p_mmhg / prm.mmhg_atm
    return p_atm

x0 = [1]
u0 = [Fin]
x0 = fsolve(dh_dt, x0, args=(theta, u0,T,P))

RuntimeWarning: invalid value encountered in sqrt
  return 2*L*np.sqrt(R**2-(R-h)**2)

IIUC, there was a bevy of things that needed to change. This code runs and produces an output:

import numpy as np
from scipy.optimize import fsolve
# Define the parameters
rho = 0.8 # kg/L
R = 1 # m
L = 5 # m
kp = -200
kf = 300
M = 46.068 # g/mol
prm_A = 8.13484
prm_B = 1662.48
prm_C = 238.131
mmhg_atm = 760
Vtot = np.pi*R**2*L
P = 3.3 # atm
T = 300 # K
Fin = 20000 # L/h
# Define the functions
def psat(T):
p_mmhg = 10**(prm_A-prm_B/(T+prm_C))
p_atm = p_mmhg / mmhg_atm
return p_atm
def Aire(R,L,h):
if h >= R:
return 0
else:
return 2*L*np.sqrt(R**2-(R-h)**2)
def dv_dh(R,L,h):
if h >= R:
return 0
else:
return R/np.sqrt(1-(h/R-1)**2)+L*np.sqrt((2*R-h)*h)-(L*(2*R-2*h)*(R-h))/(2*np.sqrt((2*R-h)*h))
def dh_dt(h):
return (Fin- kp*(psat(T)-P)*Aire(R,L,h)-kf*h**0.5)/rho/ dv_dh(R,L,h)
# Solve for the initial value of h
x0 = [0.5*R]
h0 = fsolve(dh_dt, x0)[0]
print(f"The initial value of h is {h0:.4f} m.")

A few notes...

• I refactored a few things, including removing some imports that weren't needed and a few variables that weren't being used
• Changed the initial value of x0 to [0.5*R], which is a better initial guess for the root-finding algorithm used by fsolve.
• Modified the Aire and dv_dh functions to handle cases where h >= R, which can cause invalid values to be passed to np.sqrt.