<jupyter_start><jupyter_text>Symbolic Python In standard mathematics we routinely write down abstract variables or concepts and manipulate them without ever assigning specific values to them. An example would be the quadratic equation$$  a x^2 + b x + c = 0 $$and its roots $x_{\pm}$: we can write down the solutions of the equation and discuss the existence, within the real numbers, of the roots, without specifying the particular values of the parameters $a, b$ and $c$.In a standard computer programming language, we can write *functions* that encapsulate the solutions of the equation, but calling those functions requires us to specify values of the parameters. In general, the value of a variable must be given before the variable can be used.However, there *do* exist *Computer Algebra Systems* that can perform manipulations in the "standard" mathematical form. Through the university you will have access to Wolfram Mathematica and Maple, which are commercial packages providing a huge range of mathematical tools. There are also freely available packages, such as SageMath and `sympy`. These are not always easy to use, as all CAS have their own formal languages that rarely perfectly match your expectations.Here we will briefly look at `sympy`, which is a pure Python CAS. `sympy` is not suitable for complex calculations, as it's far slower than the alternatives. However, it does interface very cleanly with Python, so can be used inside Python code, especially to avoid entering lengthy expressions.  sympy  Setting up Setting up `sympy` is straightforward:<jupyter_code>import sympy
sympy.init_printing()<jupyter_output><empty_output><jupyter_text>The standard `import` command is used. The `init_printing` command looks at your system to find the clearest way of displaying the output; this isn't necessary, but is helpful for understanding the results. To do *anything* in `sympy` we have to explicitly tell it if something is a variable, and what name it has. There are two commands that do this. To declare a single variable, use<jupyter_code>x = sympy.Symbol('x')<jupyter_output><empty_output><jupyter_text>To declare multiple variables at once, use<jupyter_code>y, z0 = sympy.symbols(('y', 'z_0'))<jupyter_output><empty_output><jupyter_text>Note that the "name" of the variable does not need to match the symbol with which it is displayed. We have used this with `z0` above:<jupyter_code>z0<jupyter_output><empty_output><jupyter_text>Once we have variables, we can define new variables by operating on old ones:<jupyter_code>a = x + y
b = y * z0
print("a={}. b={}.".format(a, b))
a<jupyter_output><empty_output><jupyter_text>In addition to variables, we can also define general functions. There is only one option for this:<jupyter_code>f = sympy.Function('f')<jupyter_output><empty_output><jupyter_text>In-built functions We have seen already that mathematical functions can be found in different packages. For example, the $\sin$ function appears in `math` as `math.sin`, acting on a single number. It also appears in `numpy` as `numpy.sin`, where it can act on vectors and arrays in one go. `sympy` re-implements many mathematical functions, for example as `sympy.sin`, which can act on abstract (`sympy`) variables.Whenever using `sympy` we should use `sympy` functions, as these can be manipulated and simplified. For example:<jupyter_code>c = sympy.sin(x)**2 + sympy.cos(x)**2
c
c.simplify()<jupyter_output><empty_output><jupyter_text>Note the steps taken here. `c` is an object, something that `sympy` has created. Once created it can be manipulated and simplified, using the methods on the object. It is useful to use tab completion to look at the available commands. For example,<jupyter_code>d = sympy.cosh(x)**2 - sympy.sinh(x)**2<jupyter_output><empty_output><jupyter_text>Now type `d.` and then tab, to inspect all the available methods. As before, we could do<jupyter_code>d.simplify()<jupyter_output><empty_output><jupyter_text>but there are many other options.  Solving equations Let us go back to our quadratic equation and check the solution. To define an *equation* we use the `sympy.Eq` function:<jupyter_code>a, b, c, x = sympy.symbols(('a', 'b', 'c', 'x'))
quadratic_equation = sympy.Eq(a*x**2+b*x+c, 0)
sympy.solve(quadratic_equation)<jupyter_output><empty_output><jupyter_text>What happened here? `sympy` is not smart enough to know that we wanted to solve for `x`! Instead, it solved for the first variable it encountered. Let us try again:<jupyter_code>sympy.solve(quadratic_equation, x)<jupyter_output><empty_output><jupyter_text>This is our expectation: multiple solutions, returned as a list. We can access and manipulate these results:<jupyter_code>roots = sympy.solve(quadratic_equation, x)
xplus, xminus = sympy.symbols(('x_{+}', 'x_{-}'))
xplus = roots[0]
xminus = roots[1]<jupyter_output><empty_output><jupyter_text>We can substitute in specific values for the parameters to find solutions:<jupyter_code>xplus_solution = xplus.subs([(a,1), (b,2), (c,3)])
xplus_solution<jupyter_output><empty_output><jupyter_text>We have a list of substitutions. Each substitution is given by a tuple, containing the variable to be replaced, and the expression replacing it. We do not have to substitute in numbers, as here, but could use other variables:<jupyter_code>xminus_solution = xminus.subs([(b,a), (c,a+z0)])
xminus_solution
xminus_solution.simplify()<jupyter_output><empty_output><jupyter_text>We can use similar syntax to solve *systems* of equations, such as$$ \begin{aligned}  x + 2 y &= 0, \\ xy & = z_0. \end{aligned} $$<jupyter_code>eq1 = sympy.Eq(x+2*y, 0)
eq2 = sympy.Eq(x*y, z0)
sympy.solve([eq1, eq2], [x, y])<jupyter_output><empty_output><jupyter_text>Differentiation and integration  Differentiation There is a standard function for differentiation, `diff`:<jupyter_code>expression = x**2*sympy.sin(sympy.log(x))
sympy.diff(expression, x)<jupyter_output><empty_output><jupyter_text>A parameter can control how many times to differentiate:<jupyter_code>sympy.diff(expression, x, 3)<jupyter_output><empty_output><jupyter_text>Partial differentiation with respect to multiple variables can also be performed by increasing the number of arguments:<jupyter_code>expression2 = x*sympy.cos(y**2 + x)
sympy.diff(expression2, x, 2, y, 3)<jupyter_output><empty_output><jupyter_text>There is also a function representing an *unevaluated* derivative:<jupyter_code>sympy.Derivative(expression2, x, 2, y, 3)<jupyter_output><empty_output><jupyter_text>These can be useful for display, building up a calculation in stages, simplification, or when the derivative cannot be evaluated. It can be explicitly evaluated using the `doit` function:<jupyter_code>sympy.Derivative(expression2, x, 2, y, 3).doit()<jupyter_output><empty_output><jupyter_text>Integration Integration uses the `integrate` function. This can calculate either definite or indefinite integrals, but will *not* include the integration constant.<jupyter_code>integrand=sympy.log(x)**2
sympy.integrate(integrand, x)
sympy.integrate(integrand, (x, 1, 10))<jupyter_output><empty_output><jupyter_text>The definite integral is specified by passing a tuple, with the variable to be integrated (here `x`) and the lower and upper limits (which can be expressions).Note that `sympy` includes an "infinity" object `oo` (two `o`'s), which can be used in the limits of integration:<jupyter_code>sympy.integrate(sympy.exp(-x), (x, 0, sympy.oo))<jupyter_output><empty_output><jupyter_text>Multiple integration for higher dimensional integrals can be performed:<jupyter_code>sympy.integrate(sympy.exp(-(x+y))*sympy.cos(x)*sympy.sin(y), x, y)
sympy.integrate(sympy.exp(-(x+y))*sympy.cos(x)*sympy.sin(y), 
                (x, 0, sympy.pi), (y, 0, sympy.pi))<jupyter_output><empty_output><jupyter_text>Again, there is an unevaluated integral:<jupyter_code>sympy.Integral(integrand, x)
sympy.Integral(integrand, (x, 1, 10))<jupyter_output><empty_output><jupyter_text>Again, the `doit` method will explicitly evaluate the result where possible.  Differential equations Defining and solving differential equations uses the pattern from the previous sections. We'll use the same example problem as in the `scipy` case, $$  \frac{\text{d} y}{\text{d} t} = e^{-t} - y, \qquad y(0) = 1. $$First we define that $y$ is a function, currently unknown, and $t$ is a variable.<jupyter_code>y = sympy.Function('y')
t = sympy.Symbol('t')<jupyter_output><empty_output><jupyter_text>`y` is a general function, and can be a function of anything at this point (any number of variables with any name). To use it consistently, we *must* refer to it explicitly as a function of $t$ everywhere. For example,<jupyter_code>y(t)<jupyter_output><empty_output><jupyter_text>We then define the differential equation. `sympy.Eq` defines the equation, and `diff` differentiates:<jupyter_code>ode = sympy.Eq(y(t).diff(t), sympy.exp(-t) - y(t))
ode<jupyter_output><empty_output><jupyter_text>Here we have used `diff` as a method applied to the function. As `sympy` can't differentiate $y(t)$ (as it doesn't have an explicit value), it leaves it unevaluated.We can now use the `dsolve` function to get the solution to the ODE. The syntax is very similar to the `solve` function used above:<jupyter_code>sympy.dsolve(ode, y(t))<jupyter_output><empty_output><jupyter_text>This is simple enough to solve, but we'll use symbolic methods to find the constant, by setting $t = 0$ and $y(t) = y(0) = 1$.<jupyter_code>general_solution = sympy.dsolve(ode, y(t))
value = general_solution.subs([(t,0), (y(0), 1)])
value<jupyter_output><empty_output><jupyter_text>We then find the specific solution of the ODE.<jupyter_code>ode_solution = general_solution.subs([(value.rhs,value.lhs)])
ode_solution<jupyter_output><empty_output><jupyter_text>Plotting `sympy` provides an interface to `matplotlib` so that expressions can be directly plotted. For example,<jupyter_code>%matplotlib inline
from matplotlib import rcParams
rcParams['figure.figsize']=(12,9)
sympy.plot(sympy.sin(x));<jupyter_output><empty_output><jupyter_text>We can explicitly set limits, for example<jupyter_code>sympy.plot(sympy.exp(-x)*sympy.sin(x**2), (x, 0, 1));<jupyter_output><empty_output><jupyter_text>We can plot the solution to the differential equation computed above:<jupyter_code>sympy.plot(ode_solution.rhs, xlim=(0, 1), ylim=(0.7, 1.05));<jupyter_output><empty_output><jupyter_text>This can be *visually* compared to the previous result. However, we would often like a more precise comparison, which requires numerically evaluating the solution to the ODE at specific points.  lambdify At the end of a symbolic calculation using `sympy` we will have a result that is often long and complex, and that is needed in another part of another code. We could type the appropriate expression in by hand, but this is tedious and error prone. A better way is to make the computer do it.The example we use here is the solution to the ODE above. We have solved it symbolically, and the result is straightforward. We can also solve it numerically using `scipy`. We want to compare the two.First, let us compute the `scipy` numerical result:<jupyter_code>from numpy import exp
from scipy.integrate import odeint
import numpy

def dydt(y, t):
    """
    Defining the ODE dy/dt = e^{-t} - y.
    
    Parameters
    ----------
    
    y : real
        The value of y at time t (the current numerical approximation)
    t : real
        The current time t
        
    Returns
    -------
    
    dydt : real
        The RHS function defining the ODE.
    """
    
    return exp(-t) - y

t_scipy = numpy.linspace(0.0, 1.0)
y0 = [1.0]

y_scipy = odeint(dydt, y0, t_scipy)<jupyter_output><empty_output><jupyter_text>We want to evaluate our `sympy` solution at the same points as our `scipy` solution, in order to do a direct comparison. In order to do that, we want to construct a function that computes our `sympy` solution, without typing it in. That is what `lambdify` is for: it creates a function from a sympy expression.First let us get the expression explicitly:<jupyter_code>ode_expression = ode_solution.rhs
ode_expression<jupyter_output><empty_output><jupyter_text>Then we construct the function using `lambdify`:<jupyter_code>from sympy.utilities.lambdify import lambdify

ode_function = lambdify((t,), ode_expression, modules='numpy')<jupyter_output><empty_output><jupyter_text>The first argument to `lambdify` is a tuple containing the arguments of the function to be created. In this case that's just `t`, the time(s) at which we want to evaluate the expression. The second argument to `lambdify` is the expression that we want converted into a function. The third argument, which is optional, tells `lambdify` that where possible it should use `numpy` functions. This means that we call the function using `numpy` arrays, it will calculate using `numpy` array expressions, doing the whole calculation in a single call.We now have a function that we can directly call:<jupyter_code>print("sympy solution at t=0: {}".format(ode_function(0.0)))
print("sympy solution at t=0.5: {}".format(ode_function(0.5)))<jupyter_output>sympy solution at t=0: 1.0
sympy solution at t=0.5: 0.9097959895689501<jupyter_text>And we can directly apply this function to the times at which the `scipy` solution is constructed, for comparison:<jupyter_code>y_sympy = ode_function(t_scipy)<jupyter_output><empty_output><jupyter_text>Now we can use `matplotlib` to plot both on the same figure:<jupyter_code>from matplotlib import pyplot
pyplot.plot(t_scipy, y_scipy[:,0], 'b-', label='scipy')
pyplot.plot(t_scipy, y_sympy, 'k--', label='sympy')
pyplot.xlabel(r'$t$')
pyplot.ylabel(r'$y$')
pyplot.legend(loc='upper right')
pyplot.show()<jupyter_output><empty_output><jupyter_text>We see good visual agreement everywhere. But how accurate is it?Now that we have `numpy` arrays explicitly containing the solutions, we can manipulate these to see the differences between solutions:<jupyter_code>pyplot.semilogy(t_scipy, numpy.abs(y_scipy[:,0]-y_sympy))
pyplot.xlabel(r'$t$')
pyplot.ylabel('Difference in solutions');<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10048794286118379, 'fraction_numerical': 0.007920231949649954, 'mean_word_length': 6.0009900990099005, 'pattern_counts': {'":': 0, '<': 160, '<?xml version=': 0, '>': 160, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '40', 'dir': 'jupyter-structured-clean-dedup', 'id': '529489', 'license': 'MIT', 'n_tokens_mistral': 4567, 'n_tokens_neox': 4425, 'n_words': 1608, 'path': '07-sympy.ipynb', 'repo_name': 'LaGuer/maths-with-python', 'size': '14141'}