<jupyter_start><jupyter_text>¡ Haciendo mates con Python !  ** <NAME> ** y ** <NAME> ** para [PyConES16](http://2016.es.pycon.org/) - Almería -- ** [15 min] **  0. Índice> 1. Introducción.> 2. Mi tabla de multiplicar. > 3. Cambio de unidades.. > 4. Resolviendo ecuaciones con Python. > 5. Caso práctico 1: Resolviendo la ecuación del calor con Pyomo. > 6. Caso práctico 2: Resolviendo la ecuación de ondas con Pyomo.   1. Introducción Las ** matemáticas ** en si mismas se han convertido en la herrameinta fundamental en cualquier área de conocimiento relacionada con las ciencias puras, ingenierías, ciencias sociales... Pero también lo son en la vida cotidiana. Participamos de ellas cuando compramos una barra de pan o cuando planificamos eventos u otras actividades.  Por este motivo es interesante pensar que podamos utilizar lenguajes de programación para resolver problemas de una manera sencilla y sistemática, como veremos a continuación.  Cuando hablamos de matemáticas es inevitable pensar en números. Números que dependiendo del programa que estemos desarrollando deberán tener unas características u otras. Python distingue principalmente entre números enteros y números de punto flotante. Los primeros no contienen parte decimal, si bien si los segundos. Como se puede observar, Python los asocia con diferente tipología:<jupyter_code>print 'Esto es un núemro entero:'
print 3
type(3)
print 'Esto es un número de punto flotante:'
print 3.0
type(3.0)<jupyter_output>Esto es un número de punto flotante:
3.0<jupyter_text>Python permite el cambio automático entre el tipo de número:<jupyter_code>a = 2.0 
print ' Valor de a = ' , a 
b =  int(a)
print ' Valor de a entero = ' , b<jupyter_output>Valor de a =  2.0
 Valor de a entero =  2<jupyter_text>Y viceversa:<jupyter_code>print ' Valor de b = ' , b 
a = float(b) 
print ' Valor de b flotante = ' , a<jupyter_output>Valor de b =  2
 Valor de b flotante =  2.0<jupyter_text>Hay que llevar cuidado pues estos comandos pueden llevar a confusión. Solo se está cambiando la tipología del número, pero no se está haciendo un redondeo cuando utilizamos el comanto int. Veámoslo.<jupyter_code>try:
    a = int ( input ( 'Introduce un número:  '))
except NameError:
    print( 'Has introducido un carácter no válido ')

print ' El número que has escogido es ' , a<jupyter_output>Introduce un número:  1/10
 El número que has escogido es  0<jupyter_text>Pero las matemáticas también son operaciones. Entre las operaciones básicas se encuentran la suma, sustracción, producto y división.<jupyter_code>print ' Esto es una suma: 2 + 2 = ' , 2 + 2 
print ' Esto es una resta: 235 - 23.0 = ' , 235 - 23.0 
 print ' Esto es una multiplicación: 2.0 * 3 = ' , 2.0 * 3  
print ' Esto es una división: 100 / 5 = ' , 100 / 5<jupyter_output>Esto es una división: 100 / 5 =  20<jupyter_text>Python permite también trabajar con fracciones:<jupyter_code>from fractions import Fraction 
fraccion = Fraction ( 3 , 4)
print fraccion 
print fraccion , '+' , 1 , '=' , fraccion + 1<jupyter_output>3/4 + 1 = 7/4<jupyter_text>Pero si sumamos un número flotante obtendremos un número también flotante:<jupyter_code>print fraccion , '+', 1.0 , '=', fraccion + 1.0<jupyter_output>3/4 + 1.0 = 1.75<jupyter_text>Los números complejos también tienen lugar dentro de Python:<jupyter_code>a = 2 + 3j 
print a , type (a)<jupyter_output>(2+3j) <type 'complex'><jupyter_text>Llegados a este punto ya podemos comenzar a elaborar pequeñas funciones que nos sean de utilidad. Por ejemplo, si queremos saber si un número es factor de otro basta con dividir uno entre otro y comprobar si el resto es 0:<jupyter_code>def factor ( a , b) : 
    
    if b % a == 0 : 
        return True
    else:
        return False 
factor( 4 , 1024)
factor( 3 , 1000 )<jupyter_output><empty_output><jupyter_text>Incluso puede que estemos interesados en conocer cuáles son los factores de un número cualquiera:<jupyter_code>def factor (b):
    
    for i in range ( 1 , b + 1):
        
        if b % i == 0: 
            
            print(i) 
factor ( 100 )<jupyter_output>1
2
4
5
10
20
25
50
100<jupyter_text>Nótese que en este caso solo podemos introducir números enteros y positivos si queremos que el programa funcione correctamente:<jupyter_code>factor ( 100.0 )<jupyter_output><empty_output><jupyter_text>Para solucionarlo podemos adaptar nuestro código:<jupyter_code>def factor_v2 (b) :
    
    for i in range ( 1 , b + 1 ) : 
        
        if b % i == 0: 
            
            print(i) 
        
if __name__ == '__main__': 
    b = input (' Introduce un número para conocer sus factores: ')
    b = float (b) 
    
    if b > 0 and b.is_integer(): 
        factor( int(b) )
        
    else:
        print 'Introduce un entero positivo, por favor '<jupyter_output>Introduce un número para conocer sus factores: 100.0
1
2
4
5
10
20
25
50
100<jupyter_text>2. Mi tabla de multiplicar.  *Ejemplo extraído del libro: "Doing Math with Python - <NAME>* Podemos elaborar programas más complejos que nos ayuden, por ejemplo, a reproducir automáticamente cualquier tabla de multiplicar:<jupyter_code>def multi_table( a , b ):
    
    for i in range( 1 , b +1 ):
        
        c = a * i 
        
        print c 
multi_table( 1 , 10 )
def multi_table_v2( a , b ): 
    
    for i in range( 1 , b + 1 ):
        
        c = a * i 
        
        print '{0} x {1} = {2}' . format( a , i , c)
multi_table_v2(1 , 10)<jupyter_output>1 x 1 = 1
1 x 2 = 2
1 x 3 = 3
1 x 4 = 4
1 x 5 = 5
1 x 6 = 6
1 x 7 = 7
1 x 8 = 8
1 x 9 = 9
1 x 10 = 10<jupyter_text>3. Cambio de unidades.  Del mismo modo podemos elaborar un programa que nos ralice un cambio de unidades, como por ejemplo, de grados Celsius a grados Kelvin.<jupyter_code>def menu():
    print ' 1. De Celsius a Kelvin '
    print ' 2. De Kelvin a Celsius '
    
def Celsius_a_Kelvin(): 
    
    C = float( input( ' Introduce la temperatura en Celsius: ')) 
    K = C + 273.15 
    
    print ' Temperatura en grados kelvin: {0:2f}'.format(K) 
    
def Kelvin_a_Celsius(): 
    
    K = float( input( ' Introduce la temperatura en Kelvin: '))
    C = K - 273.15 
    
    print ' Temperatura en grados Celsius: {0:2f}'.format(C)
    
if __name__ == '__main__' :
    
    menu()
    
    Escoja = input( ' ¿Qué quieres hacer?: ')
    
    if Escoja == 1 : 
        
        Celsius_a_Kelvin() 
        
    if Escoja == 2 : 
        
        Kelvin_a_Celsius()<jupyter_output>1. De Celsius a Kelvin 
 2. De Kelvin a Celsius 
 ¿Qué quieres hacer?: 1
 Introduce la temperatura en Celsius: 0
 Temperatura en grados kelvin: 273.150000<jupyter_text>4. Resolviendo ecuaciones con Python  Existe una variedad de maneras para resolver ecuaciones en Python. Una de ellas es implementar un programa que relice todas las operaciones, como si resolvieramos la ecuación a mano:<jupyter_code>def raiz( a , b , c):
    
    D = ( b ** 2 - 4 * a * c ) ** 0.5 
    
    x1 = ( - b + D ) / ( 2 * a ) 
    x2 = ( - b - D ) / ( 2 * a ) 
    
    print 'raíz 1 --> x_1 = {0:.2f}'.format(x1) 
    print 'raíz 2 --Z x_2 = {0:.2f}'.format(x2) 

    return x1 , x2 

if __name__ == '__main__' : 
    
    print ' ax^" + bx + c = 0 '
    
    a = input('Introduce el valor de a:  ')
    b = input('Introduce el valor de b:  ')
    c = input('Introduce el valor de c:  ') 
    
    [x1 , x2] = raiz(a , b , c)
def raiz( a , b , c):
    
    D = ( b ** 2 - 4 * a * c ) ** 0.5 
    
    x1 = ( - b + D ) / ( 2 * a ) 
    x2 = ( - b - D ) / ( 2 * a ) 
    
    print ('raíz 1 --> x_1 = {0:.2f}').format(x1) 
    print ('raíz 2 --Z x_2 = {0:.2f}').format(x2) 

    return x1 , x2 

if __name__ == '__main__' : 
    
    print (' ax^" + bx + c = 0 ')
    
    a = input('Introduce el valor de a:  ')
    b = input('Introduce el valor de b:  ')
    c = input('Introduce el valor de c:  ') 
    
    try: 
        [x1 , x2] = raiz(a , b , c)
        
    except ValueError:
        print (' La ecuación tiene raíces complejas, no la puedo resolver')<jupyter_output>ax^" + bx + c = 0 
Introduce el valor de a:  1
Introduce el valor de b:  1
Introduce el valor de c:  1
 La ecuación tiene raíces complejas, no la puedo resolver<jupyter_text>Otra opción para resolver ecuaciones es utilizando el cálculo simbólico con SymPy:<jupyter_code>from sympy import Symbol , solve 
a = Symbol( 'a' ) 
type(a)
expresion = a - 5 
solve(expresion)<jupyter_output><empty_output><jupyter_text>Podemos resolver también ecuaciones cuadráticas, como la vista al comienzo del apartado:<jupyter_code>x = Symbol ( 'x' )
ecuacion = x**2 + 5*x + 4 
solve(ecuacion)<jupyter_output><empty_output><jupyter_text>En este caso si podemos obtener la solución en números complejos:<jupyter_code>ecuacion = x**2 + x + 1 
solve(ecuacion)<jupyter_output><empty_output><jupyter_text>Se pueden resolver incluso sistemas de ecuaciones:<jupyter_code>y = Symbol ( 'y' )
z = Symbol ( 'z' )
ecuacion1 = 2*x + 3*y - z + 6
ecuacion2 = 3*x + 2*y -12 
ecuacion3 = z - 5 
solve((ecuacion1 , ecuacion2, ecuacion3), dict = True ) 
import sympy as sy
from sympy import init_session
init_session(use_latex=True)
f = sy.sqrt( 1 - (abs(x) - 1) ** 2)
g = acos(1 - abs(x)) - pi 
p1 = plot(f ,(x , -2 , 2) , show = False , line_color= 'r')
p2 = plot(g ,(x , -2 , 2) , show = False , line_color= 'r')
p1.extend(p2)
p1.show()
from sympy.plotting import plot3d
a = x**3 + y**3
b = plot3d(a)
b.show()<jupyter_output><empty_output><jupyter_text>5. Ecuaciones diferenciales (Pyomo) ** Ecuación del calor ** $$ \pi ^2 \frac{\partial T}{\partial t} = \frac{\partial^2 T}{\partial t^2} $$Condiciones iniciales:$$ T(x,0) = sin(\pi x) $$Contorno:$$ T(0,t) = 0 $$$$ \pi e^{-t}+\frac{\partial T}{\partial t}(L,t) = 0 $$$$ 0 \leq t \leq 2 $$$$ 0 \leq l \leq L=1 $$<jupyter_code># Required imports
from pyomo.environ import *
from pyomo.dae import *
from pyomo.dae.plugins.finitedifference import Finite_Difference_Transformation

m = ConcreteModel()

m.pi = Param(initialize=3.1416)

#variables del problema
m.t = ContinuousSet(bounds=(0, 2))
m.x = ContinuousSet(bounds=(0, 1))
m.u = Var(m.x, m.t)

#variables derivadas respecto a las variables independientes x t
m.dudx  = DerivativeVar(m.u, wrt=m.x)
m.dudx2 = DerivativeVar(m.u, wrt=(m.x,m.x))
m.dudt  = DerivativeVar(m.u, wrt=m.t)<jupyter_output><empty_output><jupyter_text>$$ \pi ^2 \frac{\partial T}{\partial t} = \frac{\partial^2 T}{\partial t^2} $$<jupyter_code>def _pde(m, i, j):
    if i == 0 or i == 1 or j == 0:
        return Constraint.Skip
    return m.pi**2*m.dudt[i,j] == m.dudx2[i,j]
m.pde = Constraint(m.x, m.t, rule=_pde)<jupyter_output><empty_output><jupyter_text>$$ T(x,0) = sin(\pi x) $$<jupyter_code>def _initcon(m, i):
    if i == 0 or i == 1:
        return Constraint.Skip
    return m.u[i,0] == sin(m.pi*i)
m.initcon = Constraint(m.x,  rule=_initcon)<jupyter_output><empty_output><jupyter_text>$$ T(0,t) = 0 $$$$ \pi e^{-t}+\frac{\partial T}{\partial t}(L,t) = 0 $$<jupyter_code>def _lowerbound(m, j):
    return m.u[0,j] == 0
m.lowerbound = Constraint(m.t, rule=_lowerbound)

def _upperbound(m, j):
    return m.pi*exp(-j)+m.dudx[1,j] == 0
m.upperbound = Constraint(m.t, rule=_upperbound)<jupyter_output><empty_output><jupyter_text>**Función objetivo**<jupyter_code>m.obj = Objective(expr=1)<jupyter_output><empty_output><jupyter_text>$$ 0 \leq t \leq 2 $$$$ 0 \leq l \leq L=1 $$<jupyter_code>Nx = 60 #puntos de la coordenada espacial x
Nt = 25 #puntos de la coordenada temporal t

#discretizar puntos
discretize = Finite_Difference_Transformation()
disc = discretize.apply(m,nfe=Nx,wrt=m.x,scheme='BACKWARD')
disc = discretize.apply(m,nfe=Nt,wrt=m.t,scheme='BACKWARD')<jupyter_output><empty_output><jupyter_text>**Resolver**<jupyter_code>solver = SolverFactory('ipopt')
results = solver.solve(m, tee=True)<jupyter_output>******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
 Ipopt is released as open source code under the Eclipse Public License (EPL).
         For more information visit http://projects.coin-or.org/Ipopt
******************************************************************************

This is Ipopt version 3.12.4, running with linear solver mumps.
NOTE: Other linear solvers might be more efficient (see Ipopt documentation).

Number of nonzeros in equality constraint Jacobian...:    18452
Number of nonzeros in inequality constraint Jacobian.:        0
Number of nonzeros in Lagrangian Hessian.............:        0

Total number of variables............................:     6230
                     variables with only lower bounds:        0
                variables with lower and upper bounds:        0
                     variables with only upper bounds:        0
Total number of equ[...]<jupyter_text>**Representación de la solución**<jupyter_code>import numpy as np

x = sorted(disc.x)
t = sorted(disc.t)
u = sorted(disc.u)

U = np.zeros((Nx+1, Nt+1))
for i in range(Nx+1):
    ix = x[i]
    for j in range(Nt+1):
        it = t[j]
        U[i,j] = value(disc.u[ix,it])
#representacion grafica
from mpl_toolkits.mplot3d.axes3d import Axes3D
import matplotlib.pyplot as plt

from matplotlib import cm

fig = plt.figure(figsize=plt.figaspect(0.5))

ax = fig.add_subplot(1, 2, 1, projection='3d')
tt, xx = np.meshgrid(t, x)
#figura 1
surf = ax.plot_surface(xx, tt, U, rstride=1, cstride=1, cmap=cm.coolwarm,
                       linewidth=0, antialiased=False)
ax.set_zlim3d(0, 1.01)
fig.colorbar(surf, shrink=0.5, aspect=5)
ax.set_xlabel('x')
ax.set_ylabel('t')
ax.set_zlabel('u')
# figura 2
ax = fig.add_subplot(1, 2, 2, projection='3d')
ax.plot_wireframe(xx, tt, U)
ax.set_zlim3d(0, 1.01)
ax.set_xlabel('x')
ax.set_ylabel('t')
ax.set_zlabel('u')
plt.show()<jupyter_output><empty_output><jupyter_text>**Ecuaciones de ondas** $$ 2^2 \frac{\partial^2 u}{\partial x^2} = \frac{\partial^2 u}{\partial t^2} $$Condiciones iniciales:$$ u(x,0) = sin(\pi x)+sin(2\pi x) $$Contorno:$$ u(0,y) = 0 $$$$ u(L,t) = 0 $$$$ 0 \leq t \leq 5 $$$$ 0 \leq l \leq L = 1 $$<jupyter_code>! python ondas.py
! python ondasAnimate.py<jupyter_output>******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
 Ipopt is released as open source code under the Eclipse Public License (EPL).
         For more information visit http://projects.coin-or.org/Ipopt
******************************************************************************

This is Ipopt version 3.12.4, running with linear solver mumps.
NOTE: Other linear solvers might be more efficient (see Ipopt documentation).

Number of nonzeros in equality constraint Jacobian...:   108053
Number of nonzeros in inequality constraint Jacobian.:        0
Number of nonzeros in Lagrangian Hessian.............:        0

Total number of variables............................:    32878
                     variables with only lower bounds:        0
                variables with lower and upper bounds:        0
                     variables with only upper bounds:        0
Total number of equ[...]
{'fraction_non_alphanumeric': 0.14142819289878114, 'fraction_numerical': 0.030802861685214625, 'mean_word_length': 3.276770538243626, 'pattern_counts': {'":': 0, '<': 102, '<?xml version=': 0, '>': 111, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '28', 'dir': 'jupyter-structured-clean-dedup', 'id': '392563', 'license': 'MIT', 'n_tokens_mistral': 5919, 'n_tokens_neox': 5572, 'n_words': 1716, 'path': 'Mates_Python.ipynb', 'repo_name': 'Planelles20/PyConES-HaciendoMates', 'size': '15130'}