import sympy
sympy.init_printing()

"""
Explanation: 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:
End of explanation
"""


x = sympy.Symbol('x')

"""
Explanation: 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
End of explanation
"""


y, z0 = sympy.symbols(('y', 'z_0'))

"""
Explanation: To declare multiple variables at once, use
End of explanation
"""


z0

"""
Explanation: 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:
End of explanation
"""


a = x + y
b = y * z0
print("a={}. b={}.".format(a, b))

a

"""
Explanation: Once we have variables, we can define new variables by operating on old ones:
End of explanation
"""


f = sympy.Function('f')

"""
Explanation: In addition to variables, we can also define general functions. There is only one option for this:
End of explanation
"""


c = sympy.sin(x)**2 + sympy.cos(x)**2

c

c.simplify()

"""
Explanation: 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:
End of explanation
"""


d = sympy.cosh(x)**2 - sympy.sinh(x)**2

"""
Explanation: 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,
End of explanation
"""


d.simplify()

"""
Explanation: Now type d. and then tab, to inspect all the available methods. As before, we could do
End of explanation
"""


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)

"""
Explanation: 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:
End of explanation
"""


sympy.solve(quadratic_equation, x)

"""
Explanation: 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:
End of explanation
"""


roots = sympy.solve(quadratic_equation, x)
xplus, xminus = sympy.symbols(('x_{+}', 'x_{-}'))
xplus = roots[0]
xminus = roots[1]

"""
Explanation: This is our expectation: multiple solutions, returned as a list. We can access and manipulate these results:
End of explanation
"""


xplus_solution = xplus.subs([(a,1), (b,2), (c,3)])
xplus_solution

"""
Explanation: We can substitute in specific values for the parameters to find solutions:
End of explanation
"""


xminus_solution = xminus.subs([(b,a), (c,a+z0)])
xminus_solution

xminus_solution.simplify()

"""
Explanation: 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:
End of explanation
"""


eq1 = sympy.Eq(x+2*y, 0)
eq2 = sympy.Eq(x*y, z0)
sympy.solve([eq1, eq2], [x, y])

"""
Explanation: We can use similar syntax to solve systems of equations, such as
$$ \begin{aligned}  x + 2 y &= 0, \ xy & = z_0. \end{aligned} $$
End of explanation
"""


expression = x**2*sympy.sin(sympy.log(x))
sympy.diff(expression, x)

"""
Explanation: Differentiation and integration
Differentiation
There is a standard function for differentiation, diff:
End of explanation
"""


sympy.diff(expression, x, 3)

"""
Explanation: A parameter can control how many times to differentiate:
End of explanation
"""


expression2 = x*sympy.cos(y**2 + x)
sympy.diff(expression2, x, 2, y, 3)

"""
Explanation: Partial differentiation with respect to multiple variables can also be performed by increasing the number of arguments:
End of explanation
"""


sympy.Derivative(expression2, x, 2, y, 3)

"""
Explanation: There is also a function representing an unevaluated derivative:
End of explanation
"""


sympy.Derivative(expression2, x, 2, y, 3).doit()

"""
Explanation: 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:
End of explanation
"""


integrand=sympy.log(x)**2
sympy.integrate(integrand, x)

sympy.integrate(integrand, (x, 1, 10))

"""
Explanation: Integration
Integration uses the integrate function. This can calculate either definite or indefinite integrals, but will not include the integration constant.
End of explanation
"""


sympy.integrate(sympy.exp(-x), (x, 0, sympy.oo))

"""
Explanation: 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:
End of explanation
"""


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))

"""
Explanation: Multiple integration for higher dimensional integrals can be performed:
End of explanation
"""


sympy.Integral(integrand, x)

sympy.Integral(integrand, (x, 1, 10))

"""
Explanation: Again, there is an unevaluated integral:
End of explanation
"""


y = sympy.Function('y')
t = sympy.Symbol('t')

"""
Explanation: 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.
End of explanation
"""


y(t)

"""
Explanation: 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,
End of explanation
"""


ode = sympy.Eq(y(t).diff(t), sympy.exp(-t) - y(t))
ode

"""
Explanation: We then define the differential equation. sympy.Eq defines the equation, and diff differentiates:
End of explanation
"""


sympy.dsolve(ode, y(t))

"""
Explanation: 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:
End of explanation
"""


general_solution = sympy.dsolve(ode, y(t))
value = general_solution.subs([(t,0), (y(0), 1)])
value

"""
Explanation: 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$.
End of explanation
"""


ode_solution = general_solution.subs([(value.rhs,value.lhs)])
ode_solution

"""
Explanation: We then find the specific solution of the ODE.
End of explanation
"""


%matplotlib inline
from matplotlib import rcParams
rcParams['figure.figsize']=(12,9)

sympy.plot(sympy.sin(x));

"""
Explanation: Plotting
sympy provides an interface to matplotlib so that expressions can be directly plotted. For example,
End of explanation
"""


sympy.plot(sympy.exp(-x)*sympy.sin(x**2), (x, 0, 1));

"""
Explanation: We can explicitly set limits, for example
End of explanation
"""


sympy.plot(ode_solution.rhs, xlim=(0, 1), ylim=(0.7, 1.05));

"""
Explanation: We can plot the solution to the differential equation computed above:
End of explanation
"""


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)

"""
Explanation: 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:
End of explanation
"""


ode_expression = ode_solution.rhs
ode_expression

"""
Explanation: 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:
End of explanation
"""


from sympy.utilities.lambdify import lambdify

ode_function = lambdify((t,), ode_expression, modules='numpy')

"""
Explanation: Then we construct the function using lambdify:
End of explanation
"""


print("sympy solution at t=0: {}".format(ode_function(0.0)))
print("sympy solution at t=0.5: {}".format(ode_function(0.5)))

"""
Explanation: 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:
End of explanation
"""


y_sympy = ode_function(t_scipy)

"""
Explanation: And we can directly apply this function to the times at which the scipy solution is constructed, for comparison:
End of explanation
"""


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()

"""
Explanation: Now we can use matplotlib to plot both on the same figure:
End of explanation
"""


pyplot.semilogy(t_scipy, numpy.abs(y_scipy[:,0]-y_sympy))
pyplot.xlabel(r'$t$')
pyplot.ylabel('Difference in solutions');

"""
Explanation: 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:
End of explanation
"""


