import sympy as smp
from scipy.misc import derivative
import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
from matplotlib import style

Solving equations

x, y, z, a, b, c, d, e, f = smp.symbols('x y z a b c d e f', real=True)

f1 = smp.Eq(a*x+b,c)

f1

$\displaystyle a x + b = c$

smp.solve(f1, x)

$\displaystyle \left[ \frac{- b + c}{a}\right]$

smp.init_printing()

smp.solve(f1, x)

$\displaystyle \left[ \frac{- b + c}{a}\right]$

smp.solve(f1, x, set = True)

$\displaystyle \left( \left[ x\right], \ \left\{\left( \frac{- b + c}{a},\right)\right\}\right)$

f2=smp.Eq(a*x**2+b*x+c,0)

f2

$\displaystyle a x^{2} + b x + c = 0$

smp.solve(f, x, set = True)

$\displaystyle \left( \left[ x\right], \ \left\{\left( \frac{- b + \sqrt{- 4 a c + b^{2}}}{2 a},\right), \left( - \frac{b + \sqrt{- 4 a c + b^{2}}}{2 a},\right)\right\}\right)$

f3=smp.Eq(a*x**3+b*x**2+c*x+d,0)
f3

$\displaystyle a x^{3} + b x^{2} + c x + d = 0$

smp.solve(f, x, set = True)

$\displaystyle \left( \left[ x\right], \ \left\{\left( - \frac{- \frac{3 c}{a} + \frac{b^{2}}{a^{2}}}{3 \sqrt[3]{\frac{\sqrt{- 4 \left(- \frac{3 c}{a} + \frac{b^{2}}{a^{2}}\right)^{3} + \left(\frac{27 d}{a} - \frac{9 b c}{a^{2}} + \frac{2 b^{3}}{a^{3}}\right)^{2}}}{2} + \frac{27 d}{2 a} - \frac{9 b c}{2 a^{2}} + \frac{b^{3}}{a^{3}}}} - \frac{\sqrt[3]{\frac{\sqrt{- 4 \left(- \frac{3 c}{a} + \frac{b^{2}}{a^{2}}\right)^{3} + \left(\frac{27 d}{a} - \frac{9 b c}{a^{2}} + \frac{2 b^{3}}{a^{3}}\right)^{2}}}{2} + \frac{27 d}{2 a} - \frac{9 b c}{2 a^{2}} + \frac{b^{3}}{a^{3}}}}{3} - \frac{b}{3 a},\right), \left( - \frac{- \frac{3 c}{a} + \frac{b^{2}}{a^{2}}}{3 \left(- \frac{1}{2} - \frac{\sqrt{3} i}{2}\right) \sqrt[3]{\frac{\sqrt{- 4 \left(- \frac{3 c}{a} + \frac{b^{2}}{a^{2}}\right)^{3} + \left(\frac{27 d}{a} - \frac{9 b c}{a^{2}} + \frac{2 b^{3}}{a^{3}}\right)^{2}}}{2} + \frac{27 d}{2 a} - \frac{9 b c}{2 a^{2}} + \frac{b^{3}}{a^{3}}}} - \frac{\left(- \frac{1}{2} - \frac{\sqrt{3} i}{2}\right) \sqrt[3]{\frac{\sqrt{- 4 \left(- \frac{3 c}{a} + \frac{b^{2}}{a^{2}}\right)^{3} + \left(\frac{27 d}{a} - \frac{9 b c}{a^{2}} + \frac{2 b^{3}}{a^{3}}\right)^{2}}}{2} + \frac{27 d}{2 a} - \frac{9 b c}{2 a^{2}} + \frac{b^{3}}{a^{3}}}}{3} - \frac{b}{3 a},\right), \left( - \frac{- \frac{3 c}{a} + \frac{b^{2}}{a^{2}}}{3 \left(- \frac{1}{2} + \frac{\sqrt{3} i}{2}\right) \sqrt[3]{\frac{\sqrt{- 4 \left(- \frac{3 c}{a} + \frac{b^{2}}{a^{2}}\right)^{3} + \left(\frac{27 d}{a} - \frac{9 b c}{a^{2}} + \frac{2 b^{3}}{a^{3}}\right)^{2}}}{2} + \frac{27 d}{2 a} - \frac{9 b c}{2 a^{2}} + \frac{b^{3}}{a^{3}}}} - \frac{\left(- \frac{1}{2} + \frac{\sqrt{3} i}{2}\right) \sqrt[3]{\frac{\sqrt{- 4 \left(- \frac{3 c}{a} + \frac{b^{2}}{a^{2}}\right)^{3} + \left(\frac{27 d}{a} - \frac{9 b c}{a^{2}} + \frac{2 b^{3}}{a^{3}}\right)^{2}}}{2} + \frac{27 d}{2 a} - \frac{9 b c}{2 a^{2}} + \frac{b^{3}}{a^{3}}}}{3} - \frac{b}{3 a},\right)\right\}\right)$

help(smp.solve)

Help on function solve in module sympy.solvers.solvers:

solve(f, *symbols, **flags)
    Algebraically solves equations and systems of equations.
    
    Explanation
    ===========
    
    Currently supported:
        - polynomial
        - transcendental
        - piecewise combinations of the above
        - systems of linear and polynomial equations
        - systems containing relational expressions
    
    Examples
    ========
    
    The output varies according to the input and can be seen by example:
    
        >>> from sympy import solve, Poly, Eq, Function, exp
        >>> from sympy.abc import x, y, z, a, b
        >>> f = Function('f')
    
    Boolean or univariate Relational:
    
        >>> solve(x < 3)
        (-oo < x) & (x < 3)
    
    
    To always get a list of solution mappings, use flag dict=True:
    
        >>> solve(x - 3, dict=True)
        [{x: 3}]
        >>> sol = solve([x - 3, y - 1], dict=True)
        >>> sol
        [{x: 3, y: 1}]
        >>> sol[0][x]
        3
        >>> sol[0][y]
        1
    
    
    To get a list of *symbols* and set of solution(s) use flag set=True:
    
        >>> solve([x**2 - 3, y - 1], set=True)
        ([x, y], {(-sqrt(3), 1), (sqrt(3), 1)})
    
    
    Single expression and single symbol that is in the expression:
    
        >>> solve(x - y, x)
        [y]
        >>> solve(x - 3, x)
        [3]
        >>> solve(Eq(x, 3), x)
        [3]
        >>> solve(Poly(x - 3), x)
        [3]
        >>> solve(x**2 - y**2, x, set=True)
        ([x], {(-y,), (y,)})
        >>> solve(x**4 - 1, x, set=True)
        ([x], {(-1,), (1,), (-I,), (I,)})
    
    Single expression with no symbol that is in the expression:
    
        >>> solve(3, x)
        []
        >>> solve(x - 3, y)
        []
    
    Single expression with no symbol given. In this case, all free *symbols*
    will be selected as potential *symbols* to solve for. If the equation is
    univariate then a list of solutions is returned; otherwise - as is the case
    when *symbols* are given as an iterable of length greater than 1 - a list of
    mappings will be returned:
    
        >>> solve(x - 3)
        [3]
        >>> solve(x**2 - y**2)
        [{x: -y}, {x: y}]
        >>> solve(z**2*x**2 - z**2*y**2)
        [{x: -y}, {x: y}, {z: 0}]
        >>> solve(z**2*x - z**2*y**2)
        [{x: y**2}, {z: 0}]
    
    When an object other than a Symbol is given as a symbol, it is
    isolated algebraically and an implicit solution may be obtained.
    This is mostly provided as a convenience to save you from replacing
    the object with a Symbol and solving for that Symbol. It will only
    work if the specified object can be replaced with a Symbol using the
    subs method:
    
    >>> solve(f(x) - x, f(x))
    [x]
    >>> solve(f(x).diff(x) - f(x) - x, f(x).diff(x))
    [x + f(x)]
    >>> solve(f(x).diff(x) - f(x) - x, f(x))
    [-x + Derivative(f(x), x)]
    >>> solve(x + exp(x)**2, exp(x), set=True)
    ([exp(x)], {(-sqrt(-x),), (sqrt(-x),)})
    
    >>> from sympy import Indexed, IndexedBase, Tuple, sqrt
    >>> A = IndexedBase('A')
    >>> eqs = Tuple(A[1] + A[2] - 3, A[1] - A[2] + 1)
    >>> solve(eqs, eqs.atoms(Indexed))
    {A[1]: 1, A[2]: 2}
    
        * To solve for a symbol implicitly, use implicit=True:
    
            >>> solve(x + exp(x), x)
            [-LambertW(1)]
            >>> solve(x + exp(x), x, implicit=True)
            [-exp(x)]
    
        * It is possible to solve for anything that can be targeted with
          subs:
    
            >>> solve(x + 2 + sqrt(3), x + 2)
            [-sqrt(3)]
            >>> solve((x + 2 + sqrt(3), x + 4 + y), y, x + 2)
            {y: -2 + sqrt(3), x + 2: -sqrt(3)}
    
        * Nothing heroic is done in this implicit solving so you may end up
          with a symbol still in the solution:
    
            >>> eqs = (x*y + 3*y + sqrt(3), x + 4 + y)
            >>> solve(eqs, y, x + 2)
            {y: -sqrt(3)/(x + 3), x + 2: -2*x/(x + 3) - 6/(x + 3) + sqrt(3)/(x + 3)}
            >>> solve(eqs, y*x, x)
            {x: -y - 4, x*y: -3*y - sqrt(3)}
    
        * If you attempt to solve for a number remember that the number
          you have obtained does not necessarily mean that the value is
          equivalent to the expression obtained:
    
            >>> solve(sqrt(2) - 1, 1)
            [sqrt(2)]
            >>> solve(x - y + 1, 1)  # /!\ -1 is targeted, too
            [x/(y - 1)]
            >>> [_.subs(z, -1) for _ in solve((x - y + 1).subs(-1, z), 1)]
            [-x + y]
    
        * To solve for a function within a derivative, use ``dsolve``.
    
    Single expression and more than one symbol:
    
        * When there is a linear solution:
    
            >>> solve(x - y**2, x, y)
            [(y**2, y)]
            >>> solve(x**2 - y, x, y)
            [(x, x**2)]
            >>> solve(x**2 - y, x, y, dict=True)
            [{y: x**2}]
    
        * When undetermined coefficients are identified:
    
            * That are linear:
    
                >>> solve((a + b)*x - b + 2, a, b)
                {a: -2, b: 2}
    
            * That are nonlinear:
    
                >>> solve((a + b)*x - b**2 + 2, a, b, set=True)
                ([a, b], {(-sqrt(2), sqrt(2)), (sqrt(2), -sqrt(2))})
    
        * If there is no linear solution, then the first successful
          attempt for a nonlinear solution will be returned:
    
            >>> solve(x**2 - y**2, x, y, dict=True)
            [{x: -y}, {x: y}]
            >>> solve(x**2 - y**2/exp(x), x, y, dict=True)
            [{x: 2*LambertW(-y/2)}, {x: 2*LambertW(y/2)}]
            >>> solve(x**2 - y**2/exp(x), y, x)
            [(-x*sqrt(exp(x)), x), (x*sqrt(exp(x)), x)]
    
    Iterable of one or more of the above:
    
        * Involving relationals or bools:
    
            >>> solve([x < 3, x - 2])
            Eq(x, 2)
            >>> solve([x > 3, x - 2])
            False
    
        * When the system is linear:
    
            * With a solution:
    
                >>> solve([x - 3], x)
                {x: 3}
                >>> solve((x + 5*y - 2, -3*x + 6*y - 15), x, y)
                {x: -3, y: 1}
                >>> solve((x + 5*y - 2, -3*x + 6*y - 15), x, y, z)
                {x: -3, y: 1}
                >>> solve((x + 5*y - 2, -3*x + 6*y - z), z, x, y)
                {x: 2 - 5*y, z: 21*y - 6}
    
            * Without a solution:
    
                >>> solve([x + 3, x - 3])
                []
    
        * When the system is not linear:
    
            >>> solve([x**2 + y -2, y**2 - 4], x, y, set=True)
            ([x, y], {(-2, -2), (0, 2), (2, -2)})
    
        * If no *symbols* are given, all free *symbols* will be selected and a
          list of mappings returned:
    
            >>> solve([x - 2, x**2 + y])
            [{x: 2, y: -4}]
            >>> solve([x - 2, x**2 + f(x)], {f(x), x})
            [{x: 2, f(x): -4}]
    
        * If any equation does not depend on the symbol(s) given, it will be
          eliminated from the equation set and an answer may be given
          implicitly in terms of variables that were not of interest:
    
            >>> solve([x - y, y - 3], x)
            {x: y}
    
    **Additional Examples**
    
    ``solve()`` with check=True (default) will run through the symbol tags to
    elimate unwanted solutions. If no assumptions are included, all possible
    solutions will be returned:
    
        >>> from sympy import Symbol, solve
        >>> x = Symbol("x")
        >>> solve(x**2 - 1)
        [-1, 1]
    
    By using the positive tag, only one solution will be returned:
    
        >>> pos = Symbol("pos", positive=True)
        >>> solve(pos**2 - 1)
        [1]
    
    Assumptions are not checked when ``solve()`` input involves
    relationals or bools.
    
    When the solutions are checked, those that make any denominator zero
    are automatically excluded. If you do not want to exclude such solutions,
    then use the check=False option:
    
        >>> from sympy import sin, limit
        >>> solve(sin(x)/x)  # 0 is excluded
        [pi]
    
    If check=False, then a solution to the numerator being zero is found: x = 0.
    In this case, this is a spurious solution since $\sin(x)/x$ has the well
    known limit (without dicontinuity) of 1 at x = 0:
    
        >>> solve(sin(x)/x, check=False)
        [0, pi]
    
    In the following case, however, the limit exists and is equal to the
    value of x = 0 that is excluded when check=True:
    
        >>> eq = x**2*(1/x - z**2/x)
        >>> solve(eq, x)
        []
        >>> solve(eq, x, check=False)
        [0]
        >>> limit(eq, x, 0, '-')
        0
        >>> limit(eq, x, 0, '+')
        0
    
    **Disabling High-Order Explicit Solutions**
    
    When solving polynomial expressions, you might not want explicit solutions
    (which can be quite long). If the expression is univariate, ``CRootOf``
    instances will be returned instead:
    
        >>> solve(x**3 - x + 1)
        [-1/((-1/2 - sqrt(3)*I/2)*(3*sqrt(69)/2 + 27/2)**(1/3)) - (-1/2 -
        sqrt(3)*I/2)*(3*sqrt(69)/2 + 27/2)**(1/3)/3, -(-1/2 +
        sqrt(3)*I/2)*(3*sqrt(69)/2 + 27/2)**(1/3)/3 - 1/((-1/2 +
        sqrt(3)*I/2)*(3*sqrt(69)/2 + 27/2)**(1/3)), -(3*sqrt(69)/2 +
        27/2)**(1/3)/3 - 1/(3*sqrt(69)/2 + 27/2)**(1/3)]
        >>> solve(x**3 - x + 1, cubics=False)
        [CRootOf(x**3 - x + 1, 0),
         CRootOf(x**3 - x + 1, 1),
         CRootOf(x**3 - x + 1, 2)]
    
    If the expression is multivariate, no solution might be returned:
    
        >>> solve(x**3 - x + a, x, cubics=False)
        []
    
    Sometimes solutions will be obtained even when a flag is False because the
    expression could be factored. In the following example, the equation can
    be factored as the product of a linear and a quadratic factor so explicit
    solutions (which did not require solving a cubic expression) are obtained:
    
        >>> eq = x**3 + 3*x**2 + x - 1
        >>> solve(eq, cubics=False)
        [-1, -1 + sqrt(2), -sqrt(2) - 1]
    
    **Solving Equations Involving Radicals**
    
    Because of SymPy's use of the principle root, some solutions
    to radical equations will be missed unless check=False:
    
        >>> from sympy import root
        >>> eq = root(x**3 - 3*x**2, 3) + 1 - x
        >>> solve(eq)
        []
        >>> solve(eq, check=False)
        [1/3]
    
    In the above example, there is only a single solution to the
    equation. Other expressions will yield spurious roots which
    must be checked manually; roots which give a negative argument
    to odd-powered radicals will also need special checking:
    
        >>> from sympy import real_root, S
        >>> eq = root(x, 3) - root(x, 5) + S(1)/7
        >>> solve(eq)  # this gives 2 solutions but misses a 3rd
        [CRootOf(7*x**5 - 7*x**3 + 1, 1)**15,
        CRootOf(7*x**5 - 7*x**3 + 1, 2)**15]
        >>> sol = solve(eq, check=False)
        >>> [abs(eq.subs(x,i).n(2)) for i in sol]
        [0.48, 0.e-110, 0.e-110, 0.052, 0.052]
    
    The first solution is negative so ``real_root`` must be used to see that it
    satisfies the expression:
    
        >>> abs(real_root(eq.subs(x, sol[0])).n(2))
        0.e-110
    
    If the roots of the equation are not real then more care will be
    necessary to find the roots, especially for higher order equations.
    Consider the following expression:
    
        >>> expr = root(x, 3) - root(x, 5)
    
    We will construct a known value for this expression at x = 3 by selecting
    the 1-th root for each radical:
    
        >>> expr1 = root(x, 3, 1) - root(x, 5, 1)
        >>> v = expr1.subs(x, -3)
    
    The ``solve`` function is unable to find any exact roots to this equation:
    
        >>> eq = Eq(expr, v); eq1 = Eq(expr1, v)
        >>> solve(eq, check=False), solve(eq1, check=False)
        ([], [])
    
    The function ``unrad``, however, can be used to get a form of the equation
    for which numerical roots can be found:
    
        >>> from sympy.solvers.solvers import unrad
        >>> from sympy import nroots
        >>> e, (p, cov) = unrad(eq)
        >>> pvals = nroots(e)
        >>> inversion = solve(cov, x)[0]
        >>> xvals = [inversion.subs(p, i) for i in pvals]
    
    Although ``eq`` or ``eq1`` could have been used to find ``xvals``, the
    solution can only be verified with ``expr1``:
    
        >>> z = expr - v
        >>> [xi.n(chop=1e-9) for xi in xvals if abs(z.subs(x, xi).n()) < 1e-9]
        []
        >>> z1 = expr1 - v
        >>> [xi.n(chop=1e-9) for xi in xvals if abs(z1.subs(x, xi).n()) < 1e-9]
        [-3.0]
    
    Parameters
    ==========
    
    f :
        - a single Expr or Poly that must be zero
        - an Equality
        - a Relational expression
        - a Boolean
        - iterable of one or more of the above
    
    symbols : (object(s) to solve for) specified as
        - none given (other non-numeric objects will be used)
        - single symbol
        - denested list of symbols
          (e.g., ``solve(f, x, y)``)
        - ordered iterable of symbols
          (e.g., ``solve(f, [x, y])``)
    
    flags :
        dict=True (default is False)
            Return list (perhaps empty) of solution mappings.
        set=True (default is False)
            Return list of symbols and set of tuple(s) of solution(s).
        exclude=[] (default)
            Do not try to solve for any of the free symbols in exclude;
            if expressions are given, the free symbols in them will
            be extracted automatically.
        check=True (default)
            If False, do not do any testing of solutions. This can be
            useful if you want to include solutions that make any
            denominator zero.
        numerical=True (default)
            Do a fast numerical check if *f* has only one symbol.
        minimal=True (default is False)
            A very fast, minimal testing.
        warn=True (default is False)
            Show a warning if ``checksol()`` could not conclude.
        simplify=True (default)
            Simplify all but polynomials of order 3 or greater before
            returning them and (if check is not False) use the
            general simplify function on the solutions and the
            expression obtained when they are substituted into the
            function which should be zero.
        force=True (default is False)
            Make positive all symbols without assumptions regarding sign.
        rational=True (default)
            Recast Floats as Rational; if this option is not used, the
            system containing Floats may fail to solve because of issues
            with polys. If rational=None, Floats will be recast as
            rationals but the answer will be recast as Floats. If the
            flag is False then nothing will be done to the Floats.
        manual=True (default is False)
            Do not use the polys/matrix method to solve a system of
            equations, solve them one at a time as you might "manually."
        implicit=True (default is False)
            Allows ``solve`` to return a solution for a pattern in terms of
            other functions that contain that pattern; this is only
            needed if the pattern is inside of some invertible function
            like cos, exp, ect.
        particular=True (default is False)
            Instructs ``solve`` to try to find a particular solution to a linear
            system with as many zeros as possible; this is very expensive.
        quick=True (default is False)
            When using particular=True, use a fast heuristic to find a
            solution with many zeros (instead of using the very slow method
            guaranteed to find the largest number of zeros possible).
        cubics=True (default)
            Return explicit solutions when cubic expressions are encountered.
        quartics=True (default)
            Return explicit solutions when quartic expressions are encountered.
        quintics=True (default)
            Return explicit solutions (if possible) when quintic expressions
            are encountered.
    
    See Also
    ========
    
    rsolve: For solving recurrence relationships
    dsolve: For solving differential equations

eq1=smp.Eq(a*x+b*y, c)
eq1

$\displaystyle a x + b y = c$

eq2=smp.Eq(d*x+e*y, f)
eq2

$\displaystyle d x + e y = f$

smp.solve([eq1, eq2], (x, y))

$\displaystyle \left\{ x : \frac{- b f + c e}{a e - b d}, \ y : \frac{a f - c d}{a e - b d}\right\}$

smp.init_printing()

smp.solve([eq1, eq2], (x, y))

$\displaystyle \left\{ x : \frac{- b f + c e}{a e - b d}, \ y : \frac{a f - c d}{a e - b d}\right\}$

smp.solve([eq1, eq2], (x, y), set = True)

$\displaystyle \left( \left[ x, \ y\right], \ \left\{\left( \frac{- b f + c e}{a e - b d}, \ \frac{a f - c d}{a e - b d}\right)\right\}\right)$

a, b, c, d, e, f = [1, 2, 3, 4, 5, 6]

eq1=smp.Eq(a*x+b*y, c)
eq1

$\displaystyle x + 2 y = 3$

eq2=smp.Eq(d*x+e*y, f)
eq2

$\displaystyle 4 x + 5 y = 6$

smp.solve([eq1, eq2], (x, y))

$\displaystyle \left\{ x : -1, \ y : 2\right\}$

A = np.array([[1, 2], [4, 5]])
A

array([[1, 2],
       [4, 5]])

c=np.array([3, 6])
c

array([3, 6])

np.linalg.solve(A, c)

array([-1.,  2.])

x, y, z, a11, a12, a13, a21, a22, a23, a31, a32, a33, b1, b2, b3 = smp.symbols('x y z a11, a12, a13, a21, a22, a23, a31, a32, a33, b1, b2, b3', real=True)

eq1=smp.Eq(a11*x+a12*y+a13*z, b1)
eq1

$\displaystyle a_{11} x + a_{12} y + a_{13} z = b_{1}$

eq2=smp.Eq(a21*x+a22*y+a23*z, b2)
eq2

$\displaystyle a_{21} x + a_{22} y + a_{23} z = b_{2}$

eq3=smp.Eq(a31*x+a32*y+a33*z, b3)
eq3

$\displaystyle a_{31} x + a_{32} y + a_{33} z = b_{3}$

smp.solve([eq1, eq2, eq3], (x, y, z), set = True)

$\displaystyle \left( \left[ x, \ y, \ z\right], \ \left\{\left( \frac{a_{12} a_{23} b_{3} - a_{12} a_{33} b_{2} - a_{13} a_{22} b_{3} + a_{13} a_{32} b_{2} + a_{22} a_{33} b_{1} - a_{23} a_{32} b_{1}}{a_{11} a_{22} a_{33} - a_{11} a_{23} a_{32} - a_{12} a_{21} a_{33} + a_{12} a_{23} a_{31} + a_{13} a_{21} a_{32} - a_{13} a_{22} a_{31}}, \ \frac{- a_{11} a_{23} b_{3} + a_{11} a_{33} b_{2} + a_{13} a_{21} b_{3} - a_{13} a_{31} b_{2} - a_{21} a_{33} b_{1} + a_{23} a_{31} b_{1}}{a_{11} a_{22} a_{33} - a_{11} a_{23} a_{32} - a_{12} a_{21} a_{33} + a_{12} a_{23} a_{31} + a_{13} a_{21} a_{32} - a_{13} a_{22} a_{31}}, \ \frac{a_{11} a_{22} b_{3} - a_{11} a_{32} b_{2} - a_{12} a_{21} b_{3} + a_{12} a_{31} b_{2} + a_{21} a_{32} b_{1} - a_{22} a_{31} b_{1}}{a_{11} a_{22} a_{33} - a_{11} a_{23} a_{32} - a_{12} a_{21} a_{33} + a_{12} a_{23} a_{31} + a_{13} a_{21} a_{32} - a_{13} a_{22} a_{31}}\right)\right\}\right)$

a11, a12, a13, b1, a21, a22, a23, b2, a31, a32, a33, b3 =[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12]

eq1=smp.Eq(a11*x+a12*y+a13*z, b1)
eq1

$\displaystyle x + 2 y + 3 z = 4$

eq2=smp.Eq(a21*x+a22*y+a23*z, b2)
eq2

$\displaystyle 5 x + 6 y + 7 z = 8$

eq3=smp.Eq(a31*x+a32*y+a33*z, b3)
eq3

$\displaystyle 9 x + 10 y + 11 z = 12$

smp.solve([eq1, eq2, eq3], (x, y, z), set = True)

$\displaystyle \left( \left[ x, \ y\right], \ \left\{\left( z - 2, \ 3 - 2 z\right)\right\}\right)$

A=np.array([[1, 2, 3, ], [5, 6, 7], [9, 10, 11]])b

A

array([[ 1,  2,  3],
       [ 5,  6,  7],
       [ 9, 10, 11]])

b=np.array([4, 8, 12])
b

array([ 4,  8, 12])

np.linalg.solve(A, b)

array([-0.75,  0.5 ,  1.25])

1.25-2

$\displaystyle -0.75$

3-2*1.25

$\displaystyle 0.5$

x, y, z, A11, A12, A13, A21, A22, A23, A31, A32, A33, B1, B2, B3 = smp.symbols('x y z A11, A12, A13, A21, A22, A23, A31, A32, A33, B1, B2, B3', real=True)

eqns=[A11*x+A22*y+A13*z-B1, A21*x+A22*y+A23*z-B2, a31*x+a32*y+a33*z-B3]
A, B = smp.linear_eq_to_matrix(eqns, x, y, z)
A, B

$\displaystyle \left( \left[\begin{matrix}A_{11} & A_{22} & A_{13}\\A_{21} & A_{22} & A_{23}\\a_{31} & a_{32} & a_{33}\end{matrix}\right], \ \left[\begin{matrix}B_{1}\\B_{2}\\B_{3}\end{matrix}\right]\right)$

smp.linsolve((A, B), [x,y,z])

$\displaystyle \left\{\left( \frac{A_{13} A_{22} B_{3} - A_{13} B_{2} a_{32} - A_{22} A_{23} B_{3} - A_{22} B_{1} a_{33} + A_{22} B_{2} a_{33} + A_{23} B_{1} a_{32}}{- A_{11} A_{22} a_{33} + A_{11} A_{23} a_{32} - A_{13} A_{21} a_{32} + A_{13} A_{22} a_{31} + A_{21} A_{22} a_{33} - A_{22} A_{23} a_{31}}, \ \frac{A_{11} A_{23} B_{3} - A_{11} B_{2} a_{33} - A_{13} A_{21} B_{3} + A_{13} B_{2} a_{31} + A_{21} B_{1} a_{33} - A_{23} B_{1} a_{31}}{- A_{11} A_{22} a_{33} + A_{11} A_{23} a_{32} - A_{13} A_{21} a_{32} + A_{13} A_{22} a_{31} + A_{21} A_{22} a_{33} - A_{22} A_{23} a_{31}}, \ \frac{- A_{11} A_{22} B_{3} + A_{11} B_{2} a_{32} + A_{21} A_{22} B_{3} - A_{21} B_{1} a_{32} + A_{22} B_{1} a_{31} - A_{22} B_{2} a_{31}}{- A_{11} A_{22} a_{33} + A_{11} A_{23} a_{32} - A_{13} A_{21} a_{32} + A_{13} A_{22} a_{31} + A_{21} A_{22} a_{33} - A_{22} A_{23} a_{31}}\right)\right\}$

A = smp.Matrix([[1,2,3], [5,6,7], [9,10,11]])
B = smp.Matrix([4, 8, 12])

A, B

$\displaystyle \left( \left[\begin{matrix}1 & 2 & 3\\5 & 6 & 7\\9 & 10 & 11\end{matrix}\right], \ \left[\begin{matrix}4\\8\\12\end{matrix}\right]\right)$

smp.linsolve((A, B), [x,y,z])

$\displaystyle \left\{\left( z - 2, \ 3 - 2 z, \ z\right)\right\}$

smp.linsolve((A, B), [x, y, z])

$\displaystyle \left\{\left( z - 2, \ 3 - 2 z, \ z\right)\right\}$

aug = smp.Matrix([[1,2,3, 4], [5,6,7, 8], [9,10,11, 12]])
aug

$\displaystyle \left[\begin{matrix}1 & 2 & 3 & 4\\5 & 6 & 7 & 8\\9 & 10 & 11 & 12\end{matrix}\right]$

smp.linsolve(aug, x,y,z)

$\displaystyle \left\{\left( \frac{3}{10}, \ \frac{2}{5}, \ 0\right)\right\}$

smp.linsolve(aug, x, y, z)

$\displaystyle \left\{\left( \frac{3}{10}, \ \frac{2}{5}, \ 0\right)\right\}$

Eqns = [3*x + 2*y - z - 1, 2*x - 2*y + 4*z + 2, - x + y/2 - z]

smp.linsolve(Eqns, x, y, z)

$\displaystyle \left\{\left( 1, \ -2, \ -2\right)\right\}$

aug = smp.Matrix([[2, 1, 3, 1], [2, 6, 8, 3], [6, 8, 18, 5]])

aug

$\displaystyle \left[\begin{matrix}2 & 1 & 3 & 1\\2 & 6 & 8 & 3\\6 & 8 & 18 & 5\end{matrix}\right]$

smp.linsolve(aug, x, y, z)

$\displaystyle \left\{\left( \frac{3}{10}, \ \frac{2}{5}, \ 0\right)\right\}$

Derivatives

Derivatives by smp.diff(f, x, n), where f is a function of variable x, f(x), and n is the order of derivatives

Example 1:

x = smp.symbols('x', real = True)

f = x

f

$\displaystyle x$

dfdx = smp.diff(f, x)

dfdx

$\displaystyle 1$

d2fdx2=smp.diff(f, x, 2)

d2fdx2

$\displaystyle 0$

Example 2

f = x ** 2

dfdx = smp.diff(f,x)
dfdx

$\displaystyle 2 x$

d2fdx2=smp.diff(f, x, 2)

d2fdx2

$\displaystyle 2$

Example 3

f = smp.sin(2*x)
f

$\displaystyle \sin{\left(2 x \right)}$

dfdx = smp.diff(f, x)
dfdx

$\displaystyle 2 \cos{\left(2 x \right)}$

d2fdx2 = smp.diff(f, x, 2)
d2fdx2

$\displaystyle - 4 \sin{\left(2 x \right)}$

d3fdx3 = smp.diff(f, x, 3)
d3fdx3

$\displaystyle - 8 \cos{\left(2 x \right)}$

Example 4

x, a = smp.symbols('x a', real = True)

f = a * smp.exp(a * x)
f

$\displaystyle a e^{a x}$

dfdx = smp.diff(f, x)
dfdx

$\displaystyle a^{2} e^{a x}$

dfdx = smp.diff(f, a)
dfdx

$\displaystyle a x e^{a x} + e^{a x}$

d3fdx3 = smp.diff(f, a, 3)
d3fdx3

$\displaystyle x^{2} \left(a x + 3\right) e^{a x}$

Example 5: $$f(x) = e^{-a\sin(x^2)}\cdot \sin(b^x) \cdot \ln(c\sin^2(x)/x)$$

x, a, b, c = smp.symbols('x a b c', real=True)

f = smp.exp(-a*smp.sin(x**2)) * smp.sin(b**x) * smp.log(c*smp.sin(x)**2 /x)

f

$\displaystyle e^{- a \sin{\left(x^{2} \right)}} \log{\left(\frac{c \sin^{2}{\left(x \right)}}{x} \right)} \sin{\left(b^{x} \right)}$

dfdx = smp.diff(f, x)
dfdx

$\displaystyle - 2 a x e^{- a \sin{\left(x^{2} \right)}} \log{\left(\frac{c \sin^{2}{\left(x \right)}}{x} \right)} \sin{\left(b^{x} \right)} \cos{\left(x^{2} \right)} + b^{x} e^{- a \sin{\left(x^{2} \right)}} \log{\left(b \right)} \log{\left(\frac{c \sin^{2}{\left(x \right)}}{x} \right)} \cos{\left(b^{x} \right)} + \frac{x \left(\frac{2 c \sin{\left(x \right)} \cos{\left(x \right)}}{x} - \frac{c \sin^{2}{\left(x \right)}}{x^{2}}\right) e^{- a \sin{\left(x^{2} \right)}} \sin{\left(b^{x} \right)}}{c \sin^{2}{\left(x \right)}}$

d4fdx4 = smp.diff(f, x, 4)
d4fdx4

$\displaystyle \left(\frac{24 a b^{x} x \left(b^{x} \sin{\left(b^{x} \right)} - \cos{\left(b^{x} \right)}\right) \left(2 \cos{\left(x \right)} - \frac{\sin{\left(x \right)}}{x}\right) \log{\left(b \right)}^{2} \cos{\left(x^{2} \right)}}{\sin{\left(x \right)}} + 8 a b^{x} x \left(b^{2 x} \cos{\left(b^{x} \right)} + 3 b^{x} \sin{\left(b^{x} \right)} - \cos{\left(b^{x} \right)}\right) \log{\left(b \right)}^{3} \log{\left(\frac{c \sin^{2}{\left(x \right)}}{x} \right)} \cos{\left(x^{2} \right)} + \frac{24 a b^{x} x \left(\frac{2 \left(2 \cos{\left(x \right)} - \frac{\sin{\left(x \right)}}{x}\right) \cos{\left(x \right)}}{\sin{\left(x \right)}} + \frac{2 \left(\sin^{2}{\left(x \right)} - \cos^{2}{\left(x \right)} + \frac{2 \sin{\left(x \right)} \cos{\left(x \right)}}{x} - \frac{\sin^{2}{\left(x \right)}}{x^{2}}\right)}{\sin{\left(x \right)}} - \frac{2 \cos{\left(x \right)} - \frac{\sin{\left(x \right)}}{x}}{x}\right) \log{\left(b \right)} \cos{\left(b^{x} \right)} \cos{\left(x^{2} \right)}}{\sin{\left(x \right)}} + 16 a b^{x} x \left(- 2 a^{2} x^{2} \cos^{3}{\left(x^{2} \right)} - 6 a x^{2} \sin{\left(x^{2} \right)} \cos{\left(x^{2} \right)} + 3 a \cos^{2}{\left(x^{2} \right)} + 2 x^{2} \cos{\left(x^{2} \right)} + 3 \sin{\left(x^{2} \right)}\right) \log{\left(b \right)} \log{\left(\frac{c \sin^{2}{\left(x \right)}}{x} \right)} \cos{\left(b^{x} \right)} - 12 a b^{x} \left(b^{x} \sin{\left(b^{x} \right)} - \cos{\left(b^{x} \right)}\right) \left(2 a x^{2} \cos^{2}{\left(x^{2} \right)} + 2 x^{2} \sin{\left(x^{2} \right)} - \cos{\left(x^{2} \right)}\right) \log{\left(b \right)}^{2} \log{\left(\frac{c \sin^{2}{\left(x \right)}}{x} \right)} + \frac{24 a b^{x} \left(2 \cos{\left(x \right)} - \frac{\sin{\left(x \right)}}{x}\right) \left(2 a x^{2} \cos^{2}{\left(x^{2} \right)} + 2 x^{2} \sin{\left(x^{2} \right)} - \cos{\left(x^{2} \right)}\right) \log{\left(b \right)} \cos{\left(b^{x} \right)}}{\sin{\left(x \right)}} + \frac{16 a x \left(2 \cos{\left(x \right)} - \frac{\sin{\left(x \right)}}{x}\right) \left(- 2 a^{2} x^{2} \cos^{3}{\left(x^{2} \right)} - 6 a x^{2} \sin{\left(x^{2} \right)} \cos{\left(x^{2} \right)} + 3 a \cos^{2}{\left(x^{2} \right)} + 2 x^{2} \cos{\left(x^{2} \right)} + 3 \sin{\left(x^{2} \right)}\right) \sin{\left(b^{x} \right)}}{\sin{\left(x \right)}} + \frac{16 a x \left(- \frac{3 \left(2 \cos{\left(x \right)} - \frac{\sin{\left(x \right)}}{x}\right) \cos^{2}{\left(x \right)}}{\sin^{2}{\left(x \right)}} - \frac{4 \left(\sin^{2}{\left(x \right)} - \cos^{2}{\left(x \right)} + \frac{2 \sin{\left(x \right)} \cos{\left(x \right)}}{x} - \frac{\sin^{2}{\left(x \right)}}{x^{2}}\right) \cos{\left(x \right)}}{\sin^{2}{\left(x \right)}} + \frac{4 \sin{\left(x \right)} \cos{\left(x \right)} - \frac{3 \sin^{2}{\left(x \right)}}{x} + \frac{3 \cos^{2}{\left(x \right)}}{x} - \frac{6 \sin{\left(x \right)} \cos{\left(x \right)}}{x^{2}} + \frac{3 \sin^{2}{\left(x \right)}}{x^{3}}}{\sin{\left(x \right)}} - 2 \cos{\left(x \right)} + \frac{2 \left(2 \cos{\left(x \right)} - \frac{\sin{\left(x \right)}}{x}\right) \cos{\left(x \right)}}{x \sin{\left(x \right)}} + \frac{2 \left(\sin^{2}{\left(x \right)} - \cos^{2}{\left(x \right)} + \frac{2 \sin{\left(x \right)} \cos{\left(x \right)}}{x} - \frac{\sin^{2}{\left(x \right)}}{x^{2}}\right)}{x \sin{\left(x \right)}} + \frac{\sin{\left(x \right)}}{x}\right) \sin{\left(b^{x} \right)} \cos{\left(x^{2} \right)}}{\sin{\left(x \right)}} - \frac{12 a \left(2 a x^{2} \cos^{2}{\left(x^{2} \right)} + 2 x^{2} \sin{\left(x^{2} \right)} - \cos{\left(x^{2} \right)}\right) \left(\frac{2 \left(2 \cos{\left(x \right)} - \frac{\sin{\left(x \right)}}{x}\right) \cos{\left(x \right)}}{\sin{\left(x \right)}} + \frac{2 \left(\sin^{2}{\left(x \right)} - \cos^{2}{\left(x \right)} + \frac{2 \sin{\left(x \right)} \cos{\left(x \right)}}{x} - \frac{\sin^{2}{\left(x \right)}}{x^{2}}\right)}{\sin{\left(x \right)}} - \frac{2 \cos{\left(x \right)} - \frac{\sin{\left(x \right)}}{x}}{x}\right) \sin{\left(b^{x} \right)}}{\sin{\left(x \right)}} + 4 a \left(4 a^{3} x^{4} \cos^{4}{\left(x^{2} \right)} + 24 a^{2} x^{4} \sin{\left(x^{2} \right)} \cos^{2}{\left(x^{2} \right)} - 12 a^{2} x^{2} \cos^{3}{\left(x^{2} \right)} + 12 a x^{4} \sin^{2}{\left(x^{2} \right)} - 16 a x^{4} \cos^{2}{\left(x^{2} \right)} - 36 a x^{2} \sin{\left(x^{2} \right)} \cos{\left(x^{2} \right)} + 3 a \cos^{2}{\left(x^{2} \right)} - 4 x^{4} \sin{\left(x^{2} \right)} + 12 x^{2} \cos{\left(x^{2} \right)} + 3 \sin{\left(x^{2} \right)}\right) \log{\left(\frac{c \sin^{2}{\left(x \right)}}{x} \right)} \sin{\left(b^{x} \right)} + \frac{6 b^{x} \left(b^{x} \sin{\left(b^{x} \right)} - \cos{\left(b^{x} \right)}\right) \left(\frac{2 \left(2 \cos{\left(x \right)} - \frac{\sin{\left(x \right)}}{x}\right) \cos{\left(x \right)}}{\sin{\left(x \right)}} + \frac{2 \left(\sin^{2}{\left(x \right)} - \cos^{2}{\left(x \right)} + \frac{2 \sin{\left(x \right)} \cos{\left(x \right)}}{x} - \frac{\sin^{2}{\left(x \right)}}{x^{2}}\right)}{\sin{\left(x \right)}} - \frac{2 \cos{\left(x \right)} - \frac{\sin{\left(x \right)}}{x}}{x}\right) \log{\left(b \right)}^{2}}{\sin{\left(x \right)}} - \frac{4 b^{x} \left(2 \cos{\left(x \right)} - \frac{\sin{\left(x \right)}}{x}\right) \left(b^{2 x} \cos{\left(b^{x} \right)} + 3 b^{x} \sin{\left(b^{x} \right)} - \cos{\left(b^{x} \right)}\right) \log{\left(b \right)}^{3}}{\sin{\left(x \right)}} + b^{x} \left(b^{3 x} \sin{\left(b^{x} \right)} - 6 b^{2 x} \cos{\left(b^{x} \right)} - 7 b^{x} \sin{\left(b^{x} \right)} + \cos{\left(b^{x} \right)}\right) \log{\left(b \right)}^{4} \log{\left(\frac{c \sin^{2}{\left(x \right)}}{x} \right)} - \frac{8 b^{x} \left(- \frac{3 \left(2 \cos{\left(x \right)} - \frac{\sin{\left(x \right)}}{x}\right) \cos^{2}{\left(x \right)}}{\sin^{2}{\left(x \right)}} - \frac{4 \left(\sin^{2}{\left(x \right)} - \cos^{2}{\left(x \right)} + \frac{2 \sin{\left(x \right)} \cos{\left(x \right)}}{x} - \frac{\sin^{2}{\left(x \right)}}{x^{2}}\right) \cos{\left(x \right)}}{\sin^{2}{\left(x \right)}} + \frac{4 \sin{\left(x \right)} \cos{\left(x \right)} - \frac{3 \sin^{2}{\left(x \right)}}{x} + \frac{3 \cos^{2}{\left(x \right)}}{x} - \frac{6 \sin{\left(x \right)} \cos{\left(x \right)}}{x^{2}} + \frac{3 \sin^{2}{\left(x \right)}}{x^{3}}}{\sin{\left(x \right)}} - 2 \cos{\left(x \right)} + \frac{2 \left(2 \cos{\left(x \right)} - \frac{\sin{\left(x \right)}}{x}\right) \cos{\left(x \right)}}{x \sin{\left(x \right)}} + \frac{2 \left(\sin^{2}{\left(x \right)} - \cos^{2}{\left(x \right)} + \frac{2 \sin{\left(x \right)} \cos{\left(x \right)}}{x} - \frac{\sin^{2}{\left(x \right)}}{x^{2}}\right)}{x \sin{\left(x \right)}} + \frac{\sin{\left(x \right)}}{x}\right) \log{\left(b \right)} \cos{\left(b^{x} \right)}}{\sin{\left(x \right)}} - \frac{2 \left(\frac{8 \left(2 \cos{\left(x \right)} - \frac{\sin{\left(x \right)}}{x}\right) \cos{\left(x \right)}}{\sin{\left(x \right)}} + \frac{12 \left(2 \cos{\left(x \right)} - \frac{\sin{\left(x \right)}}{x}\right) \cos^{3}{\left(x \right)}}{\sin^{3}{\left(x \right)}} + \frac{6 \left(\sin^{2}{\left(x \right)} - \cos^{2}{\left(x \right)} + \frac{2 \sin{\left(x \right)} \cos{\left(x \right)}}{x} - \frac{\sin^{2}{\left(x \right)}}{x^{2}}\right)}{\sin{\left(x \right)}} + \frac{18 \left(\sin^{2}{\left(x \right)} - \cos^{2}{\left(x \right)} + \frac{2 \sin{\left(x \right)} \cos{\left(x \right)}}{x} - \frac{\sin^{2}{\left(x \right)}}{x^{2}}\right) \cos^{2}{\left(x \right)}}{\sin^{3}{\left(x \right)}} - \frac{6 \left(4 \sin{\left(x \right)} \cos{\left(x \right)} - \frac{3 \sin^{2}{\left(x \right)}}{x} + \frac{3 \cos^{2}{\left(x \right)}}{x} - \frac{6 \sin{\left(x \right)} \cos{\left(x \right)}}{x^{2}} + \frac{3 \sin^{2}{\left(x \right)}}{x^{3}}\right) \cos{\left(x \right)}}{\sin^{2}{\left(x \right)}} - \frac{4 \left(\sin^{2}{\left(x \right)} - \cos^{2}{\left(x \right)} + \frac{4 \sin{\left(x \right)} \cos{\left(x \right)}}{x} - \frac{3 \sin^{2}{\left(x \right)}}{x^{2}} + \frac{3 \cos^{2}{\left(x \right)}}{x^{2}} - \frac{6 \sin{\left(x \right)} \cos{\left(x \right)}}{x^{3}} + \frac{3 \sin^{2}{\left(x \right)}}{x^{4}}\right)}{\sin{\left(x \right)}} - \frac{3 \left(2 \cos{\left(x \right)} - \frac{\sin{\left(x \right)}}{x}\right)}{x} - \frac{9 \left(2 \cos{\left(x \right)} - \frac{\sin{\left(x \right)}}{x}\right) \cos^{2}{\left(x \right)}}{x \sin^{2}{\left(x \right)}} - \frac{12 \left(\sin^{2}{\left(x \right)} - \cos^{2}{\left(x \right)} + \frac{2 \sin{\left(x \right)} \cos{\left(x \right)}}{x} - \frac{\sin^{2}{\left(x \right)}}{x^{2}}\right) \cos{\left(x \right)}}{x \sin^{2}{\left(x \right)}} + \frac{3 \left(4 \sin{\left(x \right)} \cos{\left(x \right)} - \frac{3 \sin^{2}{\left(x \right)}}{x} + \frac{3 \cos^{2}{\left(x \right)}}{x} - \frac{6 \sin{\left(x \right)} \cos{\left(x \right)}}{x^{2}} + \frac{3 \sin^{2}{\left(x \right)}}{x^{3}}\right)}{x \sin{\left(x \right)}}\right) \sin{\left(b^{x} \right)}}{\sin{\left(x \right)}}\right) e^{- a \sin{\left(x^{2} \right)}}$

Finding an analytical expression for the derivative after substituting the values of x, a, b, and c into the derivative expression

d4fdx4.subs([(x,4),(a,1),(b,2),(c,3)])

$\displaystyle \left(512 \left(255 \cos{\left(16 \right)} + 48 \sin{\left(16 \right)}\right) \log{\left(2 \right)}^{3} \log{\left(\frac{3 \sin^{2}{\left(4 \right)}}{4} \right)} \cos{\left(16 \right)} + 1024 \left(32 \cos{\left(16 \right)} - 96 \sin{\left(16 \right)} \cos{\left(16 \right)} + 3 \sin{\left(16 \right)} + 3 \cos^{2}{\left(16 \right)} - 32 \cos^{3}{\left(16 \right)}\right) \log{\left(2 \right)} \log{\left(\frac{3 \sin^{2}{\left(4 \right)}}{4} \right)} \cos{\left(16 \right)} + \frac{384 \left(2 \cos{\left(4 \right)} - \frac{\sin{\left(4 \right)}}{4}\right) \left(32 \sin{\left(16 \right)} - \cos{\left(16 \right)} + 32 \cos^{2}{\left(16 \right)}\right) \log{\left(2 \right)} \cos{\left(16 \right)}}{\sin{\left(4 \right)}} - 192 \left(16 \sin{\left(16 \right)} - \cos{\left(16 \right)}\right) \left(32 \sin{\left(16 \right)} - \cos{\left(16 \right)} + 32 \cos^{2}{\left(16 \right)}\right) \log{\left(2 \right)}^{2} \log{\left(\frac{3 \sin^{2}{\left(4 \right)}}{4} \right)} + 4 \left(- 4093 \cos^{2}{\left(16 \right)} + 6144 \sin{\left(16 \right)} \cos^{2}{\left(16 \right)} + 192 \cos{\left(16 \right)} - 576 \sin{\left(16 \right)} \cos{\left(16 \right)} - 192 \cos^{3}{\left(16 \right)} + 3072 \sin^{2}{\left(16 \right)} - 1021 \sin{\left(16 \right)} + 1024 \cos^{4}{\left(16 \right)}\right) \log{\left(\frac{3 \sin^{2}{\left(4 \right)}}{4} \right)} \sin{\left(16 \right)} + 16 \left(3984 \sin{\left(16 \right)} - 1535 \cos{\left(16 \right)}\right) \log{\left(2 \right)}^{4} \log{\left(\frac{3 \sin^{2}{\left(4 \right)}}{4} \right)} + \frac{96 \left(16 \sin{\left(16 \right)} - \cos{\left(16 \right)}\right) \left(\frac{2 \left(2 \cos{\left(4 \right)} - \frac{\sin{\left(4 \right)}}{4}\right) \cos{\left(4 \right)}}{\sin{\left(4 \right)}} + \frac{2 \left(- \cos^{2}{\left(4 \right)} + \frac{\sin{\left(4 \right)} \cos{\left(4 \right)}}{2} + \frac{15 \sin^{2}{\left(4 \right)}}{16}\right)}{\sin{\left(4 \right)}} + \frac{\sin{\left(4 \right)}}{16} - \frac{\cos{\left(4 \right)}}{2}\right) \log{\left(2 \right)}^{2}}{\sin{\left(4 \right)}} - \frac{128 \left(\frac{- \frac{45 \sin^{2}{\left(4 \right)}}{64} + \frac{3 \cos^{2}{\left(4 \right)}}{4} + \frac{29 \sin{\left(4 \right)} \cos{\left(4 \right)}}{8}}{\sin{\left(4 \right)}} + \frac{\left(2 \cos{\left(4 \right)} - \frac{\sin{\left(4 \right)}}{4}\right) \cos{\left(4 \right)}}{2 \sin{\left(4 \right)}} + \frac{- \cos^{2}{\left(4 \right)} + \frac{\sin{\left(4 \right)} \cos{\left(4 \right)}}{2} + \frac{15 \sin^{2}{\left(4 \right)}}{16}}{2 \sin{\left(4 \right)}} + \frac{\sin{\left(4 \right)}}{4} - 2 \cos{\left(4 \right)} - \frac{4 \left(- \cos^{2}{\left(4 \right)} + \frac{\sin{\left(4 \right)} \cos{\left(4 \right)}}{2} + \frac{15 \sin^{2}{\left(4 \right)}}{16}\right) \cos{\left(4 \right)}}{\sin^{2}{\left(4 \right)}} - \frac{3 \left(2 \cos{\left(4 \right)} - \frac{\sin{\left(4 \right)}}{4}\right) \cos^{2}{\left(4 \right)}}{\sin^{2}{\left(4 \right)}}\right) \log{\left(2 \right)} \cos{\left(16 \right)}}{\sin{\left(4 \right)}} + \frac{64 \left(\frac{- \frac{45 \sin^{2}{\left(4 \right)}}{64} + \frac{3 \cos^{2}{\left(4 \right)}}{4} + \frac{29 \sin{\left(4 \right)} \cos{\left(4 \right)}}{8}}{\sin{\left(4 \right)}} + \frac{\left(2 \cos{\left(4 \right)} - \frac{\sin{\left(4 \right)}}{4}\right) \cos{\left(4 \right)}}{2 \sin{\left(4 \right)}} + \frac{- \cos^{2}{\left(4 \right)} + \frac{\sin{\left(4 \right)} \cos{\left(4 \right)}}{2} + \frac{15 \sin^{2}{\left(4 \right)}}{16}}{2 \sin{\left(4 \right)}} + \frac{\sin{\left(4 \right)}}{4} - 2 \cos{\left(4 \right)} - \frac{4 \left(- \cos^{2}{\left(4 \right)} + \frac{\sin{\left(4 \right)} \cos{\left(4 \right)}}{2} + \frac{15 \sin^{2}{\left(4 \right)}}{16}\right) \cos{\left(4 \right)}}{\sin^{2}{\left(4 \right)}} - \frac{3 \left(2 \cos{\left(4 \right)} - \frac{\sin{\left(4 \right)}}{4}\right) \cos^{2}{\left(4 \right)}}{\sin^{2}{\left(4 \right)}}\right) \sin{\left(16 \right)} \cos{\left(16 \right)}}{\sin{\left(4 \right)}} - \frac{2 \left(\frac{12 \left(2 \cos{\left(4 \right)} - \frac{\sin{\left(4 \right)}}{4}\right) \cos^{3}{\left(4 \right)}}{\sin^{3}{\left(4 \right)}} + \frac{8 \left(2 \cos{\left(4 \right)} - \frac{\sin{\left(4 \right)}}{4}\right) \cos{\left(4 \right)}}{\sin{\left(4 \right)}} + \frac{18 \left(- \cos^{2}{\left(4 \right)} + \frac{\sin{\left(4 \right)} \cos{\left(4 \right)}}{2} + \frac{15 \sin^{2}{\left(4 \right)}}{16}\right) \cos^{2}{\left(4 \right)}}{\sin^{3}{\left(4 \right)}} + \frac{6 \left(- \cos^{2}{\left(4 \right)} + \frac{\sin{\left(4 \right)} \cos{\left(4 \right)}}{2} + \frac{15 \sin^{2}{\left(4 \right)}}{16}\right)}{\sin{\left(4 \right)}} + \frac{3 \left(- \frac{45 \sin^{2}{\left(4 \right)}}{64} + \frac{3 \cos^{2}{\left(4 \right)}}{4} + \frac{29 \sin{\left(4 \right)} \cos{\left(4 \right)}}{8}\right)}{4 \sin{\left(4 \right)}} + \frac{3 \sin{\left(4 \right)}}{16} - \frac{3 \cos{\left(4 \right)}}{2} - \frac{3 \left(- \cos^{2}{\left(4 \right)} + \frac{\sin{\left(4 \right)} \cos{\left(4 \right)}}{2} + \frac{15 \sin^{2}{\left(4 \right)}}{16}\right) \cos{\left(4 \right)}}{\sin^{2}{\left(4 \right)}} - \frac{9 \left(2 \cos{\left(4 \right)} - \frac{\sin{\left(4 \right)}}{4}\right) \cos^{2}{\left(4 \right)}}{4 \sin^{2}{\left(4 \right)}} - \frac{4 \left(- \frac{13 \cos^{2}{\left(4 \right)}}{16} + \frac{29 \sin{\left(4 \right)} \cos{\left(4 \right)}}{32} + \frac{211 \sin^{2}{\left(4 \right)}}{256}\right)}{\sin{\left(4 \right)}} - \frac{6 \left(- \frac{45 \sin^{2}{\left(4 \right)}}{64} + \frac{3 \cos^{2}{\left(4 \right)}}{4} + \frac{29 \sin{\left(4 \right)} \cos{\left(4 \right)}}{8}\right) \cos{\left(4 \right)}}{\sin^{2}{\left(4 \right)}}\right) \sin{\left(16 \right)}}{\sin{\left(4 \right)}} - \frac{12 \left(32 \sin{\left(16 \right)} - \cos{\left(16 \right)} + 32 \cos^{2}{\left(16 \right)}\right) \left(\frac{2 \left(2 \cos{\left(4 \right)} - \frac{\sin{\left(4 \right)}}{4}\right) \cos{\left(4 \right)}}{\sin{\left(4 \right)}} + \frac{2 \left(- \cos^{2}{\left(4 \right)} + \frac{\sin{\left(4 \right)} \cos{\left(4 \right)}}{2} + \frac{15 \sin^{2}{\left(4 \right)}}{16}\right)}{\sin{\left(4 \right)}} + \frac{\sin{\left(4 \right)}}{16} - \frac{\cos{\left(4 \right)}}{2}\right) \sin{\left(16 \right)}}{\sin{\left(4 \right)}} + \frac{64 \left(2 \cos{\left(4 \right)} - \frac{\sin{\left(4 \right)}}{4}\right) \left(32 \cos{\left(16 \right)} - 96 \sin{\left(16 \right)} \cos{\left(16 \right)} + 3 \sin{\left(16 \right)} + 3 \cos^{2}{\left(16 \right)} - 32 \cos^{3}{\left(16 \right)}\right) \sin{\left(16 \right)}}{\sin{\left(4 \right)}} + \frac{1536 \left(\frac{2 \left(2 \cos{\left(4 \right)} - \frac{\sin{\left(4 \right)}}{4}\right) \cos{\left(4 \right)}}{\sin{\left(4 \right)}} + \frac{2 \left(- \cos^{2}{\left(4 \right)} + \frac{\sin{\left(4 \right)} \cos{\left(4 \right)}}{2} + \frac{15 \sin^{2}{\left(4 \right)}}{16}\right)}{\sin{\left(4 \right)}} + \frac{\sin{\left(4 \right)}}{16} - \frac{\cos{\left(4 \right)}}{2}\right) \log{\left(2 \right)} \cos^{2}{\left(16 \right)}}{\sin{\left(4 \right)}} + \frac{1536 \left(16 \sin{\left(16 \right)} - \cos{\left(16 \right)}\right) \left(2 \cos{\left(4 \right)} - \frac{\sin{\left(4 \right)}}{4}\right) \log{\left(2 \right)}^{2} \cos{\left(16 \right)}}{\sin{\left(4 \right)}} - \frac{64 \left(2 \cos{\left(4 \right)} - \frac{\sin{\left(4 \right)}}{4}\right) \left(255 \cos{\left(16 \right)} + 48 \sin{\left(16 \right)}\right) \log{\left(2 \right)}^{3}}{\sin{\left(4 \right)}}\right) e^{- \sin{\left(16 \right)}}$

Evaluating the derivatives numerically to find the final value

d4fdx4.subs([(x,4),(a,1),(b,2),(c,3)]).evalf()

$\displaystyle -73035.8044625845$

Making ready the numerical derivatives function for plotting

d4fdx4_f = smp.lambdify((x,a,b,c), d4fdx4)

help(smp.lambdify)

Help on function lambdify in module sympy.utilities.lambdify:

lambdify(args: Iterable, expr, modules=None, printer=None, use_imps=True, dummify=False)
    Convert a SymPy expression into a function that allows for fast
    numeric evaluation.
    
    .. warning::
       This function uses ``exec``, and thus shouldn't be used on
       unsanitized input.
    
    .. versionchanged:: 1.7.0
       Passing a set for the *args* parameter is deprecated as sets are
       unordered. Use an ordered iterable such as a list or tuple.
    
    Explanation
    ===========
    
    For example, to convert the SymPy expression ``sin(x) + cos(x)`` to an
    equivalent NumPy function that numerically evaluates it:
    
    >>> from sympy import sin, cos, symbols, lambdify
    >>> import numpy as np
    >>> x = symbols('x')
    >>> expr = sin(x) + cos(x)
    >>> expr
    sin(x) + cos(x)
    >>> f = lambdify(x, expr, 'numpy')
    >>> a = np.array([1, 2])
    >>> f(a)
    [1.38177329 0.49315059]
    
    The primary purpose of this function is to provide a bridge from SymPy
    expressions to numerical libraries such as NumPy, SciPy, NumExpr, mpmath,
    and tensorflow. In general, SymPy functions do not work with objects from
    other libraries, such as NumPy arrays, and functions from numeric
    libraries like NumPy or mpmath do not work on SymPy expressions.
    ``lambdify`` bridges the two by converting a SymPy expression to an
    equivalent numeric function.
    
    The basic workflow with ``lambdify`` is to first create a SymPy expression
    representing whatever mathematical function you wish to evaluate. This
    should be done using only SymPy functions and expressions. Then, use
    ``lambdify`` to convert this to an equivalent function for numerical
    evaluation. For instance, above we created ``expr`` using the SymPy symbol
    ``x`` and SymPy functions ``sin`` and ``cos``, then converted it to an
    equivalent NumPy function ``f``, and called it on a NumPy array ``a``.
    
    Parameters
    ==========
    
    args : List[Symbol]
        A variable or a list of variables whose nesting represents the
        nesting of the arguments that will be passed to the function.
    
        Variables can be symbols, undefined functions, or matrix symbols.
    
        >>> from sympy import Eq
        >>> from sympy.abc import x, y, z
    
        The list of variables should match the structure of how the
        arguments will be passed to the function. Simply enclose the
        parameters as they will be passed in a list.
    
        To call a function like ``f(x)`` then ``[x]``
        should be the first argument to ``lambdify``; for this
        case a single ``x`` can also be used:
    
        >>> f = lambdify(x, x + 1)
        >>> f(1)
        2
        >>> f = lambdify([x], x + 1)
        >>> f(1)
        2
    
        To call a function like ``f(x, y)`` then ``[x, y]`` will
        be the first argument of the ``lambdify``:
    
        >>> f = lambdify([x, y], x + y)
        >>> f(1, 1)
        2
    
        To call a function with a single 3-element tuple like
        ``f((x, y, z))`` then ``[(x, y, z)]`` will be the first
        argument of the ``lambdify``:
    
        >>> f = lambdify([(x, y, z)], Eq(z**2, x**2 + y**2))
        >>> f((3, 4, 5))
        True
    
        If two args will be passed and the first is a scalar but
        the second is a tuple with two arguments then the items
        in the list should match that structure:
    
        >>> f = lambdify([x, (y, z)], x + y + z)
        >>> f(1, (2, 3))
        6
    
    expr : Expr
        An expression, list of expressions, or matrix to be evaluated.
    
        Lists may be nested.
        If the expression is a list, the output will also be a list.
    
        >>> f = lambdify(x, [x, [x + 1, x + 2]])
        >>> f(1)
        [1, [2, 3]]
    
        If it is a matrix, an array will be returned (for the NumPy module).
    
        >>> from sympy import Matrix
        >>> f = lambdify(x, Matrix([x, x + 1]))
        >>> f(1)
        [[1]
        [2]]
    
        Note that the argument order here (variables then expression) is used
        to emulate the Python ``lambda`` keyword. ``lambdify(x, expr)`` works
        (roughly) like ``lambda x: expr``
        (see :ref:`lambdify-how-it-works` below).
    
    modules : str, optional
        Specifies the numeric library to use.
    
        If not specified, *modules* defaults to:
    
        - ``["scipy", "numpy"]`` if SciPy is installed
        - ``["numpy"]`` if only NumPy is installed
        - ``["math", "mpmath", "sympy"]`` if neither is installed.
    
        That is, SymPy functions are replaced as far as possible by
        either ``scipy`` or ``numpy`` functions if available, and Python's
        standard library ``math``, or ``mpmath`` functions otherwise.
    
        *modules* can be one of the following types:
    
        - The strings ``"math"``, ``"mpmath"``, ``"numpy"``, ``"numexpr"``,
          ``"scipy"``, ``"sympy"``, or ``"tensorflow"``. This uses the
          corresponding printer and namespace mapping for that module.
        - A module (e.g., ``math``). This uses the global namespace of the
          module. If the module is one of the above known modules, it will
          also use the corresponding printer and namespace mapping
          (i.e., ``modules=numpy`` is equivalent to ``modules="numpy"``).
        - A dictionary that maps names of SymPy functions to arbitrary
          functions
          (e.g., ``{'sin': custom_sin}``).
        - A list that contains a mix of the arguments above, with higher
          priority given to entries appearing first
          (e.g., to use the NumPy module but override the ``sin`` function
          with a custom version, you can use
          ``[{'sin': custom_sin}, 'numpy']``).
    
    dummify : bool, optional
        Whether or not the variables in the provided expression that are not
        valid Python identifiers are substituted with dummy symbols.
    
        This allows for undefined functions like ``Function('f')(t)`` to be
        supplied as arguments. By default, the variables are only dummified
        if they are not valid Python identifiers.
    
        Set ``dummify=True`` to replace all arguments with dummy symbols
        (if ``args`` is not a string) - for example, to ensure that the
        arguments do not redefine any built-in names.
    
    
    Examples
    ========
    
    >>> from sympy.utilities.lambdify import implemented_function
    >>> from sympy import sqrt, sin, Matrix
    >>> from sympy import Function
    >>> from sympy.abc import w, x, y, z
    
    >>> f = lambdify(x, x**2)
    >>> f(2)
    4
    >>> f = lambdify((x, y, z), [z, y, x])
    >>> f(1,2,3)
    [3, 2, 1]
    >>> f = lambdify(x, sqrt(x))
    >>> f(4)
    2.0
    >>> f = lambdify((x, y), sin(x*y)**2)
    >>> f(0, 5)
    0.0
    >>> row = lambdify((x, y), Matrix((x, x + y)).T, modules='sympy')
    >>> row(1, 2)
    Matrix([[1, 3]])
    
    ``lambdify`` can be used to translate SymPy expressions into mpmath
    functions. This may be preferable to using ``evalf`` (which uses mpmath on
    the backend) in some cases.
    
    >>> f = lambdify(x, sin(x), 'mpmath')
    >>> f(1)
    0.8414709848078965
    
    Tuple arguments are handled and the lambdified function should
    be called with the same type of arguments as were used to create
    the function:
    
    >>> f = lambdify((x, (y, z)), x + y)
    >>> f(1, (2, 4))
    3
    
    The ``flatten`` function can be used to always work with flattened
    arguments:
    
    >>> from sympy.utilities.iterables import flatten
    >>> args = w, (x, (y, z))
    >>> vals = 1, (2, (3, 4))
    >>> f = lambdify(flatten(args), w + x + y + z)
    >>> f(*flatten(vals))
    10
    
    Functions present in ``expr`` can also carry their own numerical
    implementations, in a callable attached to the ``_imp_`` attribute. This
    can be used with undefined functions using the ``implemented_function``
    factory:
    
    >>> f = implemented_function(Function('f'), lambda x: x+1)
    >>> func = lambdify(x, f(x))
    >>> func(4)
    5
    
    ``lambdify`` always prefers ``_imp_`` implementations to implementations
    in other namespaces, unless the ``use_imps`` input parameter is False.
    
    Usage with Tensorflow:
    
    >>> import tensorflow as tf
    >>> from sympy import Max, sin, lambdify
    >>> from sympy.abc import x
    
    >>> f = Max(x, sin(x))
    >>> func = lambdify(x, f, 'tensorflow')
    
    After tensorflow v2, eager execution is enabled by default.
    If you want to get the compatible result across tensorflow v1 and v2
    as same as this tutorial, run this line.
    
    >>> tf.compat.v1.enable_eager_execution()
    
    If you have eager execution enabled, you can get the result out
    immediately as you can use numpy.
    
    If you pass tensorflow objects, you may get an ``EagerTensor``
    object instead of value.
    
    >>> result = func(tf.constant(1.0))
    >>> print(result)
    tf.Tensor(1.0, shape=(), dtype=float32)
    >>> print(result.__class__)
    <class 'tensorflow.python.framework.ops.EagerTensor'>
    
    You can use ``.numpy()`` to get the numpy value of the tensor.
    
    >>> result.numpy()
    1.0
    
    >>> var = tf.Variable(2.0)
    >>> result = func(var) # also works for tf.Variable and tf.Placeholder
    >>> result.numpy()
    2.0
    
    And it works with any shape array.
    
    >>> tensor = tf.constant([[1.0, 2.0], [3.0, 4.0]])
    >>> result = func(tensor)
    >>> result.numpy()
    [[1. 2.]
     [3. 4.]]
    
    Notes
    =====
    
    - For functions involving large array calculations, numexpr can provide a
      significant speedup over numpy. Please note that the available functions
      for numexpr are more limited than numpy but can be expanded with
      ``implemented_function`` and user defined subclasses of Function. If
      specified, numexpr may be the only option in modules. The official list
      of numexpr functions can be found at:
      https://numexpr.readthedocs.io/en/latest/user_guide.html#supported-functions
    
    - In previous versions of SymPy, ``lambdify`` replaced ``Matrix`` with
      ``numpy.matrix`` by default. As of SymPy 1.0 ``numpy.array`` is the
      default. To get the old default behavior you must pass in
      ``[{'ImmutableDenseMatrix':  numpy.matrix}, 'numpy']`` to the
      ``modules`` kwarg.
    
      >>> from sympy import lambdify, Matrix
      >>> from sympy.abc import x, y
      >>> import numpy
      >>> array2mat = [{'ImmutableDenseMatrix': numpy.matrix}, 'numpy']
      >>> f = lambdify((x, y), Matrix([x, y]), modules=array2mat)
      >>> f(1, 2)
      [[1]
       [2]]
    
    - In the above examples, the generated functions can accept scalar
      values or numpy arrays as arguments.  However, in some cases
      the generated function relies on the input being a numpy array:
    
      >>> from sympy import Piecewise
      >>> from sympy.testing.pytest import ignore_warnings
      >>> f = lambdify(x, Piecewise((x, x <= 1), (1/x, x > 1)), "numpy")
    
      >>> with ignore_warnings(RuntimeWarning):
      ...     f(numpy.array([-1, 0, 1, 2]))
      [-1.   0.   1.   0.5]
    
      >>> f(0)
      Traceback (most recent call last):
          ...
      ZeroDivisionError: division by zero
    
      In such cases, the input should be wrapped in a numpy array:
    
      >>> with ignore_warnings(RuntimeWarning):
      ...     float(f(numpy.array([0])))
      0.0
    
      Or if numpy functionality is not required another module can be used:
    
      >>> f = lambdify(x, Piecewise((x, x <= 1), (1/x, x > 1)), "math")
      >>> f(0)
      0
    
    .. _lambdify-how-it-works:
    
    How it works
    ============
    
    When using this function, it helps a great deal to have an idea of what it
    is doing. At its core, lambdify is nothing more than a namespace
    translation, on top of a special printer that makes some corner cases work
    properly.
    
    To understand lambdify, first we must properly understand how Python
    namespaces work. Say we had two files. One called ``sin_cos_sympy.py``,
    with
    
    .. code:: python
    
        # sin_cos_sympy.py
    
        from sympy import sin, cos
    
        def sin_cos(x):
            return sin(x) + cos(x)
    
    
    and one called ``sin_cos_numpy.py`` with
    
    .. code:: python
    
        # sin_cos_numpy.py
    
        from numpy import sin, cos
    
        def sin_cos(x):
            return sin(x) + cos(x)
    
    The two files define an identical function ``sin_cos``. However, in the
    first file, ``sin`` and ``cos`` are defined as the SymPy ``sin`` and
    ``cos``. In the second, they are defined as the NumPy versions.
    
    If we were to import the first file and use the ``sin_cos`` function, we
    would get something like
    
    >>> from sin_cos_sympy import sin_cos # doctest: +SKIP
    >>> sin_cos(1) # doctest: +SKIP
    cos(1) + sin(1)
    
    On the other hand, if we imported ``sin_cos`` from the second file, we
    would get
    
    >>> from sin_cos_numpy import sin_cos # doctest: +SKIP
    >>> sin_cos(1) # doctest: +SKIP
    1.38177329068
    
    In the first case we got a symbolic output, because it used the symbolic
    ``sin`` and ``cos`` functions from SymPy. In the second, we got a numeric
    result, because ``sin_cos`` used the numeric ``sin`` and ``cos`` functions
    from NumPy. But notice that the versions of ``sin`` and ``cos`` that were
    used was not inherent to the ``sin_cos`` function definition. Both
    ``sin_cos`` definitions are exactly the same. Rather, it was based on the
    names defined at the module where the ``sin_cos`` function was defined.
    
    The key point here is that when function in Python references a name that
    is not defined in the function, that name is looked up in the "global"
    namespace of the module where that function is defined.
    
    Now, in Python, we can emulate this behavior without actually writing a
    file to disk using the ``exec`` function. ``exec`` takes a string
    containing a block of Python code, and a dictionary that should contain
    the global variables of the module. It then executes the code "in" that
    dictionary, as if it were the module globals. The following is equivalent
    to the ``sin_cos`` defined in ``sin_cos_sympy.py``:
    
    >>> import sympy
    >>> module_dictionary = {'sin': sympy.sin, 'cos': sympy.cos}
    >>> exec('''
    ... def sin_cos(x):
    ...     return sin(x) + cos(x)
    ... ''', module_dictionary)
    >>> sin_cos = module_dictionary['sin_cos']
    >>> sin_cos(1)
    cos(1) + sin(1)
    
    and similarly with ``sin_cos_numpy``:
    
    >>> import numpy
    >>> module_dictionary = {'sin': numpy.sin, 'cos': numpy.cos}
    >>> exec('''
    ... def sin_cos(x):
    ...     return sin(x) + cos(x)
    ... ''', module_dictionary)
    >>> sin_cos = module_dictionary['sin_cos']
    >>> sin_cos(1)
    1.38177329068
    
    So now we can get an idea of how ``lambdify`` works. The name "lambdify"
    comes from the fact that we can think of something like ``lambdify(x,
    sin(x) + cos(x), 'numpy')`` as ``lambda x: sin(x) + cos(x)``, where
    ``sin`` and ``cos`` come from the ``numpy`` namespace. This is also why
    the symbols argument is first in ``lambdify``, as opposed to most SymPy
    functions where it comes after the expression: to better mimic the
    ``lambda`` keyword.
    
    ``lambdify`` takes the input expression (like ``sin(x) + cos(x)``) and
    
    1. Converts it to a string
    2. Creates a module globals dictionary based on the modules that are
       passed in (by default, it uses the NumPy module)
    3. Creates the string ``"def func({vars}): return {expr}"``, where ``{vars}`` is the
       list of variables separated by commas, and ``{expr}`` is the string
       created in step 1., then ``exec``s that string with the module globals
       namespace and returns ``func``.
    
    In fact, functions returned by ``lambdify`` support inspection. So you can
    see exactly how they are defined by using ``inspect.getsource``, or ``??`` if you
    are using IPython or the Jupyter notebook.
    
    >>> f = lambdify(x, sin(x) + cos(x))
    >>> import inspect
    >>> print(inspect.getsource(f))
    def _lambdifygenerated(x):
        return (sin(x) + cos(x))
    
    This shows us the source code of the function, but not the namespace it
    was defined in. We can inspect that by looking at the ``__globals__``
    attribute of ``f``:
    
    >>> f.__globals__['sin']
    <ufunc 'sin'>
    >>> f.__globals__['cos']
    <ufunc 'cos'>
    >>> f.__globals__['sin'] is numpy.sin
    True
    
    This shows us that ``sin`` and ``cos`` in the namespace of ``f`` will be
    ``numpy.sin`` and ``numpy.cos``.
    
    Note that there are some convenience layers in each of these steps, but at
    the core, this is how ``lambdify`` works. Step 1 is done using the
    ``LambdaPrinter`` printers defined in the printing module (see
    :mod:`sympy.printing.lambdarepr`). This allows different SymPy expressions
    to define how they should be converted to a string for different modules.
    You can change which printer ``lambdify`` uses by passing a custom printer
    in to the ``printer`` argument.
    
    Step 2 is augmented by certain translations. There are default
    translations for each module, but you can provide your own by passing a
    list to the ``modules`` argument. For instance,
    
    >>> def mysin(x):
    ...     print('taking the sin of', x)
    ...     return numpy.sin(x)
    ...
    >>> f = lambdify(x, sin(x), [{'sin': mysin}, 'numpy'])
    >>> f(1)
    taking the sin of 1
    0.8414709848078965
    
    The globals dictionary is generated from the list by merging the
    dictionary ``{'sin': mysin}`` and the module dictionary for NumPy. The
    merging is done so that earlier items take precedence, which is why
    ``mysin`` is used above instead of ``numpy.sin``.
    
    If you want to modify the way ``lambdify`` works for a given function, it
    is usually easiest to do so by modifying the globals dictionary as such.
    In more complicated cases, it may be necessary to create and pass in a
    custom printer.
    
    Finally, step 3 is augmented with certain convenience operations, such as
    the addition of a docstring.
    
    Understanding how ``lambdify`` works can make it easier to avoid certain
    gotchas when using it. For instance, a common mistake is to create a
    lambdified function for one module (say, NumPy), and pass it objects from
    another (say, a SymPy expression).
    
    For instance, say we create
    
    >>> from sympy.abc import x
    >>> f = lambdify(x, x + 1, 'numpy')
    
    Now if we pass in a NumPy array, we get that array plus 1
    
    >>> import numpy
    >>> a = numpy.array([1, 2])
    >>> f(a)
    [2 3]
    
    But what happens if you make the mistake of passing in a SymPy expression
    instead of a NumPy array:
    
    >>> f(x + 1)
    x + 2
    
    This worked, but it was only by accident. Now take a different lambdified
    function:
    
    >>> from sympy import sin
    >>> g = lambdify(x, x + sin(x), 'numpy')
    
    This works as expected on NumPy arrays:
    
    >>> g(a)
    [1.84147098 2.90929743]
    
    But if we try to pass in a SymPy expression, it fails
    
    >>> try:
    ...     g(x + 1)
    ... # NumPy release after 1.17 raises TypeError instead of
    ... # AttributeError
    ... except (AttributeError, TypeError):
    ...     raise AttributeError() # doctest: +IGNORE_EXCEPTION_DETAIL
    Traceback (most recent call last):
    ...
    AttributeError:
    
    Now, let's look at what happened. The reason this fails is that ``g``
    calls ``numpy.sin`` on the input expression, and ``numpy.sin`` does not
    know how to operate on a SymPy object. **As a general rule, NumPy
    functions do not know how to operate on SymPy expressions, and SymPy
    functions do not know how to operate on NumPy arrays. This is why lambdify
    exists: to provide a bridge between SymPy and NumPy.**
    
    However, why is it that ``f`` did work? That's because ``f`` doesn't call
    any functions, it only adds 1. So the resulting function that is created,
    ``def _lambdifygenerated(x): return x + 1`` does not depend on the globals
    namespace it is defined in. Thus it works, but only by accident. A future
    version of ``lambdify`` may remove this behavior.
    
    Be aware that certain implementation details described here may change in
    future versions of SymPy. The API of passing in custom modules and
    printers will not change, but the details of how a lambda function is
    created may change. However, the basic idea will remain the same, and
    understanding it will be helpful to understanding the behavior of
    lambdify.
    
    **In general: you should create lambdified functions for one module (say,
    NumPy), and only pass it input types that are compatible with that module
    (say, NumPy arrays).** Remember that by default, if the ``module``
    argument is not provided, ``lambdify`` creates functions using the NumPy
    and SciPy namespaces.

x = np.linspace(1,2,100)

y = d4fdx4_f(x, a=1, b=2, c=3)

plt.style.use(['seaborn', 'notebook', 'grid'])

plt.plot(x,y, c='orange')
plt.ylabel('$d^4 f / dx^4$', fontsize=24)
plt.xlabel('$x$', fontsize=24)

Text(0.5, 0, '$x$')

Numerical Case

x, y = np.loadtxt('data/sample_data1.txt')

plt.plot(x, y, c = 'green')

[<matplotlib.lines.Line2D at 0x7f87c893ec10>]

plt.plot(x, y, 'o', c = 'green')

[<matplotlib.lines.Line2D at 0x7f87c8ba5be0>]

plt.scatter(x, y, c = 'green')

<matplotlib.collections.PathCollection at 0x7f87c8e86a60>

plt.plot(x, y, 'o-', c = 'green')

[<matplotlib.lines.Line2D at 0x7f87c8e54a60>]

plt.plot(x, y, 'o--', c = 'green')

[<matplotlib.lines.Line2D at 0x7f87e61e89a0>]

1) The basic way, which works fine if the data is smooth but not if the data is noisy

dydx = np.gradient(y,x)

plt.plot(x,y, 'o--', label='$y(x)$', c = 'g', zorder = 0)
plt.plot(x,dydx, 'o--', label='$y^\prime (x)$', c='orange', zorder = 1)
plt.legend()

<matplotlib.legend.Legend at 0x7f87c76e3b20>

x, y = np.loadtxt('data/sample_data2.txt')
dydx = np.gradient(y,x)

plt.plot(x, y, 'o--', c='m', alpha = 0.7)

[<matplotlib.lines.Line2D at 0x7f87c8fe4b50>]

For the noisey data, this method does not work, beacuse noises are amplified in the derivative.

fig, ax = plt.subplots(1, 2, figsize=(10,3))
ax[0].plot(x,y, label='$y(x)$', c = 'orange')
ax[1].plot(x,dydx, label='$y\'(x)$', color='r')
ax[0].legend()
ax[1].legend()
plt.show()

fig, ax = plt.subplots(1, 2, figsize=(10,3))
ax[0].plot(x,y, label='$y(x)$', c = 'orange')
ax[1].plot(x,dydx, label='$y\'(x)$', color='r')
[a.legend() for a in ax]
plt.show()

Thus, let us introduce the second method:

2. The Clever Way

Smooth data then take derivative. Consider the following new covid cases per day data

x, y = np.loadtxt('data/coviddata.txt')
dydx = np.gradient(y,x)

Taking the derivative naively gives a bad result

fig, ax = plt.subplots(1, 2, figsize=(10,3))
ax[0].plot(x,y, label='$y(x)$')
ax[1].plot(x,dydx, label='$y\'(x)$', color='r')
[a.legend() for a in ax]
plt.show()

Smooth the data by convolving it with a rectangle

filt = np.ones(15)/15

np.ones(5)

array([1., 1., 1., 1., 1.])

filt

array([0.06666667, 0.06666667, 0.06666667, 0.06666667, 0.06666667,
       0.06666667, 0.06666667, 0.06666667, 0.06666667, 0.06666667,
       0.06666667, 0.06666667, 0.06666667, 0.06666667, 0.06666667])

1/15

$\displaystyle 0.0666666666666667$

y_smooth = np.convolve(y, filt, mode='valid')
dysdx = np.gradient(y_smooth, x[7:-7])

x[7:-7]

array([  7.,   8.,   9.,  10.,  11.,  12.,  13.,  14.,  15.,  16.,  17.,
        18.,  19.,  20.,  21.,  22.,  23.,  24.,  25.,  26.,  27.,  28.,
        29.,  30.,  31.,  32.,  33.,  34.,  35.,  36.,  37.,  38.,  39.,
        40.,  41.,  42.,  43.,  44.,  45.,  46.,  47.,  48.,  49.,  50.,
        51.,  52.,  53.,  54.,  55.,  56.,  57.,  58.,  59.,  60.,  61.,
        62.,  63.,  64.,  65.,  66.,  67.,  68.,  69.,  70.,  71.,  72.,
        73.,  74.,  75.,  76.,  77.,  78.,  79.,  80.,  81.,  82.,  83.,
        84.,  85.,  86.,  87.,  88.,  89.,  90.,  91.,  92.,  93.,  94.,
        95.,  96.,  97.,  98.,  99., 100., 101., 102., 103., 104., 105.,
       106., 107., 108., 109., 110., 111., 112., 113., 114., 115., 116.,
       117., 118., 119., 120., 121., 122., 123., 124., 125., 126., 127.,
       128., 129., 130., 131., 132., 133., 134., 135., 136., 137., 138.,
       139., 140., 141., 142., 143., 144., 145., 146., 147., 148., 149.,
       150., 151., 152., 153., 154., 155., 156., 157., 158., 159., 160.,
       161., 162., 163., 164., 165., 166., 167., 168., 169., 170., 171.,
       172., 173., 174., 175., 176., 177., 178., 179., 180., 181., 182.,
       183., 184., 185., 186., 187., 188., 189., 190., 191., 192., 193.,
       194., 195., 196., 197., 198., 199., 200., 201., 202., 203., 204.,
       205., 206., 207., 208., 209., 210., 211., 212., 213., 214., 215.,
       216., 217., 218., 219., 220., 221., 222., 223., 224., 225., 226.,
       227., 228., 229., 230., 231., 232., 233., 234., 235., 236., 237.,
       238., 239., 240., 241., 242., 243., 244., 245., 246., 247., 248.,
       249., 250., 251., 252., 253., 254., 255., 256., 257., 258., 259.,
       260., 261., 262., 263., 264., 265., 266., 267., 268., 269., 270.,
       271., 272., 273., 274., 275., 276., 277., 278., 279., 280., 281.,
       282., 283., 284., 285., 286., 287., 288., 289., 290., 291., 292.,
       293., 294., 295., 296., 297., 298., 299., 300., 301., 302., 303.,
       304., 305., 306., 307., 308., 309., 310., 311., 312., 313., 314.,
       315., 316., 317., 318., 319., 320., 321., 322., 323., 324., 325.,
       326., 327., 328., 329., 330., 331., 332., 333., 334., 335., 336.,
       337., 338., 339., 340., 341., 342., 343., 344., 345., 346., 347.,
       348., 349., 350., 351., 352., 353., 354., 355., 356., 357., 358.,
       359., 360., 361., 362., 363., 364., 365., 366., 367., 368., 369.,
       370., 371., 372., 373., 374., 375., 376., 377., 378., 379., 380.,
       381., 382., 383., 384., 385., 386., 387., 388., 389., 390., 391.,
       392., 393., 394., 395., 396., 397., 398., 399., 400., 401., 402.,
       403., 404., 405., 406., 407., 408., 409., 410., 411., 412., 413.,
       414., 415., 416., 417., 418., 419., 420., 421., 422., 423., 424.,
       425., 426., 427., 428., 429., 430., 431., 432., 433., 434., 435.,
       436., 437., 438., 439., 440., 441., 442., 443., 444., 445., 446.,
       447., 448., 449., 450., 451., 452., 453., 454., 455., 456., 457.,
       458., 459., 460., 461., 462., 463., 464., 465., 466., 467., 468.,
       469., 470., 471., 472., 473., 474., 475., 476., 477., 478., 479.,
       480., 481., 482., 483., 484., 485., 486., 487., 488., 489., 490.,
       491., 492., 493., 494., 495., 496., 497., 498., 499., 500., 501.,
       502., 503., 504., 505., 506., 507., 508., 509., 510., 511., 512.,
       513., 514., 515., 516., 517., 518., 519., 520., 521., 522., 523.,
       524., 525., 526., 527., 528., 529., 530., 531., 532., 533., 534.,
       535., 536., 537., 538., 539., 540., 541., 542., 543., 544.])

y_smooth

array([4.00000000e-01, 3.33333333e-01, 3.33333333e-01, 3.33333333e-01,
       3.33333333e-01, 3.33333333e-01, 5.33333333e-01, 5.33333333e-01,
       5.33333333e-01, 3.33333333e-01, 3.33333333e-01, 3.33333333e-01,
       7.33333333e-01, 6.00000000e-01, 6.00000000e-01, 6.00000000e-01,
       6.00000000e-01, 6.00000000e-01, 6.00000000e-01, 6.00000000e-01,
       1.60000000e+00, 1.40000000e+00, 2.60000000e+00, 4.00000000e+00,
       4.20000000e+00, 5.00000000e+00, 5.60000000e+00, 7.20000000e+00,
       1.04000000e+01, 1.20000000e+01, 1.52000000e+01, 2.08000000e+01,
       2.66000000e+01, 3.58666667e+01, 4.58666667e+01, 5.44666667e+01,
       6.57333333e+01, 7.71333333e+01, 9.04000000e+01, 1.00600000e+02,
       1.10600000e+02, 1.22533333e+02, 1.36733333e+02, 1.47133333e+02,
       1.57733333e+02, 1.67733333e+02, 1.71333333e+02, 1.75533333e+02,
       1.77733333e+02, 1.78333333e+02, 1.77933333e+02, 1.74066667e+02,
       1.67266667e+02, 1.55600000e+02, 1.55000000e+02, 1.51000000e+02,
       1.44666667e+02, 1.37266667e+02, 1.29666667e+02, 1.21666667e+02,
       1.11866667e+02, 1.08266667e+02, 1.06266667e+02, 9.98000000e+01,
       9.68000000e+01, 9.22000000e+01, 8.98000000e+01, 8.82000000e+01,
       9.40000000e+01, 9.48666667e+01, 9.40666667e+01, 9.98000000e+01,
       1.05066667e+02, 1.05466667e+02, 1.06466667e+02, 1.14133333e+02,
       1.14333333e+02, 1.12333333e+02, 1.13133333e+02, 1.12200000e+02,
       1.12200000e+02, 1.10666667e+02, 1.11133333e+02, 1.07333333e+02,
       9.98666667e+01, 9.86666667e+01, 8.95333333e+01, 7.80666667e+01,
       7.42666667e+01, 7.14666667e+01, 6.24000000e+01, 6.12666667e+01,
       5.68666667e+01, 5.46666667e+01, 5.02000000e+01, 4.72666667e+01,
       4.56666667e+01, 4.36666667e+01, 4.11333333e+01, 4.25333333e+01,
       3.93333333e+01, 3.81333333e+01, 3.72666667e+01, 3.60666667e+01,
       3.46666667e+01, 3.44666667e+01, 3.16000000e+01, 3.24000000e+01,
       3.12000000e+01, 3.10000000e+01, 3.17333333e+01, 3.24666667e+01,
       3.04000000e+01, 2.97333333e+01, 2.61333333e+01, 2.66000000e+01,
       2.66000000e+01, 2.60666667e+01, 2.76666667e+01, 2.98666667e+01,
       3.02000000e+01, 3.16000000e+01, 3.20666667e+01, 3.09333333e+01,
       3.05333333e+01, 3.24666667e+01, 3.09333333e+01, 3.15333333e+01,
       3.11333333e+01, 3.17333333e+01, 3.33333333e+01, 3.34000000e+01,
       3.60000000e+01, 3.66666667e+01, 3.48666667e+01, 3.44666667e+01,
       3.34666667e+01, 3.22000000e+01, 3.25333333e+01, 3.38000000e+01,
       3.15333333e+01, 3.21333333e+01, 3.36000000e+01, 3.48666667e+01,
       3.45333333e+01, 3.30666667e+01, 3.50000000e+01, 3.48000000e+01,
       3.51333333e+01, 3.93333333e+01, 4.20000000e+01, 4.34000000e+01,
       4.42666667e+01, 4.70666667e+01, 5.08000000e+01, 5.49333333e+01,
       5.98666667e+01, 6.21333333e+01, 6.55333333e+01, 7.00000000e+01,
       7.38000000e+01, 7.69333333e+01, 7.97333333e+01, 8.10000000e+01,
       7.96666667e+01, 7.94666667e+01, 7.98666667e+01, 8.78666667e+01,
       8.76666667e+01, 9.00000000e+01, 9.26000000e+01, 9.54000000e+01,
       9.87333333e+01, 9.95333333e+01, 1.00866667e+02, 1.02866667e+02,
       1.07133333e+02, 1.09133333e+02, 1.12466667e+02, 1.16133333e+02,
       1.24533333e+02, 1.35733333e+02, 1.41733333e+02, 1.56133333e+02,
       1.69333333e+02, 1.72733333e+02, 1.79200000e+02, 1.81133333e+02,
       1.88933333e+02, 1.98866667e+02, 2.12466667e+02, 2.25666667e+02,
       2.32133333e+02, 2.35733333e+02, 2.40000000e+02, 2.36466667e+02,
       2.40866667e+02, 2.46666667e+02, 2.49066667e+02, 2.47933333e+02,
       2.53000000e+02, 2.48933333e+02, 2.58133333e+02, 2.59666667e+02,
       2.72466667e+02, 2.75866667e+02, 2.73733333e+02, 2.81066667e+02,
       2.87066667e+02, 2.97666667e+02, 3.13000000e+02, 3.17600000e+02,
       3.25133333e+02, 3.24333333e+02, 3.20133333e+02, 3.19333333e+02,
       3.38200000e+02, 3.47200000e+02, 3.61600000e+02, 3.57800000e+02,
       3.54400000e+02, 3.57533333e+02, 3.54933333e+02, 3.52733333e+02,
       3.59600000e+02, 3.51200000e+02, 3.46600000e+02, 3.40000000e+02,
       3.37400000e+02, 3.39800000e+02, 3.45333333e+02, 3.38466667e+02,
       3.38866667e+02, 3.27266667e+02, 3.26266667e+02, 3.29333333e+02,
       3.22133333e+02, 3.24933333e+02, 3.33333333e+02, 3.29466667e+02,
       3.49866667e+02, 3.62266667e+02, 3.61466667e+02, 3.72266667e+02,
       3.84466667e+02, 3.80666667e+02, 3.87266667e+02, 4.00666667e+02,
       4.08866667e+02, 4.24666667e+02, 4.28000000e+02, 4.61200000e+02,
       4.98400000e+02, 5.14000000e+02, 5.43000000e+02, 5.65600000e+02,
       5.79400000e+02, 6.08800000e+02, 6.33400000e+02, 6.56266667e+02,
       7.06066667e+02, 7.52266667e+02, 7.86066667e+02, 8.33666667e+02,
       8.60266667e+02, 9.05666667e+02, 9.72666667e+02, 1.00926667e+03,
       1.05906667e+03, 1.11806667e+03, 1.15953333e+03, 1.21313333e+03,
       1.30113333e+03, 1.38013333e+03, 1.44286667e+03, 1.48346667e+03,
       1.59586667e+03, 1.63766667e+03, 1.67286667e+03, 1.75206667e+03,
       1.85066667e+03, 1.86086667e+03, 1.95146667e+03, 2.00626667e+03,
       2.02333333e+03, 2.05946667e+03, 2.12846667e+03, 2.16426667e+03,
       2.18386667e+03, 2.20786667e+03, 2.26293333e+03, 2.20153333e+03,
       2.19260000e+03, 2.23740000e+03, 2.24220000e+03, 2.20780000e+03,
       2.19693333e+03, 2.16560000e+03, 2.13526667e+03, 2.12060000e+03,
       2.11580000e+03, 2.10026667e+03, 2.07946667e+03, 2.07066667e+03,
       2.07673333e+03, 2.02606667e+03, 2.00160000e+03, 2.01000000e+03,
       1.98206667e+03, 1.95506667e+03, 1.94546667e+03, 1.94353333e+03,
       1.88886667e+03, 1.85526667e+03, 1.84146667e+03, 1.81686667e+03,
       1.77480000e+03, 1.70100000e+03, 1.62440000e+03, 1.57253333e+03,
       1.54033333e+03, 1.56080000e+03, 1.56233333e+03, 1.55966667e+03,
       1.55446667e+03, 1.50886667e+03, 1.50026667e+03, 1.50286667e+03,
       1.53540000e+03, 1.55406667e+03, 1.55600000e+03, 1.55066667e+03,
       1.57560000e+03, 1.59920000e+03, 1.58566667e+03, 1.60546667e+03,
       1.58626667e+03, 1.56646667e+03, 1.55666667e+03, 1.52960000e+03,
       1.51053333e+03, 1.47933333e+03, 1.47820000e+03, 1.47660000e+03,
       1.45953333e+03, 1.44733333e+03, 1.44066667e+03, 1.40653333e+03,
       1.38360000e+03, 1.40933333e+03, 1.40873333e+03, 1.40873333e+03,
       1.39366667e+03, 1.37486667e+03, 1.34153333e+03, 1.35766667e+03,
       1.36226667e+03, 1.35440000e+03, 1.33780000e+03, 1.32240000e+03,
       1.32193333e+03, 1.29246667e+03, 1.29946667e+03, 1.32746667e+03,
       1.30993333e+03, 1.29653333e+03, 1.28240000e+03, 1.28806667e+03,
       1.26446667e+03, 1.27013333e+03, 1.26893333e+03, 1.31173333e+03,
       1.32333333e+03, 1.33433333e+03, 1.34400000e+03, 1.33840000e+03,
       1.36886667e+03, 1.39386667e+03, 1.40020000e+03, 1.42893333e+03,
       1.44566667e+03, 1.45400000e+03, 1.45440000e+03, 1.47380000e+03,
       1.52080000e+03, 1.55513333e+03, 1.55733333e+03, 1.56053333e+03,
       1.56586667e+03, 1.54740000e+03, 1.55593333e+03, 1.56180000e+03,
       1.57400000e+03, 1.60460000e+03, 1.59600000e+03, 1.59773333e+03,
       1.59560000e+03, 1.60200000e+03, 1.61300000e+03, 1.62093333e+03,
       1.65533333e+03, 1.64533333e+03, 1.65866667e+03, 1.67933333e+03,
       1.71633333e+03, 1.76100000e+03, 1.81613333e+03, 1.89793333e+03,
       1.94300000e+03, 1.99800000e+03, 2.05333333e+03, 2.11393333e+03,
       2.22200000e+03, 2.29566667e+03, 2.39093333e+03, 2.46493333e+03,
       2.55500000e+03, 2.61486667e+03, 2.66846667e+03, 2.75846667e+03,
       2.87793333e+03, 2.96073333e+03, 3.02313333e+03, 3.05886667e+03,
       3.10166667e+03, 3.12920000e+03, 3.17293333e+03, 3.21326667e+03,
       3.22300000e+03, 3.22533333e+03, 3.19773333e+03, 3.19606667e+03,
       3.16306667e+03, 3.16346667e+03, 3.14866667e+03, 3.08760000e+03,
       3.02220000e+03, 2.93666667e+03, 2.87726667e+03, 2.82446667e+03,
       2.79593333e+03, 2.76920000e+03, 2.67380000e+03, 2.66640000e+03,
       2.58720000e+03, 2.52280000e+03, 2.43260000e+03, 2.40900000e+03,
       2.36520000e+03, 2.30326667e+03, 2.22386667e+03, 2.17046667e+03,
       2.10960000e+03, 2.05960000e+03, 2.01940000e+03, 1.98080000e+03,
       1.90253333e+03, 1.84353333e+03, 1.76213333e+03, 1.71233333e+03,
       1.66273333e+03, 1.65373333e+03, 1.60280000e+03, 1.54720000e+03,
       1.48773333e+03, 1.43153333e+03, 1.36713333e+03, 1.30393333e+03,
       1.25133333e+03, 1.21513333e+03, 1.15720000e+03, 1.11160000e+03,
       1.07026667e+03, 1.01686667e+03, 9.59666667e+02, 9.19866667e+02,
       8.55066667e+02, 8.15866667e+02, 7.74666667e+02, 7.31666667e+02,
       6.91600000e+02, 6.62000000e+02, 6.45400000e+02, 6.26333333e+02,
       5.76533333e+02, 5.35266667e+02, 5.08666667e+02, 4.73800000e+02,
       4.51400000e+02, 4.42400000e+02, 4.27600000e+02, 4.10666667e+02,
       3.93400000e+02, 3.66266667e+02, 3.47666667e+02, 3.32266667e+02,
       3.18066667e+02, 2.98866667e+02, 2.82000000e+02, 2.64000000e+02,
       2.52866667e+02, 2.37133333e+02, 2.29333333e+02, 2.17200000e+02,
       2.05000000e+02, 1.87266667e+02, 1.70600000e+02, 1.59666667e+02,
       1.46333333e+02, 1.46266667e+02, 1.43866667e+02, 1.40466667e+02,
       1.35866667e+02, 1.31400000e+02, 1.28266667e+02, 1.23666667e+02,
       1.24200000e+02, 1.24133333e+02, 1.26266667e+02, 1.26066667e+02,
       1.29800000e+02, 1.36400000e+02, 1.38533333e+02, 1.46600000e+02,
       1.53266667e+02, 1.63666667e+02, 1.73266667e+02, 1.81066667e+02,
       1.88066667e+02, 2.00000000e+02, 2.21066667e+02, 2.51266667e+02,
       2.85666667e+02, 3.21333333e+02, 3.42933333e+02, 3.71333333e+02,
       3.95533333e+02, 4.26133333e+02])

y

array([1.000e+00, 0.000e+00, 0.000e+00, 0.000e+00, 0.000e+00, 0.000e+00,
       0.000e+00, 0.000e+00, 3.000e+00, 0.000e+00, 0.000e+00, 0.000e+00,
       2.000e+00, 0.000e+00, 0.000e+00, 0.000e+00, 0.000e+00, 0.000e+00,
       0.000e+00, 0.000e+00, 3.000e+00, 0.000e+00, 0.000e+00, 0.000e+00,
       0.000e+00, 0.000e+00, 6.000e+00, 0.000e+00, 0.000e+00, 0.000e+00,
       0.000e+00, 0.000e+00, 0.000e+00, 0.000e+00, 1.500e+01, 0.000e+00,
       1.800e+01, 2.100e+01, 3.000e+00, 1.200e+01, 9.000e+00, 3.000e+01,
       4.800e+01, 2.400e+01, 4.800e+01, 8.400e+01, 8.700e+01, 1.390e+02,
       1.500e+02, 1.440e+02, 1.690e+02, 1.890e+02, 2.200e+02, 1.560e+02,
       1.620e+02, 1.880e+02, 2.430e+02, 2.040e+02, 1.830e+02, 1.980e+02,
       1.380e+02, 1.500e+02, 1.720e+02, 1.590e+02, 1.380e+02, 1.110e+02,
       8.700e+01, 4.500e+01, 1.470e+02, 1.020e+02, 9.300e+01, 1.320e+02,
       9.000e+01, 6.300e+01, 5.100e+01, 8.400e+01, 1.200e+02, 7.500e+01,
       1.140e+02, 6.900e+01, 7.500e+01, 6.300e+01, 1.320e+02, 1.600e+02,
       9.000e+01, 1.790e+02, 2.110e+02, 9.600e+01, 7.800e+01, 1.660e+02,
       8.700e+01, 9.000e+01, 8.700e+01, 1.000e+02, 6.900e+01, 5.200e+01,
       7.000e+01, 7.500e+01, 4.800e+01, 7.200e+01, 4.200e+01, 3.900e+01,
       3.900e+01, 3.600e+01, 3.000e+01, 7.000e+01, 2.400e+01, 5.400e+01,
       3.300e+01, 2.500e+01, 2.800e+01, 4.000e+01, 3.700e+01, 6.900e+01,
       2.400e+01, 2.400e+01, 2.600e+01, 2.100e+01, 1.500e+01, 2.700e+01,
       2.700e+01, 3.600e+01, 3.600e+01, 3.000e+01, 3.600e+01, 3.900e+01,
       9.000e+00, 2.700e+01, 1.500e+01, 3.100e+01, 2.400e+01, 1.800e+01,
       4.500e+01, 4.800e+01, 3.200e+01, 4.800e+01, 4.300e+01, 1.900e+01,
       2.400e+01, 6.500e+01, 1.600e+01, 1.800e+01, 2.100e+01, 2.400e+01,
       5.500e+01, 2.500e+01, 5.700e+01, 5.500e+01, 2.100e+01, 2.600e+01,
       3.300e+01, 2.400e+01, 2.400e+01, 4.300e+01, 3.100e+01, 2.500e+01,
       4.000e+01, 4.000e+01, 1.900e+01, 3.300e+01, 5.400e+01, 5.400e+01,
       6.000e+01, 8.400e+01, 6.600e+01, 5.400e+01, 3.700e+01, 6.600e+01,
       9.900e+01, 9.300e+01, 9.900e+01, 7.400e+01, 9.100e+01, 8.600e+01,
       9.000e+01, 1.010e+02, 9.600e+01, 7.900e+01, 6.400e+01, 6.300e+01,
       6.000e+01, 1.570e+02, 6.300e+01, 1.340e+02, 1.320e+02, 1.410e+02,
       1.240e+02, 1.030e+02, 1.060e+02, 1.200e+02, 1.650e+02, 1.260e+02,
       1.290e+02, 1.190e+02, 1.890e+02, 2.280e+02, 2.470e+02, 2.790e+02,
       3.320e+02, 1.830e+02, 2.380e+02, 1.530e+02, 2.200e+02, 2.550e+02,
       3.240e+02, 3.630e+02, 2.230e+02, 1.830e+02, 1.830e+02, 1.360e+02,
       2.940e+02, 3.340e+02, 3.150e+02, 3.150e+02, 2.590e+02, 1.770e+02,
       2.910e+02, 2.430e+02, 4.470e+02, 3.750e+02, 3.310e+02, 3.330e+02,
       2.730e+02, 3.420e+02, 3.660e+02, 3.630e+02, 4.470e+02, 3.030e+02,
       2.520e+02, 2.470e+02, 4.600e+02, 4.260e+02, 4.590e+02, 3.900e+02,
       3.240e+02, 3.780e+02, 2.940e+02, 2.400e+02, 4.450e+02, 2.400e+02,
       2.940e+02, 3.480e+02, 2.640e+02, 2.880e+02, 3.300e+02, 3.570e+02,
       4.320e+02, 2.850e+02, 3.750e+02, 3.700e+02, 2.700e+02, 3.360e+02,
       3.660e+02, 3.870e+02, 5.460e+02, 4.800e+02, 3.360e+02, 4.260e+02,
       4.710e+02, 2.730e+02, 4.560e+02, 6.330e+02, 4.080e+02, 6.120e+02,
       4.200e+02, 7.680e+02, 8.940e+02, 6.000e+02, 8.220e+02, 8.850e+02,
       6.870e+02, 7.770e+02, 7.950e+02, 8.140e+02, 1.020e+03, 1.149e+03,
       1.140e+03, 1.122e+03, 1.011e+03, 1.101e+03, 1.773e+03, 1.443e+03,
       1.347e+03, 1.707e+03, 1.507e+03, 1.491e+03, 2.097e+03, 1.980e+03,
       1.755e+03, 1.629e+03, 2.835e+03, 1.767e+03, 1.650e+03, 2.199e+03,
       2.580e+03, 1.926e+03, 2.802e+03, 2.169e+03, 1.963e+03, 2.049e+03,
       2.526e+03, 2.634e+03, 2.274e+03, 2.115e+03, 2.455e+03, 1.914e+03,
       1.633e+03, 2.322e+03, 2.271e+03, 2.064e+03, 1.763e+03, 2.332e+03,
       1.714e+03, 1.743e+03, 1.977e+03, 2.293e+03, 2.322e+03, 2.142e+03,
       2.206e+03, 1.695e+03, 1.547e+03, 1.759e+03, 1.903e+03, 1.866e+03,
       1.920e+03, 1.734e+03, 1.512e+03, 1.210e+03, 1.536e+03, 1.608e+03,
       1.662e+03, 1.215e+03, 9.930e+02, 1.428e+03, 1.212e+03, 1.854e+03,
       1.782e+03, 1.863e+03, 1.788e+03, 1.236e+03, 1.605e+03, 1.551e+03,
       1.698e+03, 1.816e+03, 1.637e+03, 1.582e+03, 1.589e+03, 1.347e+03,
       1.225e+03, 1.509e+03, 1.566e+03, 1.485e+03, 1.716e+03, 1.382e+03,
       9.500e+02, 1.137e+03, 1.534e+03, 1.674e+03, 1.560e+03, 1.454e+03,
       1.482e+03, 1.077e+03, 1.003e+03, 1.611e+03, 1.500e+03, 1.566e+03,
       1.259e+03, 1.434e+03, 8.820e+02, 1.192e+03, 1.206e+03, 1.416e+03,
       1.425e+03, 1.329e+03, 1.447e+03, 1.040e+03, 1.182e+03, 1.423e+03,
       1.348e+03, 1.299e+03, 1.354e+03, 1.344e+03, 1.080e+03, 9.670e+02,
       1.174e+03, 1.848e+03, 1.590e+03, 1.590e+03, 1.474e+03, 1.363e+03,
       1.497e+03, 1.557e+03, 1.518e+03, 1.779e+03, 1.550e+03, 1.479e+03,
       1.350e+03, 1.371e+03, 1.672e+03, 1.689e+03, 1.881e+03, 1.638e+03,
       1.670e+03, 1.197e+03, 1.491e+03, 1.585e+03, 1.740e+03, 1.977e+03,
       1.650e+03, 1.576e+03, 1.447e+03, 1.446e+03, 1.536e+03, 1.791e+03,
       2.205e+03, 1.731e+03, 1.838e+03, 1.980e+03, 1.752e+03, 2.161e+03,
       2.412e+03, 2.967e+03, 2.653e+03, 2.475e+03, 2.406e+03, 2.356e+03,
       3.067e+03, 2.641e+03, 3.220e+03, 3.315e+03, 3.082e+03, 2.736e+03,
       2.784e+03, 3.102e+03, 3.953e+03, 3.654e+03, 3.903e+03, 3.189e+03,
       3.117e+03, 2.819e+03, 3.012e+03, 3.672e+03, 2.787e+03, 3.255e+03,
       2.901e+03, 3.057e+03, 2.241e+03, 2.790e+03, 2.880e+03, 3.037e+03,
       2.673e+03, 2.620e+03, 2.298e+03, 2.325e+03, 2.391e+03, 2.611e+03,
       2.241e+03, 2.676e+03, 2.067e+03, 1.935e+03, 1.704e+03, 1.887e+03,
       2.133e+03, 1.951e+03, 1.846e+03, 1.872e+03, 1.707e+03, 1.548e+03,
       1.722e+03, 1.812e+03, 1.437e+03, 1.356e+03, 1.455e+03, 1.320e+03,
       1.191e+03, 1.569e+03, 1.123e+03, 1.299e+03, 1.059e+03, 1.003e+03,
       9.060e+02, 7.590e+02, 7.590e+02, 1.179e+03, 9.430e+02, 7.530e+02,
       7.360e+02, 6.540e+02, 4.620e+02, 5.940e+02, 5.970e+02, 5.350e+02,
       6.810e+02, 4.140e+02, 4.020e+02, 4.620e+02, 5.100e+02, 4.730e+02,
       4.320e+02, 3.240e+02, 3.540e+02, 2.130e+02, 3.180e+02, 3.270e+02,
       3.720e+02, 3.430e+02, 2.760e+02, 2.740e+02, 1.350e+02, 1.710e+02,
       2.490e+02, 2.220e+02, 2.200e+02, 1.620e+02, 1.570e+02, 1.180e+02,
       9.600e+01, 1.360e+02, 1.440e+02, 1.060e+02, 9.300e+01, 1.120e+02,
       7.400e+01, 1.340e+02, 1.350e+02, 1.980e+02, 1.530e+02, 1.530e+02,
       1.150e+02, 8.800e+01, 1.260e+02, 9.500e+01, 1.680e+02, 1.410e+02,
       1.620e+02, 1.920e+02, 1.440e+02, 1.950e+02, 2.340e+02, 2.910e+02,
       3.420e+02, 2.700e+02, 2.580e+02, 2.940e+02, 4.040e+02, 5.790e+02,
       6.110e+02, 7.030e+02, 4.650e+02, 5.880e+02, 5.550e+02, 6.030e+02])

help(np.convolve)

Help on function convolve in module numpy:

convolve(a, v, mode='full')
    Returns the discrete, linear convolution of two one-dimensional sequences.
    
    The convolution operator is often seen in signal processing, where it
    models the effect of a linear time-invariant system on a signal [1]_.  In
    probability theory, the sum of two independent random variables is
    distributed according to the convolution of their individual
    distributions.
    
    If `v` is longer than `a`, the arrays are swapped before computation.
    
    Parameters
    ----------
    a : (N,) array_like
        First one-dimensional input array.
    v : (M,) array_like
        Second one-dimensional input array.
    mode : {'full', 'valid', 'same'}, optional
        'full':
          By default, mode is 'full'.  This returns the convolution
          at each point of overlap, with an output shape of (N+M-1,). At
          the end-points of the convolution, the signals do not overlap
          completely, and boundary effects may be seen.
    
        'same':
          Mode 'same' returns output of length ``max(M, N)``.  Boundary
          effects are still visible.
    
        'valid':
          Mode 'valid' returns output of length
          ``max(M, N) - min(M, N) + 1``.  The convolution product is only given
          for points where the signals overlap completely.  Values outside
          the signal boundary have no effect.
    
    Returns
    -------
    out : ndarray
        Discrete, linear convolution of `a` and `v`.
    
    See Also
    --------
    scipy.signal.fftconvolve : Convolve two arrays using the Fast Fourier
                               Transform.
    scipy.linalg.toeplitz : Used to construct the convolution operator.
    polymul : Polynomial multiplication. Same output as convolve, but also
              accepts poly1d objects as input.
    
    Notes
    -----
    The discrete convolution operation is defined as
    
    .. math:: (a * v)[n] = \sum_{m = -\infty}^{\infty} a[m] v[n - m]
    
    It can be shown that a convolution :math:`x(t) * y(t)` in time/space
    is equivalent to the multiplication :math:`X(f) Y(f)` in the Fourier
    domain, after appropriate padding (padding is necessary to prevent
    circular convolution).  Since multiplication is more efficient (faster)
    than convolution, the function `scipy.signal.fftconvolve` exploits the
    FFT to calculate the convolution of large data-sets.
    
    References
    ----------
    .. [1] Wikipedia, "Convolution",
        https://en.wikipedia.org/wiki/Convolution
    
    Examples
    --------
    Note how the convolution operator flips the second array
    before "sliding" the two across one another:
    
    >>> np.convolve([1, 2, 3], [0, 1, 0.5])
    array([0. , 1. , 2.5, 4. , 1.5])
    
    Only return the middle values of the convolution.
    Contains boundary effects, where zeros are taken
    into account:
    
    >>> np.convolve([1,2,3],[0,1,0.5], 'same')
    array([1. ,  2.5,  4. ])
    
    The two arrays are of the same length, so there
    is only one position where they completely overlap:
    
    >>> np.convolve([1,2,3],[0,1,0.5], 'valid')
    array([2.5])

np.convolve([1, 2, 3], [0,1,0.5])

array([0. , 1. , 2.5, 4. , 1.5])

np.convolve([1,2,3],[0,1,0.5], 'same')

array([1. , 2.5, 4. ])

np.convolve([1,2,3],[0,1,0.5], 'valid')

array([2.5])

Plot

fig, ax = plt.subplots(1, 2, figsize=(14,5))
ax[0].plot(x,y, label='$y(x)$', c='orange', zorder = 0, alpha=0.9)
ax[0].plot(x[7:-7], y_smooth, label=r'$y_{{smooth}}(x)$', c='g', zorder=1, alpha = 0.8)
ax[1].plot(x,dydx, label='$y\'(x)$', color='gold')
ax[1].plot(x[7:-7],dysdx, label='$y_{smooth}\'(x)$', color='m')
ax[1].set_ylim(-100,120)
ax[1].grid()
[a.legend() for a in ax]
[a.set_xlabel('Time [Days]') for a in ax]
ax[0].set_ylabel('Cases per Day')
ax[1].set_ylabel('$\Delta$ (Cases per Day) / $\Delta t$')
fig.tight_layout()
plt.show()

Quasi-Symbolic Case

In this case you know your function $f(x) = ...$ but the function is not given by a typical expression. For example.

$$f(u) = \text{max} \left\{ \left|e^{-x_iu^2} - y_i\right| \right\} $$

for some array of $x_i$'s and $y_i$'s

x = np.linspace(0, 1, 500)
y = np.exp(-x*2.15**2) + 0.1*np.random.randn(len(x))

plt.scatter(x,y)
plt.xlabel('$x_i$', fontsize=20)
plt.ylabel('$y_i$', fontsize=20)
plt.show()

Define function

def f(u):
    return max(np.abs(np.exp(-x*u**2) -y))

Compute the function values

u = np.linspace(0,10,40)
f_u = np.vectorize(f)(u)

Plot

plt.plot(u, f_u, 'o--')
plt.xlabel('$u$', fontsize=20)
plt.ylabel('$f(u)$', fontsize=20)
plt.show()

You could always just take the derivative of the numerical array f_u, but there is a specific derivative function better assigned for this

Compute derivative

derivative(f, 0.2, dx=1e-6)

$\displaystyle -0.304508026238715$

dfdu = np.vectorize(derivative)(f, u, dx=1e-6)

Plot

plt.plot(u, dfdu)

[<matplotlib.lines.Line2D at 0x7f87c9684d60>]

help(np.vectorize)

Help on class vectorize in module numpy:

class vectorize(builtins.object)
 |  vectorize(pyfunc, otypes=None, doc=None, excluded=None, cache=False, signature=None)
 |  
 |  vectorize(pyfunc, otypes=None, doc=None, excluded=None, cache=False,
 |            signature=None)
 |  
 |  Generalized function class.
 |  
 |  Define a vectorized function which takes a nested sequence of objects or
 |  numpy arrays as inputs and returns a single numpy array or a tuple of numpy
 |  arrays. The vectorized function evaluates `pyfunc` over successive tuples
 |  of the input arrays like the python map function, except it uses the
 |  broadcasting rules of numpy.
 |  
 |  The data type of the output of `vectorized` is determined by calling
 |  the function with the first element of the input.  This can be avoided
 |  by specifying the `otypes` argument.
 |  
 |  Parameters
 |  ----------
 |  pyfunc : callable
 |      A python function or method.
 |  otypes : str or list of dtypes, optional
 |      The output data type. It must be specified as either a string of
 |      typecode characters or a list of data type specifiers. There should
 |      be one data type specifier for each output.
 |  doc : str, optional
 |      The docstring for the function. If None, the docstring will be the
 |      ``pyfunc.__doc__``.
 |  excluded : set, optional
 |      Set of strings or integers representing the positional or keyword
 |      arguments for which the function will not be vectorized.  These will be
 |      passed directly to `pyfunc` unmodified.
 |  
 |      .. versionadded:: 1.7.0
 |  
 |  cache : bool, optional
 |      If `True`, then cache the first function call that determines the number
 |      of outputs if `otypes` is not provided.
 |  
 |      .. versionadded:: 1.7.0
 |  
 |  signature : string, optional
 |      Generalized universal function signature, e.g., ``(m,n),(n)->(m)`` for
 |      vectorized matrix-vector multiplication. If provided, ``pyfunc`` will
 |      be called with (and expected to return) arrays with shapes given by the
 |      size of corresponding core dimensions. By default, ``pyfunc`` is
 |      assumed to take scalars as input and output.
 |  
 |      .. versionadded:: 1.12.0
 |  
 |  Returns
 |  -------
 |  vectorized : callable
 |      Vectorized function.
 |  
 |  See Also
 |  --------
 |  frompyfunc : Takes an arbitrary Python function and returns a ufunc
 |  
 |  Notes
 |  -----
 |  The `vectorize` function is provided primarily for convenience, not for
 |  performance. The implementation is essentially a for loop.
 |  
 |  If `otypes` is not specified, then a call to the function with the
 |  first argument will be used to determine the number of outputs.  The
 |  results of this call will be cached if `cache` is `True` to prevent
 |  calling the function twice.  However, to implement the cache, the
 |  original function must be wrapped which will slow down subsequent
 |  calls, so only do this if your function is expensive.
 |  
 |  The new keyword argument interface and `excluded` argument support
 |  further degrades performance.
 |  
 |  References
 |  ----------
 |  .. [1] :doc:`/reference/c-api/generalized-ufuncs`
 |  
 |  Examples
 |  --------
 |  >>> def myfunc(a, b):
 |  ...     "Return a-b if a>b, otherwise return a+b"
 |  ...     if a > b:
 |  ...         return a - b
 |  ...     else:
 |  ...         return a + b
 |  
 |  >>> vfunc = np.vectorize(myfunc)
 |  >>> vfunc([1, 2, 3, 4], 2)
 |  array([3, 4, 1, 2])
 |  
 |  The docstring is taken from the input function to `vectorize` unless it
 |  is specified:
 |  
 |  >>> vfunc.__doc__
 |  'Return a-b if a>b, otherwise return a+b'
 |  >>> vfunc = np.vectorize(myfunc, doc='Vectorized `myfunc`')
 |  >>> vfunc.__doc__
 |  'Vectorized `myfunc`'
 |  
 |  The output type is determined by evaluating the first element of the input,
 |  unless it is specified:
 |  
 |  >>> out = vfunc([1, 2, 3, 4], 2)
 |  >>> type(out[0])
 |  <class 'numpy.int64'>
 |  >>> vfunc = np.vectorize(myfunc, otypes=[float])
 |  >>> out = vfunc([1, 2, 3, 4], 2)
 |  >>> type(out[0])
 |  <class 'numpy.float64'>
 |  
 |  The `excluded` argument can be used to prevent vectorizing over certain
 |  arguments.  This can be useful for array-like arguments of a fixed length
 |  such as the coefficients for a polynomial as in `polyval`:
 |  
 |  >>> def mypolyval(p, x):
 |  ...     _p = list(p)
 |  ...     res = _p.pop(0)
 |  ...     while _p:
 |  ...         res = res*x + _p.pop(0)
 |  ...     return res
 |  >>> vpolyval = np.vectorize(mypolyval, excluded=['p'])
 |  >>> vpolyval(p=[1, 2, 3], x=[0, 1])
 |  array([3, 6])
 |  
 |  Positional arguments may also be excluded by specifying their position:
 |  
 |  >>> vpolyval.excluded.add(0)
 |  >>> vpolyval([1, 2, 3], x=[0, 1])
 |  array([3, 6])
 |  
 |  The `signature` argument allows for vectorizing functions that act on
 |  non-scalar arrays of fixed length. For example, you can use it for a
 |  vectorized calculation of Pearson correlation coefficient and its p-value:
 |  
 |  >>> import scipy.stats
 |  >>> pearsonr = np.vectorize(scipy.stats.pearsonr,
 |  ...                 signature='(n),(n)->(),()')
 |  >>> pearsonr([[0, 1, 2, 3]], [[1, 2, 3, 4], [4, 3, 2, 1]])
 |  (array([ 1., -1.]), array([ 0.,  0.]))
 |  
 |  Or for a vectorized convolution:
 |  
 |  >>> convolve = np.vectorize(np.convolve, signature='(n),(m)->(k)')
 |  >>> convolve(np.eye(4), [1, 2, 1])
 |  array([[1., 2., 1., 0., 0., 0.],
 |         [0., 1., 2., 1., 0., 0.],
 |         [0., 0., 1., 2., 1., 0.],
 |         [0., 0., 0., 1., 2., 1.]])
 |  
 |  Methods defined here:
 |  
 |  __call__(self, *args, **kwargs)
 |      Return arrays with the results of `pyfunc` broadcast (vectorized) over
 |      `args` and `kwargs` not in `excluded`.
 |  
 |  __init__(self, pyfunc, otypes=None, doc=None, excluded=None, cache=False, signature=None)
 |      Initialize self.  See help(type(self)) for accurate signature.
 |  
 |  ----------------------------------------------------------------------
 |  Data descriptors defined here:
 |  
 |  __dict__
 |      dictionary for instance variables (if defined)
 |  
 |  __weakref__
 |      list of weak references to the object (if defined)

import numpy as np
import time

num = 100000000

a = np.random.randn(num)
b = np.random.randn(num)

Vectorized

start_V = time.time()
c = np.dot(a, b)
end_V = time.time()
print(c)
print('Vectorized: ' + str((end - start) * 1000) + 'ms')

-11116.25556897978
Vectorized: 32024.862051010132ms

Non-vectorized

start_NV = time.time()
c = 0
for i in range(num):
    c += a[i] * b[i]
end_NV = time.time()
print(c)
print('Non-Vectorized: ' + str((end - start) * 1000) + 'ms')
print('Non-Vectorized/Vectorized: ' + str((end_NV - start_NV)/ (end_V - start_V)) + 'ms')

-11116.25556897954
Non-Vectorized: 32024.862051010132ms
Non-Vectorized/Vectorized: 387.99510348308934ms

The vectorized function evaluates pyfunc over successive tuples of the input arrays like the python map function, except it uses the broadcasting rules of numpy.

The data type of the output of vectorized is determined by calling the function with the first element of the input. This can be avoided by specifying the otypes argument.

def myfunc(a, b):

    "Return a-b if a>b, otherwise return a+b"

    if a > b:

        return a - b

    else:

        return a + b

vfunc = np.vectorize(myfunc)

vfunc([1, 2, 3, 4], 2)

array([3, 4, 1, 2])