如何解决以下使用scipy的代码中的零除错误？

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

How to resolve zero division error for the following code which uses scipy?

问题

这是我的Python代码：

  1. import numpy as np
  2. from scipy.integrate import quad
  4. a = 0.15
  6. def A(z):
  7.     return -a * z**2
  9. def integrand(x, z, B):
  10.     return np.sqrt(-(2/x)*(3*x*A(x).derivative(n=2) - 3*x*(A(x).derivative(n=1))**2 + 6*A(x).derivative(n=1) + 2*B**4*x**3 + 2*B**2*x))
  12. def Phi(z, B):
  13.     result, _ = quad(integrand, 0, z, args=(z, B))
  14.     return result
  16. phi_values = [Phi(z, 0) for z in range(101)]

我遇到了"ZeroDivisionError: float division by zero"的错误。如何解决这个问题？我尝试使用sympy，但没有成功。

英文:

Here's my python code:

  1. import numpy as np
  2. from scipy.integrate import quad
  4. a = 0.15
  6. def A(z):
  7.     return -a * z**2
  9. def integrand(x, z, B):
  10.     return np.sqrt(-(2/x)*(3*x*A(x).derivative(n=2) - 3*x*(A(x).derivative(n=1))**2 + 6*A(x).derivative(n=1) + 2*B**4*x**3 + 2*B**2*x))
  12. def Phi(z, B):
  13.     result, _ = quad(integrand, 0, z, args=(z, B))
  14.     return result
  16. phi_values = [Phi(z, 0) for z in range(101)]

I am getting a "ZeroDivisionError: float division by zero"

How to resolve this issue?

I tried using sympy but to no avail.

答案1

得分: 1

在修复了其他可能存在的问题并进行了一些猜测后，同时将 epsilon 添加为积分下界，

  1. import numpy as np
  2. from scipy.integrate import quad
  4. a = 0.15
  6. # def A(z: float) -> float:
  7. #     return -a * z**2
  8. A = np.poly1d((-a, 0, 0))
  9. Ad1 = A.deriv(m=1)
  10. Ad2 = A.deriv(m=2)
  12. def integrand(x: float, z: float, B: float) -> float:
  13.     return np.sqrt(
  14.         -(2/x)*(
  15.             3*x*Ad2(x)
  16.             - 3*x*Ad1(x)**2
  17.             + 6*Ad1(x)
  18.             + 2*B**4*x**3
  19.             + 2*B**2*x
  20.         )
  21.     )
  23. def Phi(z: float, B: float):
  24.     epsilon = 1e-6
  25.     result, _ = quad(func=integrand, a=epsilon, b=z, args=(z, B))
  26.     return result
  28. phi_values = [Phi(z, 0) for z in range(101)]
  29. print(np.array(phi_values))
  1. [-2.32379001e-06  2.36195636e+00  4.94106004e+00  7.90800675e+00
  2.   1.13771987e+01  1.54207025e+01  2.00839857e+01  2.53965776e+01
  3. ...
  4.   3.40309939e+03  3.47405018e+03  3.54573543e+03  3.61815514e+03
  5.   3.69130933e+03]

但更明智的做法是在积分之前对 x 进行解析分离：

  1. import numpy as np
  2. from numpy.polynomial import Polynomial
  3. from scipy.integrate import quad
  5. a = 0.15
  7. # def A(z: float) -> float:
  8. #     return -a * z**2
  9. A = Polynomial((0, 0, -a))
  10. Ad1 = A.deriv(m=1)
  11. Ad2 = A.deriv(m=2)
  13. # Divide by x prior to integrand
  14. Ad1ox = Ad1 // Polynomial((0, 1))
  16. def integrand(x: float, z: float, B: float) -> float:
  17.     return np.sqrt(
  18.         -2*(
  19.             3*Ad2(x)
  20.             - 3*Ad1(x)**2
  21.             + 6*Ad1ox(x)  # 6/x*Ad1(x)
  22.             + 2*B**4*x**2
  23.             + 2*B**2
  24.         )
  25.     )
  27. def Phi(z: float, B: float):
  28.     result, _ = quad(func=integrand, a=0, b=z, args=(z, B))
  29.     return result
  31. phi_values = [Phi(z, 0) for z in range(101)]
  32. print(np.array(phi_values))

这将产生相同的结果。

英文:

Fixing your other breakage re. derivatives with some wild guessing, and adding an epsilon as your lower integration bound,

  1. import numpy as np
  2. from scipy.integrate import quad
  4. a = 0.15
  7. # def A(z: float) -> float:
  8. #     return -a * z**2
  9. A = np.poly1d((-a, 0, 0))
  10. Ad1 = A.deriv(m=1)
  11. Ad2 = A.deriv(m=2)
  14. def integrand(x: float, z: float, B: float) -> float:
  15.     return np.sqrt(
  16.         -(2/x)*(
  17.             3*x*Ad2(x)
  18.             - 3*x*Ad1(x)**2
  19.             + 6*Ad1(x)
  20.             + 2*B**4*x**3
  21.             + 2*B**2*x
  22.         )
  23.     )
  26. def Phi(z: float, B: float):
  27.     epsilon = 1e-6
  28.     result, _ = quad(func=integrand, a=epsilon, b=z, args=(z, B))
  29.     return result
  32. phi_values = [Phi(z, 0) for z in range(101)]
  33. print(np.array(phi_values))
  1. [-2.32379001e-06  2.36195636e+00  4.94106004e+00  7.90800675e+00
  2.   1.13771987e+01  1.54207025e+01  2.00839857e+01  2.53965776e+01
  3. ...
  4.   3.40309939e+03  3.47405018e+03  3.54573543e+03  3.61815514e+03
  5.   3.69130933e+03]

But the saner thing to do is analytically divide out x:

  1. import numpy as np
  2. from numpy.polynomial import Polynomial
  3. from scipy.integrate import quad
  5. a = 0.15
  8. # def A(z: float) -> float:
  9. #     return -a * z**2
  10. A = Polynomial((0, 0, -a))
  11. Ad1 = A.deriv(m=1)
  12. Ad2 = A.deriv(m=2)
  14. # Divide by x prior to integrand
  15. Ad1ox = Ad1 // Polynomial((0, 1))
  18. def integrand(x: float, z: float, B: float) -> float:
  19.     return np.sqrt(
  20.         -2*(
  21.             3*Ad2(x)
  22.             - 3*Ad1(x)**2
  23.             + 6*Ad1ox(x)  # 6/x*Ad1(x)
  24.             + 2*B**4*x**2
  25.             + 2*B**2
  26.         )
  27.     )
  30. def Phi(z: float, B: float):
  31.     result, _ = quad(func=integrand, a=0, b=z, args=(z, B))
  32.     return result
  35. phi_values = [Phi(z, 0) for z in range(101)]
  36. print(np.array(phi_values))

which yields the same results.

huangapple
  • 本文由 发表于 2023年3月31日 22:08:36
  • 转载请务必保留本文链接：https://go.coder-hub.com/75899485.html
  • divide-by-zero
  • error-handling
  • numpy
  • runtime-error
  • scipy
匿名

发表评论

匿名网友

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

确定