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

Question

Here's my python code:

import numpy as np
from scipy.integrate import quad

a = 0.15

def A(z):
    return -a * z**2

def integrand(x, z, B):
    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))

def Phi(z, B):
    result, _ = quad(integrand, 0, z, args=(z, B))
    return result

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.

Answer 1

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

import numpy as np
from scipy.integrate import quad

a = 0.15


# def A(z: float) -> float:
#     return -a * z**2
A = np.poly1d((-a, 0, 0))
Ad1 = A.deriv(m=1)
Ad2 = A.deriv(m=2)


def integrand(x: float, z: float, B: float) -> float:
    return np.sqrt(
        -(2/x)*(
            3*x*Ad2(x)
            - 3*x*Ad1(x)**2
            + 6*Ad1(x)
            + 2*B**4*x**3
            + 2*B**2*x
        )
    )


def Phi(z: float, B: float):
    epsilon = 1e-6
    result, _ = quad(func=integrand, a=epsilon, b=z, args=(z, B))
    return result


phi_values = [Phi(z, 0) for z in range(101)]
print(np.array(phi_values))

[-2.32379001e-06  2.36195636e+00  4.94106004e+00  7.90800675e+00
  1.13771987e+01  1.54207025e+01  2.00839857e+01  2.53965776e+01
...
  3.40309939e+03  3.47405018e+03  3.54573543e+03  3.61815514e+03
  3.69130933e+03]

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

import numpy as np
from numpy.polynomial import Polynomial
from scipy.integrate import quad

a = 0.15


# def A(z: float) -> float:
#     return -a * z**2
A = Polynomial((0, 0, -a))
Ad1 = A.deriv(m=1)
Ad2 = A.deriv(m=2)

# Divide by x prior to integrand
Ad1ox = Ad1 // Polynomial((0, 1))


def integrand(x: float, z: float, B: float) -> float:
    return np.sqrt(
        -2*(
            3*Ad2(x)
            - 3*Ad1(x)**2
            + 6*Ad1ox(x)  # 6/x*Ad1(x)
            + 2*B**4*x**2
            + 2*B**2
        )
    )


def Phi(z: float, B: float):
    result, _ = quad(func=integrand, a=0, b=z, args=(z, B))
    return result


phi_values = [Phi(z, 0) for z in range(101)]
print(np.array(phi_values))

which yields the same results.