Scipy fsolve无法收敛到正确的值。

huangapple
go评论103阅读模式
英文:

Scipy fsolve doesn't converge to correct value

问题

使用fsolve时，迭代的值如下所示：

  1. 1
  2. 1.0
  3. 1.0
  4. 1.0000000149011612
  5. 101.0
  6. nan
  7. 1.0000000149011612
  8. 101.0
  9. nan
  10. nan
  11. nan
  12. nan
  13. nan
  14. nan
  15. nan

代码迭代的模式似乎不受给定的P、T、kf或kp值的影响。此外，代码还生成了运行时警告。因此，您的问题是代码本身是否存在问题。

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

希望这有所帮助。如果您需要进一步的解释或帮助，请随时提出。

英文:

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
  2. 1.0
  3. 1.0
  4. 1.0000000149011612
  5. 101.0
  6. nan
  7. 1.0000000149011612
  8. 101.0
  9. nan
  10. nan
  11. nan
  12. nan
  13. nan
  14. nan
  15. 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?

  1. # -*- coding: utf-8 -*-
  2. """
  3. Created on Mon Apr  3 16:46:15 2023
  5. @author: houle
  6. """
  8. import numpy as np
  9. import math
  10. import matplotlib.pyplot as plt
  11. from scipy.integrate import solve_ivp
  12. from scipy.optimize import fsolve
  14. class parametre_antoine():
  15.     A = 8.13484
  16.     B = 1662.48
  17.     C = 238.131
  18.     mmhg_atm = 760
  19.     
  20. prm_antoine = parametre_antoine() 
  22. rho = 0.8 #kg/L
  23. Tin = 110 #C   
  24. R = 1
  25. kp = -200
  26. kf = 300
  27. M = 46.068/1000 #kg/mol
  28. L = 5
  29. M = 46.068 #g/mol
  30. Vtot = np.pi*R**2*L
  32. theta = [rho,R,kp,kf,M,L,prm_antoine,Vtot]
  34. Fin = 20000
  35. u = [Fin]
  37. P = 3.3
  38. T = 300
  40. def Aire(R,L,h):
  41.     return 2*L*np.sqrt(R**2-(R-h)**2)
  43. def dv_dh(R,L,h):
  44.     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))
  45.     return dVdh
  47. def dh_dt(h,theta,u,T,P):
  48.     [rho,R,kp,kf,M,L,prm_antoine,Vtot] = theta
  49.     [Fin] = u
  50.     dhdt =(Fin- kp*(psat(prm_antoine,T)-P)*Aire(R,L,h)-kf*h**0.5)/rho/ dv_dh(R,L,h)
  51.     print(h)
  52.     return dhdt
  54. def psat(prm,T):
  55.     #T en Celcius
  56.     #Valide entre -114 et 243C
  57.     #Retourne la tension de vapeur en atm
  58.     p_mmhg = 10**(prm.A-prm.B/(T+prm.C))
  59.     p_atm = p_mmhg / prm.mmhg_atm
  60.     return p_atm
  62. x0 = [1]
  63. u0 = [Fin]
  64. x0 = fsolve(dh_dt, x0, args=(theta, u0,T,P))
  1. RuntimeWarning: invalid value encountered in sqrt
  2.   return 2*L*np.sqrt(R**2-(R-h)**2)

答案1

得分: 0

根据您提供的内容，以下是已经翻译好的代码部分：

  1. IIUC，有很多需要改变的事情。此代码运行并生成输出：
  3.     import numpy as np
  4.     from scipy.optimize import fsolve
  5.     
  6.     # 定义参数
  7.     rho = 0.8  # kg/L
  8.     R = 1  # m
  9.     L = 5  # m
  10.     kp = -200
  11.     kf = 300
  12.     M = 46.068  # g/mol
  13.     prm_A = 8.13484
  14.     prm_B = 1662.48
  15.     prm_C = 238.131
  16.     mmhg_atm = 760
  17.     Vtot = np.pi * R ** 2 * L
  18.     P = 3.3  # atm
  19.     T = 300  # K
  20.     Fin = 20000  # L/h
  21.     
  22.     # 定义函数
  23.     def psat(T):
  24.         p_mmhg = 10 ** (prm_A - prm_B / (T + prm_C))
  25.         p_atm = p_mmhg / mmhg_atm
  26.         return p_atm
  27.     
  28.     def Aire(R, L, h):
  29.         if h >= R:
  30.             return 0
  31.         else:
  32.             return 2 * L * np.sqrt(R ** 2 - (R - h) ** 2)
  33.     
  34.     def dv_dh(R, L, h):
  35.         if h >= R:
  36.             return 0
  37.         else:
  38.             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))
  39.     
  40.     def dh_dt(h):
  41.         return (Fin - kp * (psat(T) - P) * Aire(R, L, h) - kf * h ** 0.5) / rho / dv_dh(R, L, h)
  42.     
  43.     # 解决初始值 h 的问题
  44.     x0 = [0.5 * R]
  45.     h0 = fsolve(dh_dt, x0)[0]
  46.     
  47.     print(f"The initial value of h is {h0:.4f} m.")

请注意，已经翻译的代码部分中已将 HTML 实体 > 替换为正确的 Python 表达式 >=，以确保代码的正确性。

英文:

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

  1. import numpy as np
  2. from scipy.optimize import fsolve
  3. # Define the parameters
  4. rho = 0.8 # kg/L
  5. R = 1 # m
  6. L = 5 # m
  7. kp = -200
  8. kf = 300
  9. M = 46.068 # g/mol
  10. prm_A = 8.13484
  11. prm_B = 1662.48
  12. prm_C = 238.131
  13. mmhg_atm = 760
  14. Vtot = np.pi*R**2*L
  15. P = 3.3 # atm
  16. T = 300 # K
  17. Fin = 20000 # L/h
  18. # Define the functions
  19. def psat(T):
  20. p_mmhg = 10**(prm_A-prm_B/(T+prm_C))
  21. p_atm = p_mmhg / mmhg_atm
  22. return p_atm
  23. def Aire(R,L,h):
  24. if h >= R:
  25. return 0
  26. else:
  27. return 2*L*np.sqrt(R**2-(R-h)**2)
  28. def dv_dh(R,L,h):
  29. if h >= R:
  30. return 0
  31. else:
  32. 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))
  33. def dh_dt(h):
  34. return (Fin- kp*(psat(T)-P)*Aire(R,L,h)-kf*h**0.5)/rho/ dv_dh(R,L,h)
  35. # Solve for the initial value of h
  36. x0 = [0.5*R]
  37. h0 = fsolve(dh_dt, x0)[0]
  38. 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.

huangapple
  • 本文由 发表于 2023年4月4日 09:38:39
  • 转载请务必保留本文链接：https://go.coder-hub.com/75924877.html
  • fsolve
  • python
  • scipy
匿名

发表评论

匿名网友

:?: :razz: :sad: :evil: :!: :smile: :oops: :grin: :eek: :shock: :???: :cool: :lol: :mad: :twisted: :roll: :wink: :idea: :arrow: :neutral: :cry: :mrgreen:

确定