How to use sympy - 10 common examples

def print_error_bound(name, indent, variable_name, series_length, coefficients, nextc):
    is_alternating = sign(coefficients[-1]) != sign(coefficients[-2])
    if is_alternating:
        nextc = abs(nextc)
        print(f"{indent}// Error is at most {format_real_coefficient(nextc)}*{variable_name}**{series_length} when {variable_name} >= 0")
        ulp1 = 2.220446049250313e-16
        if name == "2: Digamma at 1":
            offset = S(1e6)
        elif name == "3: Digamma at 2":
            offset = 1 + digamma(1)
        elif name == "4: Digamma asymptotic":
            offset = S(1)
        elif name == "5: Trigamma at 1":
            offset = S(1e8)
        elif name == "7: Tetragamma at 1":
            offset = S(2e12)
        elif name == "6: Trigamma asymptotic":
            offset = S(1)
        elif name == "8: Tetragamma asymptotic":
            offset = S(12**-3)
        elif name == "15: log(exp(x) - 1) / x":
            offset = S(-log(1e-3))
        else:
            offset = abs(coefficients[0])
        if offset == 0:
            offset = abs(coefficients[1])
            if offset == 0:
                offset = abs(coefficients[2])

v = {}
        v[(0, 0)] = [self._print(a) for a in e.an]
        v[(0, 1)] = [self._print(a) for a in e.aother]
        v[(1, 0)] = [self._print(b) for b in e.bm]
        v[(1, 1)] = [self._print(b) for b in e.bother]

        P = self._print(e.argument)
        P.baseline = P.height()//2

        vp = {}
        for idx in v:
            vp[idx] = self._hprint_vec(v[idx])

        for i in range(2):
            maxw = max(vp[(0, i)].width(), vp[(1, i)].width())
            for j in range(2):
                s = vp[(j, i)]
                left = (maxw - s.width()) // 2
                right = maxw - left - s.width()
                s = prettyForm(*s.left(' ' * left))
                s = prettyForm(*s.right(' ' * right))
                vp[(j, i)] = s

        D1 = prettyForm(*vp[(0, 0)].right('  ', vp[(0, 1)]))
        D1 = prettyForm(*D1.below(' '))
        D2 = prettyForm(*vp[(1, 0)].right('  ', vp[(1, 1)]))
        D = prettyForm(*D1.below(D2))

        # make sure that the argument `z' is centred vertically
        D.baseline = D.height()//2

        # insert horizontal separator

if n&(n - 1): # not a power of 2
        b += 1
        n = 2**b

    a += [S.Zero]*(n - len(a))
    for i in range(1, n):
        j = int(ibin(i, b, str=True)[::-1], 2)
        if i < j:
            a[i], a[j] = a[j], a[i]

    ang = -2*pi/n if inverse else 2*pi/n

    if dps is not None:
        ang = ang.evalf(dps + 2)

    w = [cos(ang*i) + I*sin(ang*i) for i in range(n // 2)]

    h = 2
    while h <= n:
        hf, ut = h // 2, n // h
        for i in range(0, n, h):
            for j in range(hf):
                u, v = a[i + j], expand_mul(a[i + j + hf]*w[ut * j])
                a[i + j], a[i + j + hf] = u + v, u - v
        h *= 2

    if inverse:
        a = [(x/n).evalf(dps) for x in a] if dps is not None \
                            else [x/n for x in a]

    return a

def xiu(n):
    points = []
    for k in range(n + 1):
        pt = []
        # Slight adaptation:
        # The article has points for the weight 1/sqrt(2*pi) exp(−x**2/2)
        # so divide by sqrt(2) to adapt for 1/sqrt(pi) exp(−x ** 2)
        for r in range(1, n // 2 + 1):
            alpha = (2 * r * k * pi) / (n + 1)
            pt += [cos(alpha), sin(alpha)]
        if n % 2 == 1:
            pt += [(-1) ** k / sqrt(2)]
        points.append(pt)

    points = numpy.array(points)
    weights = numpy.full(n + 1, frac(1, n + 1))
    return Enr2Scheme("Xiu", n, weights, points, 2, source)

quotient between successive terms must be a quotient of integer
    polynomials.
    """
    from sympy import Float, hypersimp, lambdify

    if prec == float('inf'):
        raise NotImplementedError('does not support inf prec')

    if start:
        expr = expr.subs(n, n + start)
    hs = hypersimp(expr, n)
    if hs is None:
        raise NotImplementedError("a hypergeometric series is required")
    num, den = hs.as_numer_denom()

    func1 = lambdify(n, num)
    func2 = lambdify(n, den)

    h, g, p = check_convergence(num, den, n)

    if h < 0:
        raise ValueError("Sum diverges like (n!)^%i" % (-h))

    term = expr.subs(n, 0)
    if not term.is_Rational:
        raise NotImplementedError("Non rational term functionality is not implemented.")

    # Direct summation if geometric or faster
    if h > 0 or (h == 0 and abs(g) > 1):
        term = (MPZ(term.p) << prec) // term.q
        s = term
        k = 1

def eval(cls, arg):
        if arg.is_Number:
            if arg is S.NaN:
                return S.NaN
            elif arg is S.Infinity:
                return S.Infinity
            elif intlike(arg):
                if arg.is_positive:
                    return factorial(arg - 1)
                else:
                    return S.ComplexInfinity
            elif arg.is_Rational:
                if arg.q == 2:
                    n = abs(arg.p) // arg.q

                    if arg.is_positive:
                        k, coeff = n, S.One
                    else:
                        n = k = n + 1

def single_beam_trapezload_test(curr_doc: CurrentDoc):
    start_knot = Knot.Knot(0, 0, 0, ElSupEnum.SUPPORT_FIXED_END.value, 0)
    end_knot = Knot.Knot(1, 1, 0, ElSupEnum.SUPPORT_ROLLER_END.value,  0)
    q = Symbol('p')
    x = Symbol('x')
    l = Symbol('l')
    lineload = [0, q * x]
    temp_prop = TempProps.TempProps(0, 0, 0)
    knot_list = [[start_knot, end_knot]]
    elementlist = []
    ele = ElementCalculation(0, start_knot, end_knot, l, lineload, temp_prop)
    start_knot.add_coupled_el(0)
    end_knot.add_coupled_el(0)
    elementlist.append(ele)
    functions, x, l_list = CalculationElement(elementlist)
    testbox.print_graphs(functions, x, l_list, knot_list)

testbench.test()

    if cli:
        r = ahkab.run(mycircuit, [symbolic_sim, ac_sim])
        E = r['symbolic'][0].as_symbol('E1')
        out_hp = sympy.limit(r['symbolic'][0]['VU1o'], E, sympy.oo, '+')
        out_bp = sympy.limit(r['symbolic'][0]['VU2o'], E, sympy.oo, '+')
        out_lp = sympy.limit(r['symbolic'][0]['VU3o'], E, sympy.oo, '+')
        out_hp = out_hp.simplify()
        out_bp = out_bp.simplify()
        out_lp = out_lp.simplify()
        print("VU1o =", out_hp)
        print("VU2o =", out_bp)
        print("VU3o =", out_lp)

        w = sympy.Symbol('w')
        out_hp = out_hp.subs({r['symbolic'][0].as_symbol('RF1'):10e3,
                              r['symbolic'][0].as_symbol('C10'):15e-9,
                              r['symbolic'][0].as_symbol('V1'):1,
                              r['symbolic'][0].as_symbol('s'):1j*w,
                              })
        out_bp = out_bp.subs({r['symbolic'][0].as_symbol('RF1'):10e3,
                              r['symbolic'][0].as_symbol('C10'):15e-9,
                              r['symbolic'][0].as_symbol('V1'):1,
                              r['symbolic'][0].as_symbol('s'):1j*w,
                              })
        out_lp = out_lp.subs({r['symbolic'][0].as_symbol('RF1'):10e3,
                              r['symbolic'][0].as_symbol('C10'):15e-9,
                              r['symbolic'][0].as_symbol('V1'):1,
                              r['symbolic'][0].as_symbol('s'):1j*w,
                              })
        out_lp = sympy.lambdify((w,), out_lp, modules='numpy')

def test_hyper():
    for x in sorted(exparg):
        test("erf", x, N(sp.erf(x)))
    for x in sorted(exparg):
        test("erfc", x, N(sp.erfc(x)))

    gamarg = FiniteSet(*(x+S(1)/12 for x in exparg))
    betarg = ProductSet(gamarg, gamarg)
    for x in sorted(gamarg):
        test("lgamma", x, N(sp.log(abs(sp.gamma(x)))))
    for x in sorted(gamarg):
        test("gamma", x, N(sp.gamma(x)))
    for x, y in sorted(betarg, key=lambda (x, y): (y, x)):
        test("beta", x, y, N(sp.beta(x, y)))

    pgamarg = FiniteSet(S(1)/12, S(1)/3, S(3)/2, 5)
    pgamargp = ProductSet(gamarg & Interval(0, oo, True), pgamarg)
    for a, x in sorted(pgamargp):
        test("pgamma", a, x, N(sp.lowergamma(a, x)))
    for a, x in sorted(pgamargp):
        test("pgammac", a, x, N(sp.uppergamma(a, x)))
    for a, x in sorted(pgamargp):
        test("pgammar", a, x, N(sp.lowergamma(a, x)/sp.gamma(a)))
    for a, x in sorted(pgamargp):
        test("pgammarc", a, x, N(sp.uppergamma(a, x)/sp.gamma(a)))
    for a, x in sorted(pgamargp):
        test("ipgammarc", a, N(sp.uppergamma(a, x)/sp.gamma(a)), x)

    pbetargp = [(a, b, x) for a, b, x in ProductSet(betarg, pgamarg)
                if a > 0 and b > 0 and x < 1]
    pbetargp.sort(key=lambda (a, b, x): (b, a, x))
    for a, b, x in pbetargp:

def substitute(self, expression: sympy.Expr, substitutions: dict):
        for key, value in substitutions.items():
            if not isinstance(value, sympy.Expr):
                substitutions[key] = sympy.sympify(value)
        return expression.subs(substitutions, simultaneous=True).doit()