Scipy fsolve未收敛到正确值

imzjd6km 于 2023-04-12 发布在其他

关注(0)|答案(1)|浏览(197)

我试图求解一个简单微分方程的初始值（dh/dt，试图在t=0时找到h）。然而，使用fsolve，迭代的值是：

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

不管我给予P、T、kf或kp什么值，代码迭代都遵循相同的模式。运行代码也会给出错误警告。所以我的问题是，问题是在代码本身吗？

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

scipy

来源：https://stackoverflow.com/questions/75924877/scipy-fsolve-doesnt-converge-to-correct-value

1条答案

按热度按时间

k2arahey1#

在IIUC中，有很多东西需要改变。这段代码运行并产生输出：

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.")

几个音符。。

我重构了一些东西，包括删除一些不需要的导入和一些不使用的变量
将x0的初始值更改为[0.5*R]，这是fsolve使用的求根算法的更好的初始猜测。
修改了Aire和dv_dh函数，以处理h〉= R的情况，这可能会导致无效值传递给np.sqrt。

赞(0）回复(0）举报 2023-04-12

我来回答

Scipy fsolve未收敛到正确值

1条答案

相关问题

热门标签

最新问答