How to use the sympy.Rational function in sympy

def _represent_base(self, basis, **options):
        j = sympify(options.get('j', Rational(1, 2)))
        # TODO: move evaluation up to represent function/implement elsewhere
        evaluate = sympify(options.get('doit'))
        size, mvals = m_values(j)
        result = zeros(size, size)
        for p in range(size):
            for q in range(size):
                me = self.matrix_element(j, mvals[p], j, mvals[q])
                if evaluate:
                    result[p, q] = me.doit()
                else:
                    result[p, q] = me
        return result

z2 = float(2*c**2 + c - 2)
                z3 = float(12*c**4 - 18*c**3 + 18*c**2 - 11*c + 2)
                z4 = float(36*c**4 - 36*c**3 + 13*c**2 - 8*c + 4)
                z5 = float(69*c**3 - 62*c**2 + 28*c - 8)
                z6 = float(34*c**4 - 46*c**3 + 34*c**2 - 13*c + 2)
                B1 = 0.924574
                B2 = (12*c*(c-1)*(3*z2-z1) - (3*z2-z1)**2)/(144*c*(3*c-2)*(c-1)**2)
                B3 = (-24*(3*c-2)*(c-1)**2)/((3*z2-z1)**2 - 12*c*(c-1)*(3*z2-z1))
                A1 = 0.0
                A2 = (-z1*(6*c**2 - 4*c + 1) + 3*z3)/((2*c+1)*z1 - 3*(c+2)*(2*c-1)**2)
                A3 = (-z1*z4 + 108*(2*c-1)*c**5 - 3*(2*c-1)*z5)/(24*z1*c*(c-1)**4 + 72*c*z6 + 72*c**6 * (2*c-13))
                cls.solution_coeffs.value = [B1, B2, B3]
                cls.stage_coeffs.value = [A1, A2, A3]
            else:
                A1, A2, A3 = 0, Rational(-5, 9), Rational(-153, 128)
                B1, B2, B3 = Rational(1, 3), Rational(15, 16), Rational(8, 15)
                cls.solution_coeffs.value = [B1, B2, B3]
                cls.stage_coeffs.value = [A1, A2, A3]
        elif cls.order == 4:
            A1, A2, A3, A4, A5 = 0, -0.4178904745, -1.192151694643, -1.697784692471, -1.514183444257
            B1, B2, B3, B4, B5 = 0.1496590219993, 0.3792103129999, 0.8229550293869, 0.6994504559488, 0.1530572479681
            cls.solution_coeffs.value = [B1, B2, B3, B4, B5]
            cls.stage_coeffs.value = [A1, A2, A3, A4, A5]
        else:
            raise NotImplementedError("Only 3rd and 4th order RK schemes are currently implemented.")
        return

def stroud_enr_5_4(n):
    r = sqrt(((n + 2) * (n + 3) + (n - 1) * (n + 3) * sqrt(2 * (n + 2))) / n)
    s = sqrt(((n + 2) * (n + 3) - (n + 3) * sqrt(2 * (n + 2))) / n)
    A = frac(4 * n + 6, (n + 2) * (n + 3))
    B = frac(n + 1, (n + 2) * (n + 3) * 2 ** n)
    data = [(A, numpy.full((1, n), 0)), (B, fsd(n, (r, 1), (s, n - 1)))]

    points, weights = untangle(data)
    return EnrScheme("Stroud Enr 5-4", n, weights, points, 5, source, 1.684e-11)

def mclaren_05():
    degree = 9

    r, s = [sqrt((5 + pm_ * sqrt(5)) / 10) for pm_ in [+1, -1]]
    u, v = [sqrt((3 - pm_ * sqrt(5)) / 6) for pm_ in [+1, -1]]

    B1 = frac(25, 840)
    B2 = frac(27, 840)

    data = {"rs0": [[B1, r, s], [B2, u, v]], "a3": [B2]}
    points, weights = expand_symmetries(data)
    theta_phi = cartesian_to_spherical_sympy(points)
    return U3Scheme("McLaren 5", weights, points, theta_phi, degree, source)

def dobrodeev_1970(n):
    assert n >= 3, f"Only works for n >= 3, not n = {n}"

    A = frac(1, 8)
    B = frac(19 - 5 * n, 20)
    alpha = 35 * n * (5 * n - 33)
    C = frac((alpha + 2114) ** 3, 700 * (alpha + 1790) * (alpha + 2600))
    D = frac(729, 1750) * frac(alpha + 2114, alpha + 2600)
    E = frac(n * (n - 1) * (n - frac(47, 10)), 3) - 2 * n * (C + D) + frac(729, 125)

    a = sqrt(frac(3, 5))
    b = a
    c = sqrt(frac(3, 5) * frac(alpha + 1790, alpha + 2114))
    data = [
        (A, fsd(n, (a, 3))),
        (B, fsd(n, (b, 2))),
        (C, fsd(n, (c, 1))),
        (D, fsd(n, (1.0, 1))),
        (E, z(n)),
    ]

    points, weights = untangle(data)
    weights *= frac(125, 729)
    return CnScheme("Dobrodeev 1970", n, weights, points, 7, _source, 8.100e-13)

def Ixx(sections):
    r"""
    from http://en.wikipedia.org/wiki/Second_moment_of_area
    \f[ J_{xx} = \frac{1}{12} \sum_{i = 1}^{n-1} ( y_i^2 + y_i y_{i+1} + y_{i+1}^2 ) a_i \f]
    \f[ J_{yy} = \frac{1}{12} \sum_{i = 1}^{n-1} ( x_i^2 + x_i x_{i+1} + x_{i+1}^2 ) a_i \f]
    \f[ J_{xy} = \frac{1}{24} \sum_{i = 1}^{n-1} ( x_i y_{i+1} + 2 x_i y_i + 2 x_{i+1} y_{i+1} + x_{i+1} y_i ) a_i \f]
    """
    h = Symbol('h')
    b = Symbol('b')
    Ixx = 0
    Iyy = 0
    Ixy = 0
    half = Rational(1, 2)
    CG = [0, 0]
    Area = 0
    o3 = Rational(1, 3)
    for i, section in enumerate(sections):
        (sign, p1, p2) = section
        #print("p%s = %s" %(i,p1))
        #print("p%s = %s" %(i+1,p2))
        #Atri = half * abs(p1[0]*p2[1])
        AA = Matrix([p1, p2])
        #print("AA = \n",AA)
        detAA = AA.det()
        msgAA = '%s' % (detAA)
        #print("msgAA = ",msgAA)
        if '-' in msgAA:
            detAA *= Rational(-1, 1)
        Atri = detAA
        #Atri = half*b*h
        #Atri = (p1[0]*p2[1] - p2[0]*p1[1])*half
        ixxi = p1[1] * p1[1] + p1[1] * p2[1] + p2[1] * p2[1]

# characteristic polynomial function
        # This is always fast, so no need for alternative
        # formulas
        p1 = np.poly(beta[:-1,:].astype(float)-np.tile(beta[-1,:].astype(float),(s,1)))
        q1 = np.poly(beta[:-1,:].astype(float))
        p = np.poly1d(p1[::-1])    # Numerator
        q = np.poly1d(q1[::-1])    # Denominator

    else: # Compute symbolically
        import sympy
        z = sympy.var('z')
        
        if explicit:
            v = 1 - alpha[:,1:].sum(1)
            alpha[:,0]=0.
            q1 = [sympy.Rational(1)]
        else:
            v = 1 - alpha.sum(1)

        alpha_star=sympy.Matrix(alpha[:-1,:])
        beta_star=sympy.Matrix(beta[:-1,:])
        I = sympy.eye(s)

        v_mp1 = v[-1]
        vstar = sympy.Matrix(v[:-1])
        alpha_mp1 = sympy.Matrix(alpha[-1,:]).T
        beta_mp1 = sympy.Matrix(beta[-1,:]).T

        if formula == 'det':
            xsym = I - alpha_star - z*beta_star + vstar/v_mp1 * (alpha_mp1+z*beta_mp1)
            p1 = sympy.simplify(xsym.det(method='berkowitz')*v_mp1)
            p1 = p1.as_poly(z).all_coeffs()

def schmid_4():
    points = numpy.array(
        [
            [0, (sqrt(3) + sqrt(15)) / 6],
            [0, (sqrt(3) - sqrt(15)) / 6],
            [+sqrt(15) / 5, (+sqrt(87) - 2 * sqrt(3)) / 15],
            [-sqrt(15) / 5, (+sqrt(87) - 2 * sqrt(3)) / 15],
            [+sqrt(15) / 5, (-sqrt(87) - 2 * sqrt(3)) / 15],
            [-sqrt(15) / 5, (-sqrt(87) - 2 * sqrt(3)) / 15],
        ]
    )
    weights = numpy.array(
        [
            frac(2, 9) - 2 * sqrt(5) / 45,
            frac(2, 9) + 2 * sqrt(5) / 45,
            frac(5, 36) + 5 * sqrt(29) / 18 / 29,
            frac(5, 36) + 5 * sqrt(29) / 18 / 29,
            frac(5, 36) - 5 * sqrt(29) / 18 / 29,
            frac(5, 36) - 5 * sqrt(29) / 18 / 29,
        ]
    )
    return C2Scheme("Schmid 4", weights, points, 4, source)

def liu_vinokur_14():
    weights = numpy.concatenate(
        [
            numpy.full(1, frac(16, 105)),
            numpy.full(4, frac(1, 280)),
            numpy.full(4, frac(81, 1400)),
            numpy.full(4, frac(64, 525)),
            numpy.full(6, frac(2, 105)),
        ]
    )
    points = numpy.concatenate(
        [
            _s4(symbolic=True),
            _r_alpha(1),
            _r_alpha(-frac(1, 3)),
            _r_alpha(frac(1, 2)),
            _r_beta(frac(1, 2)),
        ]
    )
    degree = 5

    return T3Scheme("Liu-Vinokur 14", weights, points, degree, source)

def __init__(self, part):
        # 2011 Q5a [5 lines] [2 marks]
        self.num_lines, self.num_marks = 5, 2

        self._qp = copy.deepcopy(part._qp)

        choices = [self._qp['domain'].left + sympy.Rational(i, 4) * self._qp['domain'].measure for i in [1, 2, 3]]
        self._qp['bound'] = random.choice(choices)

        self._qp['direction'] = random.choice(['left', 'right'])

        if self._qp['direction'] == 'left':
            self._qp['domain'] = sympy.Interval(self._qp['domain'].left, self._qp['bound'])
        else:
            self._qp['domain'] = sympy.Interval(self._qp['bound'], self._qp['domain'].right)