# Introducción a la Computación Científica con Python

Versión original en inglés de J.R. Johansson (robert@riken.jp) http://dml.riken.jp/~rob/

Traducido/Adaptado por [G.F. Rubilar](http://google.com/+GuillermoRubilar).

La última versión de estos [notebooks de IPython](http://ipython.org/notebook.html) está disponible en [http://github.com/gfrubi/clases-python-cientifico](http://github.com/gfrubi/clases-python-cientifico).

La última versión del original (en inglés) de estos [notebooks de IPython](http://ipython.org/notebook.html) está disponible en [http://github.com/jrjohansson/scientific-python-lectures](http://github.com/jrjohansson/scientific-python-lectures).

Los otros notebooks de esta serie están listados en [http://jrjohansson.github.com](http://jrjohansson.github.com).


```python
%matplotlib inline
```

## Introducción

Existen dos sistemas de álgebra simbólica para Python:

* [SymPy](http://sympy.org/en/index.html) - Un módulo que puede ser usado en cualquier programa Python, o bien en una sesión de IPython, que incluye poderosas herramientas para cálculo simbólico. 
* [Sage](http://www.sagemath.org/) - Sage es un sistema completo y poderoso que intenta suministrar un sistema de código abierto que compita con Mathematica and Maple. Sage no es un módulo de Python, sino un ambiente de cálculo simbólico que usa Python como su lenguaje de programación.

Sage es más poderoso que SymPy en algunos aspectos, pero ambos ofreces una lista completa de funcionalidades de cálculo simbólico. La ventaja de SymPy es que es un módulo normal de Python y se integra muy bien en un notebook de IPython. 

En esta clase veremos cómo usar SymPy en un notebook de IPython.

Para comenzar a usar SymPy en un programa Python o en un notebook, importamos el módulo `sympy`:


```python
from sympy import *
```

Para que los resultados sean formateados en $\LaTeX$ podemos usar:


```python
init_printing(use_latex=True)
```

## Variables simbólicas

En SymPy podemos crear símbolos para las variables con las que deseamos trabajar. Podemos crear un nuevo símbolo usando la clase `Symbol`:


```python
x = Symbol('x')
```


```python
(pi + x)**2
```


```python
# forma alternativa de definir (varios) símbolos
a, b, c = symbols("a, b, c")
```


```python
type(a)
```




    sympy.core.symbol.Symbol



Podemos agregar algunas propiedades a los símbolos cuando son creados:


```python
x = Symbol('x', real=True)
```


```python
x.is_imaginary
```




    False




```python
x = Symbol('x', positive=True)
```


```python
x > 0
```

### Números complejos

La unidad imaginaria es denotada por `I` en Sympy. 


```python
1+1*I
```


```python
I**2
```


```python
(x * I + 1)**2
```

### Números racionales

Existen tres tipos distintos de números en SymPy: `Real`, `Rational`, `Integer`: 


```python
r1 = Rational(4,5)
r2 = Rational(5,4)
```


```python
r1
```


```python
r1+r2
```


```python
r1/r2
```

### Evaluación numérica

SymPy usa una librería para trabajar con números con precisión arbitraria, y tiene expresiones SymPy  predefinidas para varias constantes matemáticas, tales como: `pi`, `e` y `oo` para el infinito.

Para evaluar numéricamente una expresión podemos usar la función `evalf` (o bien `N`). Ésta usa un argumento `n` que especifíca el número de cifras significativas.


```python
pi.evalf(n=50)
```


```python
y = (x + pi)**2
```


```python
N(y, 5) # equivalente a evalf
```

Cuando evaluamos numéricamente expresiones a menudo deseamos substituir un símbolo por un valor numérico. En SymPy hacemos esto usando la función `subs`:


```python
y.subs(x, 1.5)
```


```python
N(y.subs(x, 1.5))
```

La función `subs` también puede ser usada para substituir símbolos y expresiones:


```python
y.subs(x, a+pi)
```

También podemos combinar la evaluación numérica de expresiones con arreglos NumPy:


```python
import numpy
```


```python
x_vec = numpy.arange(0, 10, 0.1)
```


```python
y_vec = numpy.array([N(((x + pi)**2).subs(x, xx)) for xx in x_vec])
```


```python
from matplotlib.pyplot import *
plot(x_vec, y_vec);
```

Sin embargo, este tipo de evaluación numérica puede ser muy lenta, y existe una forma mucho más eficiente de realizar la misma tarea: Usar la función `lambdify` para "mapear" una expresión de Sympy a una función que es mucho más eficiente para la evaluación numérica:


```python
f = lambdify([x], (x + pi)**2, 'numpy')  # el primer argumento es una lista de variables de las que la función f dependerá: en este caso sólo x -> f(x)
```


```python
type(f)
```




    function




```python
y_vec = f(x_vec)  # ahora podemos pasar directamente un arreglo Numpy. Así f(x) es evaluado más eficientemente
```

La mayor eficiencia de usar funciones "lambdificadas" en lugar de usar evalación numérica directa puede ser significativa, a menudo de varios órdenes de magnitud. Aún en este sencillo ejemplo obtenemos un aumento de velocidad importante:


```python
%%timeit

y_vec = numpy.array([N(((x + pi)**2).subs(x, xx)) for xx in x_vec])
```

    100 loops, best of 3: 19.1 ms per loop



```python
%%timeit

y_vec = f(x_vec)
```

    The slowest run took 18.93 times longer than the fastest. This could mean that an intermediate result is being cached 
    1000000 loops, best of 3: 1.49 µs per loop


## Manipulaciones algebráicas

Uno de los usos principales de un sistema de cálculo simbólico es que realiza manipulaciones algebráicas de expresiones. Por ejemplo, si queremos expandir un producto, factorizar una expresión, o simplificar un resultado. En esta sección presentamos las funciones para realizar estas operaciones básicas en SymPy.

### Expand and factor

Primeros pasos en la manipulación algebráica


```python
(x+1)*(x+2)*(x+3)
```


```python
expand((x+1)*(x+2)*(x+3))
```

La función `expand` acepta varias argumentos clave con los que se puede indicar qué tipo de expansión deseamos realizar. Por ejemplo, para expandir expresiones trigonométricas, usamos el argumento clave `trig=True`:


```python
sin(a+b)
```


```python
expand(sin(a+b), trig=True)
```

Ver `help(expand)` para una descripción detallada de los distintos tipos de expansiones que la función `expand` puede realizar.

También podemos factorizar expresiones, usando la función `factor` de SymPy: 


```python
factor(x**3 + 6 * x**2 + 11*x + 6)
```

### Simplify

La función `simplify` intenta simplificar una expresión usando distinta técnicas. Existen también alternativas más específicas a la función `simplify`: `trigsimp`, `powsimp`, `logcombine`, etc. 

El uso básico de estas funciones en el siguiente:


```python
# simplify expande un producto
simplify((x+1)*(x+2)*(x+3))
```


```python
# simplify usa identidades trigonometricas
simplify(sin(a)**2 + cos(a)**2)
```


```python
simplify(cos(x)/sin(x))
```

## apart and together

Podemos también manipular expresiones simbólicas que involucran fracciones usando las funciones `apart` y `together`. La primera de estas funciones separa una fracción en sus correspondientes fracciones parciales; la segunda hace todo lo contrario.


```python
f1 = 1/((a+1)*(a+2))
```


```python
f1
```


```python
apart(f1)
```


```python
f2 = 1/(a+2) + 1/(a+3)
```


```python
f2
```


```python
together(f2)
```

Usualmente `Simplify` combina las fracciones, pero no factoriza: 


```python
simplify(f2)
```

## Cálculo

Además de realizar manipulaciones algebráicas, SimPy puede realizar operaciones de cálculo, tales como derivar y derivar expresiones.

### Derivación

Derviar es usualmente algo simple. Usamos la función `diff`. El primer argumento es una expresión que será derivada, y el segundo argumento es el símbolo respecto al cual se realizará la derivada:


```python
y
```


```python
diff(y**2, x)
```

Para calcular derivadas de orden superior podemos usar:


```python
diff(y**2, x, x)
```

o bien


```python
diff(y**2, x, 2) # hace lo mismo
```

To calculate the derivative of a multivariate expression, we can do:


```python
x, y, z = symbols("x,y,z")
```


```python
f = sin(x*y) + cos(y*z)
```

$\frac{d^3f}{dxdy^2}$


```python
diff(f, x, 1, y, 2)
```

### Integration

Integration is done in a similar fashion:


```python
f
```


```python
integrate(f, x)
```

By providing limits for the integration variable we can evaluate definite integrals:


```python
integrate(f, (x, -1, 1))
```

and also improper integrals


```python
integrate(exp(-x**2), (x, -oo, oo))
```

Remember, `oo` is the SymPy notation for inifinity.

### Sums and products

We can evaluate sums and products using the functions: 'Sum'


```python
n = Symbol("n")
```


```python
Sum(1/n**2, (n, 1, 10))
```


```python
Sum(1/n**2, (n,1, 10)).evalf()
```


```python
Sum(1/n**2, (n, 1, oo)).evalf()
```

Products work much the same way:


```python
Product(n, (n, 1, 10)) # 10!
```

### Limits

Limits can be evaluated using the `limit` function. For example, 


```python
limit(sin(x)/x, x, 0)
```

We can use 'limit' to check the result of derivation using the `diff` function:


```python
f
```


```python
diff(f, x)
```

$\displaystyle \frac{\mathrm{d}f(x,y)}{\mathrm{d}x} = \frac{f(x+h,y)-f(x,y)}{h}$


```python
h = Symbol("h")
```


```python
limit((f.subs(x, x+h) - f)/h, h, 0)
```

OK!

We can change the direction from which we approach the limiting point using the `dir` keywork argument:


```python
limit(1/x, x, 0, dir="+")
```


```python
limit(1/x, x, 0, dir="-")
```

### Series

Series expansion is also one of the most useful features of a CAS. In SymPy we can perform a series expansion of an expression using the `series` function:


```python
series(exp(x), x)
```

By default it expands the expression around $x=0$, but we can expand around any value of $x$ by explicitly include a value in the function call:


```python
series(exp(x), x, 1)
```

And we can explicitly define to which order the series expansion should be carried out:


```python
series(exp(x), x, 1, 10)
```

The series expansion includes the order of the approximation, which is very useful for keeping track of the order of validity when we do calculations with series expansions of different order:


```python
s1 = cos(x).series(x, 0, 5)
s1
```


```python
s2 = sin(x).series(x, 0, 2)
s2
```


```python
expand(s1 * s2)
```

If we want to get rid of the order information we can use the `removeO` method:


```python
expand(s1.removeO() * s2.removeO())
```

But note that this is not the correct expansion of $\cos(x)\sin(x)$ to $5$th order:


```python
(cos(x)*sin(x)).series(x, 0, 6)
```

## Linear algebra

### Matrices

Matrices are defined using the `Matrix` class:


```python
m11, m12, m21, m22 = symbols("m11, m12, m21, m22")
b1, b2 = symbols("b1, b2")
```


```python
A = Matrix([[m11, m12],[m21, m22]])
A
```


```python
b = Matrix([[b1], [b2]])
b
```

With `Matrix` class instances we can do the usual matrix algebra operations:


```python
A**2
```


```python
A * b
```

And calculate determinants and inverses, and the like:


```python
A.det()
```


```python
A.inv()
```

## Solving equations

For solving equations and systems of equations we can use the `solve` function:


```python
solve(x**2 - 1, x)
```


```python
solve(x**4 - x**2 - 1, x)
```

System of equations:


```python
solve([x + y - 1, x - y - 1], [x,y])
```

In terms of other symbolic expressions:


```python
solve([x + y - a, x - y - c], [x,y])
```

## Quantum mechanics: noncommuting variables

How about non-commuting symbols? In quantum mechanics we need to work with noncommuting operators, and SymPy has a nice support for noncommuting symbols and even a subpackage for quantum mechanics related calculations!


```python
from sympy.physics.quantum import *
```

### States

We can define symbol states, kets and bras:


```python
Ket('psi')
```


```python
Bra('psi')
```


```python
u = Ket('0')
d = Ket('1')

a, b = symbols('alpha beta', complex=True)
```


```python
phi = a * u + sqrt(1-abs(a)**2) * d; phi
```


```python
Dagger(phi)
```


```python
Dagger(phi) * d
```

Use `qapply` to distribute a mutiplication:


```python
qapply(Dagger(phi) * d)
```


```python
qapply(Dagger(phi) * u)
```

### Operators


```python
A = Operator('A')
B = Operator('B')
```

Check if they are commuting!


```python
A * B == B * A
```


```python
expand((A+B)**3)
```


```python
c = Commutator(A,B)
c
```

We can use the `doit` method to evaluate the commutator:


```python
c.doit()
```

We can mix quantum operators with C-numbers:


```python
c = Commutator(a * A, b * B)
c
```

To expand the commutator, use the `expand` method with the `commutator=True` keyword argument:


```python
c = Commutator(A+B, A*B)
c.expand(commutator=True)
```


```python
Dagger(Commutator(A, B))
```


```python
ac = AntiCommutator(A,B)
```


```python
ac.doit()
```

#### Example: Quadrature commutator

Let's look at the commutator of the electromagnetic field quadatures $x$ and $p$. We can write the quadrature operators in terms of the creation and annihilation operators as:

$\displaystyle x = (a + a^\dagger)/\sqrt{2}$

$\displaystyle p = -i(a - a^\dagger)/\sqrt{2}$



```python
X = (A + Dagger(A))/sqrt(2)
X
```


```python
P = -I * (A - Dagger(A))/sqrt(2)
P
```

Let's expand the commutator $[x,p]$


```python
Commutator(X, P).expand(commutator=True).expand(commutator=True)
```

Here we see directly that the well known commutation relation for the quadratures

$[x,p]=i$

is a directly related to

$[A, A^\dagger]=1$ 

(which SymPy does not know about, and does not simplify).

For more details on the quantum module in SymPy, see:

* http://docs.sympy.org/0.7.2/modules/physics/quantum/index.html
* http://nbviewer.ipython.org/urls/raw.github.com/ipython/ipython/master/docs/examples/notebooks/sympy_quantum_computing.ipynb

## Further reading

* http://sympy.org/en/index.html - The SymPy projects web page.
* https://github.com/sympy/sympy - The source code of SymPy.
* http://live.sympy.org - Online version of SymPy for testing and demonstrations.


19fa548bbb19483169bdbe72626653598a988964
ipynb
Jupyter Notebook
Lecture-5-Sympy.ipynb
AnIsAsPe/clases-python-cientifico
a89e3b4db18bf867388d0f2f1f199544bc483f8e
CC-BY-3.0
2015-06-03T00:38:04.000Z
2022-03-15T21:01:52.000Z
Lecture-5-Sympy.ipynb
gfrubi/clases-python-cientifico
a89e3b4db18bf867388d0f2f1f199544bc483f8e
CC-BY-3.0
Lecture-5-Sympy.ipynb
gfrubi/clases-python-cientifico
a89e3b4db18bf867388d0f2f1f199544bc483f8e
CC-BY-3.0
2015-02-25T20:37:09.000Z
2021-01-22T17:27:27.000Z
Qwen/Qwen-72B
1. YES
2. YES