Finding the root of an equation (divergence at some points) I am trying to solve this equation for $$v_{Sl}$$ in the range $$[0.001,1]$$ to obtain the values of $$\delta_l$$. I have tried Fsolve with python, which gives results where the values of $$v_{Sl}$$ are between $$[0.001,0.1]$$ when I put an initial guess of $$\delta_l$$ that is close to the solution. But, when I vary the initial guess I don't obtain any solution for the whole range.

I also tried scipy root with all the methods provided and I did not obtain any solution. When I tried with an Evolutionary Algorithm, such as Genetic Algorithm, it seems to give results, but with the limitation of the range of solution provided. For instance, if I choose the range of solution for $$\delta_l$$ to be from $$[0.001,0.065]$$, I get the solution I am actually looking for, but if I widen/narrow the range I obtain a different solution. Ideally, I would give a range of solutions that may fall in let's say $$[0.0001,0.2]$$ and the algorithm would give the same results as when I give a range of $$[0.001,0.065]$$.

From this I understand that the algorithm is not converging toward the solution of the equation ($$0$$) and I need a better tool to make it converge. Is there a better algorithm to solve this problem?

Code and Figures

As a matter of fact, the final solution I would be looking to obtain is $$v_{sg}$$. I will first obtain the $$\delta_l$$ value from solving the mentioned equation then calculate $$v_{sg}$$ from the following equation. The following graph is a reference which I need to reproduce the same graph with my code. Here is the graph that I generated. The code I wrote is below. The main problem is the range of solutions for $$\delta_l$$ that I need to specify to get the desired results, for my case here I put $$[0.001,0.065]$$, but I should not narrow it as I will not know the solution in another case (different Data). Rather I would put an acceptable range let's say $$[0.0001,0.2]$$ and would still expect to get similar results, which is not the case.

Plot of $$\delta_l$$ vs function I am solving: import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import fsolve
import math
###### Constants start here
D= 0.06  ### Diameter (m)
ROL= 997.9 #### Density liquid  Kg/m3
ROG= 1.2  ##### density gas    kg/m3
UL= 1.1*(10**(-3))### viscosit Pa.S
UG= 0.000018  ### viscosit Pa.S
G= 9.81   # gravitational acceleration  m/s2
sigma= 0.06    #### Interfacial tension  N/m
alpha= 20   #### inclination from horizontal   Degrees
SL=3.14*D
### constants ends

##### list of the values of the variable that we are solving the equation with each time
list_of_liquid_velocities_VSL= np.linspace(0.001,1,80) #va
### ends here

### lists to store the output values in
list_of_derivatives=[]
list_of_interfacial_tension=[]
list_of_interfacial_tension_gas=[]
list_of_VSG=[]   ###### main values we are looking for
### ends

##################### Main algorithm  ##############

for VSL in list_of_liquid_velocities_VSL:
AL=3.14*(D**2)*(0.01-0.01**2)

DL=4*AL/SL
A=3.14*(0.5*D)**2
VL=VSL*A/AL
if (D*ROL*VSL/(UL))<2300:
n=1
CL=16
m=0.2
CG=0.046
else:
n=0.2
CL=0.046
m=0.2
CG=0.046

desired_output = 0

def fitness_func(x,x_idx):
derivative = ((G*(ROL-ROG)*D*math.sin(alpha*3.14/180)*((1-2*x)**2-2*(x-x**2)))-((CL/16)*ROL*((D*ROL/UL)**(-n))*(VSL**(2-n))*(((x-x**2)+(1-2*x)**2)/(x-x**2)**3)))  ### main equation I am solving

fitness = (1/np.abs(derivative - desired_output))
return fitness

fitness_function = fitness_func

num_generations = 100
num_parents_mating = 4

sol_per_pop = 8
num_genes = 1
init_range_low = 0.001
init_range_high = 0.101

parent_selection_type = "sss"
keep_parents = 1

crossover_type = "single_point"

mutation_type = "random"
mutation_percent_genes = 20

gene_space = [{'low': 0.001, 'high': 0.065}] ### Range I am giving for the algorithm to look for the solution in (try to vary it and see the change in the final plot)
num_parents_mating=num_parents_mating,
fitness_func=fitness_function,
sol_per_pop=sol_per_pop,
num_genes=num_genes,
init_range_low=init_range_low,
init_range_high=init_range_high,
parent_selection_type=parent_selection_type,
keep_parents=keep_parents,
crossover_type=crossover_type,
mutation_type=mutation_type,
mutation_percent_genes=mutation_percent_genes,
gene_space= gene_space)

ga_instance.run()
x, solution_fitness, solution_idx = ga_instance.best_solution()

### x is the solution of the function above
derivative = ((G*(ROL-ROG)*D*np.sin(alpha*3.14/180)*((1-2*x)**2-2*(x-x**2)))-((CL/16)*ROL*((D*ROL/UL)**(-n))*(VSL**(2-n))*(((x-x**2)+(1-2*x)**2)/(x-x**2)**3))) ### recalculation of the equation above to double confirm (which should should be equal to 0)

list_of_derivatives.append(abs(derivative))

interfacial_tenstion=((G*(ROL-ROG)*D*np.sin(alpha*3.14/180)*(x-x**2)*(1-2*x))+((CL/32)*ROL*(D*ROL/UL)**(-n)*VSL**(2-n)*((1-2*x)/(x-x**2)**2)))  ### should be equal to interfacial_tension_gas

VSG  =  (interfacial_tenstion*((1-2*x)**4) /((0.5)*(CG)*((D*ROG/UG)**(-m))*(1+300*x)*(ROG)))**(1/(2-m)) ### VSG the main output which will be plotted against VSL in a loglog plot

list_of_interfacial_tension.append(interfacial_tenstion)

interfacial_tension_gas= ((((VSG)**(2-m))*(0.5)*(CG)*((D*ROG/UG)**(-m))*(1+300*x)*(ROG))/((1-2*x)**4))/(G*(ROL-ROG)*D)  ### this should equal to interfacial_tension
list_of_interfacial_tension_gas.append(interfacial_tension_gas)

list_of_VSG.append(VSG)

VSLexperiment= [0.02,0.06,0.1]
VSGexperiment=[13.94,14.54,19.09]

plt.plot()
plt.xscale("log")
plt.yscale("log")
plt.xlim(1,100)
plt.ylim(0.001,1)
plt.xlabel("VSG", fontsize=9)
plt.ylabel("VSL",fontsize=9)
plt.grid(True, which="both", ls=":", color='0.002', linewidth=0.4)
plt.scatter(VSGexperiment,VSLexperiment, color="Black", linewidth="0.05", label="VsgExperiment")
plt.plot (list_of_VSG,list_of_liquid_velocities_VSL, label='VBarnea', color='blue')

I tried with approach of scipy.root with the following codes:

from numpy import sin, linspace
from numpy.polynomial import Polynomial
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import fsolve
import math

###### Constants start here
D= 0.06  ### Diameter (m)
ROL= 997.9 #### Density liquid  Kg/m3
ROG= 1.2  ##### density gas    kg/m3
UL= 1.1*(10**(-3))### viscosit Pa.S
UG= 0.000018  ### viscosit Pa.S
G= 9.81   # gravitational acceleration  m/s2
sigma= 0.06    #### Interfacial tension  N/m
alpha= 20   #### inclination from horizontal   Degrees
SL=3.14*D
### constants ends

##### list of the values of the variable that we are solving the equation with each time
list_of_liquid_velocities_VSL= np.linspace(0.001,1,80) #va
### ends here

### lists to store the output values in
list_of_F_Function=[]
list_of_interfacial_tension=[]
list_of_interfacial_tension_gas=[]
list_of_x_solutions=[]
list_of_VSG=[]   ###### main values we are looking for
rooots=[]
### ends

for VSL in list_of_liquid_velocities_VSL:
def fitness_roots(VSL):
if (D*ROL*VSL/(UL)) < 2300:
n = 1
CL = 16
else:
n = 0.2
CL = 0.046

a1 = G*(ROL-ROG)*D*sin(alpha*3.14/180)
a2 = -(CL/16)*ROL*(D*ROL/UL)**(-n)*VSL**(2-n)
#result = a1*((1-2*x)**2-2*(x-x**2)) + a2*((x-x**2)+(1-2*x)**2)/(x-x**2)**3
#coeffs =[ a2, -3*a2, 3*a2, a1, -9*a1, 27*a1, -37*a1, 24*a1, -6*a1] # from mathematica cause I am lazy
coeffs=[-6*a1, +24*a1, -37*a1, 27*a1, -9*a1, a1, 3*a2, -3*a2, a2]
p = Polynomial(coeffs)
r = p.roots()
# return r # if you want all roots
return r [abs(r.imag) < 1e-15].real # if you only want real roots

rooots.append(fitness_roots(VSL))
x=(rooots)
AL=3.14*(D**2)*(x-x**2)

DL=4*AL/SL
A=3.14*(0.5*D)**2
VL=VSL*A/AL
if (D*ROL*VSL/(UL))<2300:
n=1
CL=16
m=0.2
CG=0.046
else:
n=0.2
CL=0.046
m=0.2
CG=0.046

F = ((G*(ROL-ROG)*D*np.sin(alpha*3.14/180)*((1-2*x)**2-2*(x-x**2)))-((CL/16)*ROL*((D*ROL/UL)**(-n))*(VSL**(2-n))*(((x-x**2)+(1-2*x)**2)/(x-x**2)**3))) ### should be equal to 0
list_of_F_Function.append((F))

interfacial_tenstion=((G*(ROL-ROG)*D*np.sin(alpha*3.14/180)*(x-x**2)*(1-2*x))+((CL/32)*ROL*(D*ROL/UL)**(-n)*VSL**(2-n)*((1-2*x)/(x-x**2)**2)))  ### should be equal to interfacial_tension_gas

VSG  =  (interfacial_tenstion*((1-2*x)**4) /((0.5)*(CG)*((D*ROG/UG)**(-m))*(1+300*x)*(ROG)))**(1/(2-m))

list_of_VSG.append(VSG)
list_of_x_solutions.append(x)
print('VSL: {VSL}\t {fitness_roots(VSL)}')

# list of the values of the variable that we are solving the equation with each time

for VSL in list_of_liquid_velocities_VSL:
print(f'VSL: {VSL}\t {fitness_roots(VSL)}')

plt.plot()
plt.xscale("log")
plt.yscale("log")
plt.xlim(0.01,100)
plt.ylim(0.0001,100)
plt.xlabel("VSG", fontsize=9)
plt.ylabel("VSL",fontsize=9)
plt.grid(True, which="both", ls=":", color='0.002', linewidth=0.4)

plt.plot (list_of_VSG,list_of_liquid_velocities_VSL, label='VBarnea', color='blue')

or this code

`python
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import fsolve, root, minimize
import math

from scipy.interpolate import make_interp_spline, BSpline

###### Constants start here
D= 0.06  ### Diameter (m)
ROL= 997.9 #### Density liquid  Kg/m3
ROG= 1.2  ##### density gas    kg/m3
UL= 1.1*(10**(-3))### viscosit Pa.S
UG= 0.000018  ### viscosit Pa.S
G= 9.81   # gravitational acceleration  m/s2
sigma= 0.06    #### Interfacial tension  N/m
alpha= 20   #### inclination from horizontal   Degrees
SL=3.14*D
### constants ends

##### list of the values of the variable that we are solving the equation with each time
list_of_liquid_velocities_VSL= np.linspace(0.001,1,80) #va
### ends here

### lists to store the output values in
list_of_F_Function=[]
list_of_interfacial_tension=[]
list_of_interfacial_tension_gas=[]
list_of_x_solutions=[]
list_of_VSG=[]   ###### main values we are looking for
### ends

VSLexperiment= [0.02,0.06,0.1]
VSGexperiment=[13.94,14.54,19.09]

#################  Barnea Model  ##################

for VSL in list_of_liquid_velocities_VSL:
def myF(z):
x=z
if (D*ROL*VSL/(UL)) < 2300:
n = 1
CL = 16
else:
n = 0.2
CL = 0.046

a1 = G*(ROL-ROG)*D*np.sin(alpha*3.14/180)
a2 = -(CL/16)*ROL*(D*ROL/UL)**(-n)*VSL**(2-n)
#result = a1*((1-2*x)**2-2*(x-x**2)) + a2*((x-x**2)+(1-2*x)**2)/(x-x**2)**3
#coeffs =[ a2, -3*a2, 3*a2, a1, -9*a1, 27*a1, -37*a1, 24*a1, -6*a1] # from mathematica cause I am lazy
F =-6*a1*x**8 +24*a1*x**7 -37*a1*x**6 + 27*a1*x**5 -9*a1*x**7 + a1*x**3+ 3*a2*x**2 -3*a2*x + a2
return F

solution = root(myF, [0.05], method='krylov')
x=solution.x
AL=3.14*(D**2)*(x-x**2)
DL=4*AL/SL
A=3.14*(0.5*D)**2
VL=VSL*A/AL

if (D*ROL*VSL/(UL))<2300:
n=1
CL=16
m=0.2
CG=0.046
else:
n=0.2
CL=0.046
m=0.2
CG=0.046

list_of_x_solutions.append(x)

AL=3.14*(D**2)*(x-x**2)
SL=3.14*D
DL=4*AL/SL
#DG=D-a
#AG=3.14*(0.5*D-a)**2
A=3.14*(0.5*D)**2
VL=VSL*A/AL
#VG=VSG*A/AG

F= ((G*(ROL-ROG)*D*np.sin(alpha*3.14/180)*((1-2*x)**2-2*(x-x**2)))-((CL/16)*ROL*((D*ROL/UL)**(-n))*(VSL**(2-n))*(((x-x**2)+(1-2*x)**2)/(x-x**2)**3)))
list_of_F_Function.append(F)

if (DL*ROL*VL/(UL))<2300:
n=1
CL=16
m=0.2
CG=0.046
else:
n=0.2
CL=0.046
m=0.2
CG=0.046

for VSL,x in zip(list_of_liquid_velocities_VSL,list_of_x_solutions):

interfacial_tension= ((G*(ROL-ROG)*D*math.sin(alpha*3.14/180)*(x-x**2)*(1-2*x))+((CL/32)*ROL*((D*ROL/UL)**(-n))*(VSL**(2-n))*((1-2*x)/(x-x**2)**2)))

VSG= (interfacial_tension*((1-2*x)**4) /((0.5)*(CG)*((D*ROG/UG)**(-m))*(1+300*x)*(ROG)))**(1/(2-m))
interfacial_tension_gas= ((((VSG)**(2-m))*(0.5)*(CG)*((D*ROG/UG)**(-m))*(1+300*x)*(ROG))/((1-2*x)**4))

list_of_interfacial_tension_gas.append(interfacial_tension_gas)
list_of_VSG.append(VSG)
list_of_interfacial_tension.append(interfacial_tension)

plt.plot()
plt.xscale("log")
plt.yscale("log")
plt.xlim(1,100)
plt.ylim(0.001,1)
plt.xlabel("VSG", fontsize=9)
plt.ylabel("VL",fontsize=9)
plt.grid(True, which="both", ls=":", color='0.002', linewidth=0.4)

plt.plot (list_of_VSG,list_of_liquid_velocities_VSL, label='VBarnea')
plt.scatter(VSGexperiment,VSLexperiment, color="Black", linewidth="0.05", label="VsgExperiment")
• Can you plot the function on the left for different values of VSL? It should have values where it cross the x-axis (well, the $\delta_l$ axis) and it should be relatively easy to see what the problem may be. Nov 5 '21 at 2:27
• @WolfgangBangerth I edited the post probably now is clear enough. Nov 5 '21 at 3:19
• The edit didn't make it better. I'm now totally confused which equation is the problem. Please make it clear which equation you want to, but can't solve, and what all of the variables in that are. For example, in both of the equations you currently have in the question, I have no idea what variables have fixed values and which ones are free that you try to find values for that solve the equation. Nov 5 '21 at 4:08
• the first equation is the one I am trying to solve by getting the value of delta_l that makes the equation equal to 0. In that equation I have constants such as (g , D, ROL, ROG), each time I will vary VSL and look for delta_l values that make the equation equal to 0 (I will carry VSL 80 times with this np.linspace(0.001,1,80) to obtain 80 values of delta_l that makes the equation equal to 0. The second equation gets that value of delta_l to solve for VSG with the other parameters fixed. the problem is in the first equation which is defined in the code as (def fitness_func...... derivative). Nov 5 '21 at 14:42
• I tried to reformat your question, but there is a lot more that could be done to clarify your post. For one, you can convert the images of equations into Mathjax, as I have done for several variables you wrote in your question. The code is also very hard read; you could try to change the formatting (reasonable number of spaces between lines, define intermediate quantities rather than massive multiline variables). Nov 16 '21 at 20:35

Your equation essentially reads with coefficients $$a_1$$ and $$a_2(v_{Sl})$$ $$a_1 ((1-2x)^2-2(x-x^2))+a_2 \frac{x-x^2+(1-2x)^2}{(x-x^2)^3} =0$$ After multiplying both sides with the denominator $$(x-x^2)^3$$ you obtain a polynomial of degree 8. Its roots can be easily found by calculating the eigenvalues of its companion matrix, with e.g. numpy.polynomials.roots.

So there is no need for any fancy genetic algorithm.

This is done in the attached code. I didn't get from your question if you have any more criteria to select from multiple roots or what you do if there are no roots.

from numpy import sin, linspace
from numpy.polynomial import Polynomial
# Constants
D = 0.06  # Diameter (m)
ROL = 997.9  # Density liquid  Kg/m3
ROG = 1.2  # density gas    kg/m3
UL = 1.1*(10**(-3))  # viscosit Pa.S
UG = 0.000018  # viscosit Pa.S
G = 9.81   # gravitational acceleration  m/s2
sigma = 0.06  # Interfacial tension  N/m
alpha = 20  # inclination from horizontal   Degrees
SL = 3.14*D

def fitness_roots(VSL):
if (D*ROL*VSL/(UL)) < 2300:
n = 1
CL = 16

else:
n = 0.2
CL = 0.046

a1 = G*(ROL-ROG)*D*sin(alpha*3.14/180)
a2 = -(CL/16)*ROL*(D*ROL/UL)**(-n)*VSL**(2-n)
#result = a1*((1-2*x)**2-2*(x-x**2)) + a2*((x-x**2)+(1-2*x)**2)/(x-x**2)**3
coeffs =[ a2, -3*a2, 3*a2, a1, -9*a1, 27*a1, -37*a1, 24*a1, -6*a1] # from mathematica cause I am lazy

p = Polynomial(coeffs)
r = p.roots()
# return r # if you want all roots
return r[abs(r.imag) < 1e-15].real # if you only want real roots

if __name__ == '__main__':

VSLexperiment = [0.02, 0.06, 0.1]
VSGexperiment = [13.94, 14.54, 19.09]

for VSL in VSLexperiment:
print(f'VSL: {VSL}\t {fitness_roots(VSL)}')

# list of the values of the variable that we are solving the equation with each time
list_of_liquid_velocities_VSL = linspace(0.001, 1, 80)  # va

for VSL in list_of_liquid_velocities_VSL:
print(f'VSL: {VSL}\t {fitness_roots(VSL)}')
• I did what you suggested and here is the code. However, when I recalculate the function with the roots it doesn't give me 0 at all. the code was added to the post. Nov 6 '21 at 4:50
• I don't get, why you change the order of coefficients. Of course the now calculated roots are different.
– Bort
Nov 6 '21 at 8:51