计算科学 - 求方程的根（在某些点上发散） - 吾爱随笔录

求方程的根（在某些点上发散）

计算科学优化 Python

2021-11-29 01:34:46

我正在尝试解决这个等式 $v_{Sl}$ 范围中 $[0.001,1]$ 获得的值 $\delta_l$ . 我尝试过Fsolve使用 python，它给出的结果是 $v_{Sl}$ 介于 $[0.001,0.1]$ 当我初步猜测 $\delta_l$ 这接近解决方案。但是，当我改变初始猜测时，我没有获得整个范围的任何解决方案。

我还使用提供的所有方法尝试了 scipy root，但没有获得任何解决方案。当我尝试使用进化算法（例如遗传算法）时，它似乎给出了结果，但受到所提供解决方案范围的限制。例如，如果我选择解决方案的范围 $\delta_l$ 来自 $[0.001,0.065]$ ，我得到了我真正正在寻找的解决方案，但如果我扩大/缩小范围，我会获得不同的解决方案。理想情况下，我会给出一系列可能属于的解决方案，比如说 $[0.0001,0.2]$ 并且该算法将给出与我给出一个范围时相同的结果 $[0.001,0.065]$ .

由此我了解到，该算法没有收敛到方程的解（ $0$ ) 我需要一个更好的工具来使其收敛。有没有更好的算法来解决这个问题？

代码和数字

事实上，我希望获得的最终解决方案是 $v_{sg}$ . 我将首先获得 $\delta_l$ 求解上述方程的值，然后计算 $v_{sg}$ 从以下等式。

下图是一个参考，我需要用我的代码重现相同的图。

这是我生成的图表。

我写的代码如下。主要问题是解决方案的范围 $\delta_l$ 我需要指定以获得所需的结果，对于我的情况，我把 $[0.001,0.065]$ ，但我不应该缩小范围，因为我不知道另一种情况下的解决方案（不同的数据）。相反，我会说一个可接受的范围 $[0.0001,0.2]$ 并且仍然期望得到类似的结果，但事实并非如此。

情节 $\delta_l$ vs 我正在解决的功能：

import pygad
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)
    ga_instance = pygad.GA(num_generations=num_generations,
                       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')

我尝试使用以下代码使用 scipy.root 的方法：

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[0])
    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[0]
        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")

1个回答

你的方程基本上用系数读取 $a_1$ 和 $a_2(v_{Sl})$

a_{1} ((1 - 2 x)^{2} - 2 (x - x^{2})) + a_{2} \frac{x - x^{2} + (1 - 2 x)^{2}}{(x - x^{2})^{3}} = 0

$a_1 ((1-2x)^2-2(x-x^2))+a_2 \frac{x-x^2+(1-2x)^2}{(x-x^2)^3} =0$ 两边乘以分母后

(x - x^{2})^{3}

$(x-x^2)^3$ 您将获得一个 8 次多项式。通过计算其伴随矩阵的特征值可以很容易地找到它的根，例如numpy.polynomials.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)}')

其它你可能感兴趣的问题

上一篇氢的薛定谔方程的数值解下一篇如何计算置换矩阵的行列式