• Categories
• Matriz (Matemáticas)
• Estructura de datos de matriz
• Regresión lineal
• Linealidad
• Conceptos matemáticos

Citation preview

Programacion con Python Release 0.4

Oscar Vargas Torres

02 de September de 2013

Índice general

1 Programación básica 1.1 Introducción . . . . . . . . 1.2 Python básico . . . . . . . 1.3 Expresiones y operadores . 1.4 Procedimientos y funciones 1.5 Estructuras de control . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

3 3 8 9 10 11

2 Gráficas 2.1 Gráficas simples . . . . . . 2.2 Gráficas de contorno . . . . 2.3 Gráficas en 3D . . . . . . . 2.4 Ejercicios . . . . . . . . . . 2.5 Gráficas de algunos campos

. . . . . . . . . . . . . . . . . . . . simples

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

23 23 27 28 29 30

. . . . .

. . . . .

. . . . .

3 Álgebra lineal 33 3.1 Arreglos de NumPy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 3.2 Álgebra Lineal básica con NumPy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 3.3 Ejercicios . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 4 Polinomios e interpolación 39 4.1 Introducción a la Interpolación . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 4.2 Interpolación segmentaria (splines) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 5 Ajuste de funciones 43 5.1 Regresión por mínimos cuadrados . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 6 Raíces de ecuaciones no lineales 49 6.1 Introducción: Raíces de ecuaciones no lineales . . . . . . . . . . . . . . . . . . . . . . . . . . 49 6.2 Un vistazo a la librería GNU Scientific Library (GSL) . . . . . . . . . . . . . . . . . . . . . . 54 7 Diferenciación e integración numérica 59 7.1 Diferenciación numérica básica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 7.2 Integración numérica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60

i

8 Ecuaciones Diferenciales Ordinarias 65 8.1 Ecuaciones Diferenciales Ordinarias . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 8.2 PROYECTO FINAL: Rescatando a la tripulación del Apollo 13 . . . . . . . . . . . . . . . . 70 9 Índices y tablas

73

Índice

73

ii

Programacion con Python, Release 0.4

Contents:

Índice general

1

Programacion con Python, Release 0.4

2

Índice general

CAPÍTULO

1

Programación básica

1.1 Introducción 1.1.1 ¿Por qué aprender a programar? • Porque es una habilidad que aumenta tus posibilidades de adquirir un mejor empleo. • Para realizar cálculos intensos, es necesario utilizar una computadora, pues hacerlos a mano sería muy tardado o imposible. • Ejemplos de áreas en las que puedes aplicar la programación: – Métodos numéricos: Resuelva de manera numérica (utilizando el método de Euler) el sistema masa-resorte dado por la ecuación diferencial mu′′ + βu′ + ku = F (t),

u(0) = U0 , u′ (0) = 0

Sistema masa resorte que se describe mediante una ecuación diferencial – Graficación. A continuación se muestra un histograma obtenido a partir de una distribución de probabilidad normal.

3

Programacion con Python, Release 0.4

import matplotlib.pyplot as plt import numpy as np x = np.random.randn(1000) plt.hist(x, 20) plt.grid() plt.show()

160 140 120 100 80 60 40 20 04

3

2

1

0

1

2

3

1.1.2 Python como primer lenguaje de programación Se ha escogido Python como un lenguaje de programación porque combina capacidadades notables con una sintáxis limpia, simple y compacta. Enunciamos algunas características de Python: • Python es fácil de aprender y es muy apropiado para un curso introductorio de programación. • Es muy similar a Matlab y un buen lenguaje para hacer cómputo matemático. • Es fácil de combinar con lenguajes compilados como Fortran, C y C++, que son lenguajes ampliamente usados para cómputo científico.

1.1.3 Instalación de software utilizado en el curso Se requiere: • WinPython 3.3 → Puedes descargarlo e instalarlo de http://code.google.com/p/winpython/downloads/list

4

Capítulo 1. Programación básica

Programacion con Python, Release 0.4

1.1.4 Temario 1. Programación con software matemático. (a) Conceptos básicos de programación (con Python 3.x). Ejemplo: funciones. (b) Lectura y escritura de datos simple. 2. Gráficas con software matemático. (a) Matplotlib y Numpy. (b) Gráficas simples en 2D y en 3D. import matplotlib.pyplot as plt import numpy as np x = np.arange(-2, 2, 0.01) y = np.arange(-2, 2, 0.01) X, Y = np.meshgrid(x, y) ellipses = X*X/9 + Y*Y/4 - 1 cs = plt.contour(ellipses) plt.clabel(cs) plt.show()

0.000

0.2000.000 -0.200 350

0.200

-0.400

300

-0.800

250 200 150 100

-0.600

50 00

0.200

-0.400

-0.200 50

0.000 100 150

200

250

0.000 0.200 300 350

Gráfica de contornos. 3. Álgebra lineal. (a) Matrices y vectores. (b) Sistemas de ecuaciones lineales.

1.1. Introducción

5

Programacion con Python, Release 0.4 [ Para resolver el sistema de ecuaciones lineales Ax = b dado por A =

] [ ] 3 2 18 yb= : −1 2 2

import numpy as np A = np.array([[3.0, 2.0], [-1.0, 2.0]]) b = np.array([18.0, 2.0]) from numpy.linalg import solve x = solve(A, b) print(x) # Solución: x = [4.0, 3.0]

4. Polinomios e interpolación. import numpy as np from scipy.interpolate import interp1d import matplotlib.pyplot as plt x = np.linspace(0, 10, 10) y = np.cos(-x**2/8.0) f = interp1d(x, y) f2 = interp1d(x, y, kind='cubic') xnew = np.linspace(0, 10, 40) plt.plot(x, y, 'o', xnew, f(xnew), '-', xnew, f2(xnew), '--') plt.legend(['data', 'linear', 'cubic'], loc='best') plt.show()

5. Ajuste de funciones. (a) Regresión. (b) Interpolación. 6. Raíces de ecuaciones no lineales. Se tratará de hallar las raíces de ecuaciones. Es decir, el problema es: resolver f (x) = 0 para x. Ejemplo: determine la raíz más grande de h(x) = 2x3 − 11.7x2 + 17.7x − 5 en forma gráfica. import numpy as np import matplotlib.pyplot as plt x = np.linspace(0, 4) y = 2 * x**3 - 11.7 * x**2 + 17.7 * x - 5 plt.plot(x, y) plt.grid() plt.show()

6

Capítulo 1. Programación básica

Programacion con Python, Release 0.4

8 6 4 2 0 2 4 60.0

0.5

1.0

1.5

2.0

2.5

3.0

3.5

4.0

7. Funciones de splines e interpolación no lineal. 8. Diferenciación e integración numérica. 9. Problemas de valor inicial de ecuaciones diferenciales ordinarias. El movimiento vertical de un peso sujeto a un resorte que obedece la Ley de Hooke, y que tiene un estímulo senoidal externo se puede modelar con la ecuación diferencial: y¨ + 2cy˙ + ω02 y = F0 sin(ωt)

(1.1)

donde c, ω0 y ω son números reales positivos. Considéremos el caso en que c < ω0 , que corresponde a un sistema sub-amortiguado. Además suponga que ω < ω0 . Con estas suposiciones, la solución general de (1.1) se puede expresar como y = yu + yd , donde (√ ) (√ ) yu = C1 e−ct cos t ω02 − c2 + C2 e−ct sin t ω02 − c2 (1.2) yd = F0 M (ω) sin[ωt + ϕ(ω)]

(1.3)

C1 , C2 dependen de las condiciones iniciales y M (ω) = √

1

(ω02 − ω 2 )2 + 4c2 ω 2 ( ) −2cω ϕ(ω) = arctan − π < ϕ(ω) < 0 ω02 − ω 2

1.1. Introducción

(1.4) (1.5)

7

Programacion con Python, Release 0.4

Figura 1.1: Regla trapezoidal (Integración numérica)

1.2 Python básico Para iniciar con lo más básico, te recomiendo estudiar el Tutorial de Python (no es necesario que revises todo el tutorial, pero si es importante que empieces lo más pronto posible, y que vayas practicando cada tema que se vea en clase) en http://docs.python.org/3/tutorial/index.html (este tutorial está en inglés, pero es el tutorial adecuado para la versión de Python que estaremos usando en el curso). También puedes apoyarte en una versión del tutorial en español disponible en http://docs.python.org.ar/tutorial/3/index.html. Especialmente importantes son las secciones 1 a 5 del tutorial. Asegúrate de estudiar al menos estas secciones del tutorial.

1.2.1 Práctica 1 1. Descarga WinPython de http://code.google.com/p/winpython/downloads/list y ejecuta el instalador (Verifica primero si tienes una computadora con un sistema operativo de 64 bits o de 32 bits y descarga la versión más apropiada para ti. Si tienes dudas, pregunta a tu profesor cómo averiguar esto. Otra opción que tienes es descargar la versión de 32 bits, que puede ejecutarse en los dos tipos de computadoras). 2. Estudia con tu computadora el tutorial de Python, mencionado anteriormente. En otras palabras, ejecuta en tu computadora las instrucciones que se dan en esta sección a la vez que vas leyendo el tutorial. Anota las dudas que tengas para que puedas preguntar más adelante a tu profesor.

8

Capítulo 1. Programación básica

Programacion con Python, Release 0.4

1.2.2 Práctica 2 • Visita la página http://staff.not.iac.es/~rcardenes/hg/diveintopython3-es/native-datatypes.html y estudia los tipos de datos nativos que soporta Python 3 (puedes consultar la versión original en inglés en http://getpython3.com/diveintopython3/native-datatypes.html)

1.3 Expresiones y operadores 1.3.1 Práctica 3 1. Considere el movimiento vertical de una pelota que se lanza al aire. Dicho movimiento está dado por la fórmula 1 y(t) = v0 (t) − gt2 2

(1.6)

Con ayuda del siguiente bloque de código, verifica la posición de la pelota a los 0.6 segundos de ser lanzada, si la velocidad inicial con que fue lanzada es de 6 m/s: # Program for computing the height of a ball in vertical motion. v0 = 5 # initial velocity g = 9.81 # acceleration of gravity t = 0.6 #time y = v0*t - 0.5*g*t**2 # vertical position print(y)

¿Cuáles son los operadores que utilizas para la multiplicación y para la exponenciación? 2. Considere el movimiento vertical de la pelota dado por la ecuación (1.6). ¿Cuánto tiempo le toma a la pelota alcanzar la altura yc ? 1 2 gt − v0 t + yc = 0 2 • Usando la fórmula general para encontrar las raíces de una ecuación de segundo grado, comprueba que las soluciones son: • Comprueba que las soluciones están dadas por las raíces de la ecuación

( t1 =

v0 −

)

√ v02

− 2gyc /g,

( ) √ 2 t2 = v0 + v0 − 2gyc /g

(1.7)

• Utilizando el siguiente bloque de código, averigua cuáles son los dos instantes de tiempo en que una pelota está a una altura de 0.3 m, si la velocidad inicial de la pelota es de 6 m/s: v0 = 5 g = 9.81 yc = 0.2 import math t1 = (v0 - math.sqrt(v0**2 - 2*g*yc))/g t2 = (v0 + math.sqrt(v0**2 - 2*g*yc))/g print('En t={0:.4f} s y t={1:.4f} s, la altura es {2:.4f} m.'.format(t1, t2, yc))

3. Revisa los siguientes links • http://docs.python.org/3/library/numeric.html • http://docs.python.org/3/library/math.html

1.3. Expresiones y operadores

9

Programacion con Python, Release 0.4

• Utilizando el siguiente bloque de código, indica cuáles son los valores que toma la variable v: import math as m # m is now the name of the math module v = m.sin(m.pi) from math import log as ln v = ln(5) from math import sin as s, cos as c, log as ln v = s(x)*c(x) + ln(x)

4. Con ayuda del siguiente código, calcula las dos raíces complejas de la ecuación cuadrática f (x) = ax2 + bx + c cuando a = 1, b = 2 y c = 100: >>> a = 1; b = 2; c = 100 # polynomial coefficients >>> from numpy.lib.scimath import sqrt >>> r1 = (-b + sqrt(b**2 - 4*a*c))/(2*a) >>> r2 = (-b - sqrt(b**2 - 4*a*c))/(2*a)

1.4 Procedimientos y funciones 1.4.1 Práctica 4 1. Introduce la siguiente función en WinPython (o en IDLE, emacs, o el editor de tu preferencia): def mi_funcion(param1, param2): print param1 print param2

2. Llama la función anterior con los argumentos • 3 • “Hola” ¿Cuál es el tipo de estos argumentos? Utiliza la función type() 3. Convierte de metros a unidades de longitud inglesas. Construye un programa que reciba una longitud en metros y que entonces calcule y depliegue la correspondiente longitud en pulgadas, pies y millas. Utiliza las siguientes equivalencias: • 1 pulgada = 2.54 cm • 1 pie = 12 pulgadas • 1 yarda = 3 pies • 1 milla inglesa = 1760 yardas Como verificación puede usar los siguientes datos: • 640 m = 25196.85 pulgadas = 2099.74 pies = 699.91 yardas = 0.3977 millas inglesas 4. Modifica (en caso de ser necesario) el programa del ejercicio anterior de modo que utilices funciones para resolver el problema propuesto. 5. Basándote en el siguiente código, calcula el área de un paralelogramo, de un cuadrado y de un círculo, así como el volumen de un cono utilizando funciones:

10

Capítulo 1. Programación básica

Programacion con Python, Release 0.4

from math import pi h = 5.0 # altura b = 2.0 # base r = 1.5 # radio area_paralelogramo = h*b print("El area del paralelogramo es {0:.3f}".format(area_paralelogramo)) area_cuadrado = b**2 print("El area del cuadrado es {}".format(area_cuadrado)) area_circulo = pi*r**2 print("El area del circulo es {0:.3f}".format(area_circulo)) volumen_cono = 1.0/3*pi*r**2*h print("El volumen del cono es {0:.3f}".format(volumen_cono))

1.5 Estructuras de control 1.5.1 Estructuras selectivas Estudia las páginas 29 a 31 del libro “Python para todos”. Ejemplo 1 Construya un diagrama de flujo (vea la Figura Selección simple (if)) tal que dado como dato la calificación de un alumno en un examen, escriba “aprobado” en caso de que esa calificación sea mayor o igual a 70.

Figura 1.2: Selección simple (if)

1.5. Estructuras de control

11

Programacion con Python, Release 0.4

El siguiente código imprime “Aprobado”: cal = 80 if cal >= 70: print("Aprobado")

Ejemplo 2 Construya un diagrama de flujo (vea la Figura Selección doble (if else).) tal que dado como dato la calificación de un alumno en un examen, escriba “aprobado” si su calificación es mayor o igual a 70 y “reprobado” en caso contrario.

Figura 1.3: Selección doble (if

else).

El siguiente código imprime “Reprobado”: cal = 60 if cal >= 70: print("Aprobado") else: print("Reprobado")

Ejemplo 3 Construya un diagrama de flujo tal que dado como datos la categoría y el sueldo de un trabajador, calcule el aumento correspondiente teniendo en cuenta la siguiente tabla. Imprima la categoría del trabajador y su nuevo sueldo. Categoría 1 2 3 4

Aumento 15% 10% 8% 7%

El siguiente código resuelve el ejemplo:

12

Capítulo 1. Programación básica

Programacion con Python, Release 0.4

cat = 2 sue = 1000 if cat == 1: n_sue = sue * 1.15 elif cat == 2: n_sue = sue * 1.10 elif cat == 3: n_sue = sue * 1.08 else: n_sue = sue * 1.07 print("La categoría del trabajador es {}".format(cat)) print("El nuevo sueldo es {}".format(n_sue))

¿Cuál es el salario del trabajador (considerando que el salario actual es de $1000.00) para cada categoría? Ejercicios 1. Construya un diagrama de flujo tal que dado como dato un número entero, determine e imprima si el mismo es positivo, negativo o cero. Utilice la estructura de selección multiple. Después codifique su solución en Python. 2. En una tienda efectúan un descuento a los clientes dependiendo del monto de la compra. El descuento se efectúa con base en el siguiente criterio: • Si el monto es menor que $500.00, no hay descuento. • Si el monto está comprendido entre $500.00 y $1,000.00 (inclusive), 5% de descuento. • Si el monto está comprendido entre $1,000.00 y $7,000.00 (inclusive), 7% de descuento. • Si el monto está comprendido entre $7,000.00 y $15,000.00 (inclusive), 11% de descuento. • Si el monto es mayor a $15,000.00, 18% de descuento. Construya un diagrama de flujo que ilustre la solución del ejercicio. Después codifique su solución en Python. 3. En un estudio se ha hecho un estudio sobre los pacientes registrados durante los últimos 10 años, con el objeto de hacer una aproximación de los costos de internación por paciente. Se obtuvo un promedio diario según el tipo de enfermedad que aqueja al paciente. Además se pudo determinar que en promedio todos los pacientes con edad entre 14 y 22 años implican un costo adicional del 10%. La siguiente tabla expresa los costos diarios, según el tipo de enfermedad. Tipo de enfermedad 1 2 3 4

Costo/paciente/día 25 16 20 32

Construya un programa en Python que calcule e imprima el costo total que representa un paciente (pista: considere el número de días que el paciente permanece internado).

1.5.2 Estructuras cíclicas Estudia las páginas 32 a 35 del libro “Python para todos”

1.5. Estructuras de control

13

Programacion con Python, Release 0.4

Ejemplo 1 Imprima una tabla de conversión de grados Celsius a grados Fahrenheit. Dicha tabla debe verse como la siguiente: Celsius -20 -15 -10 -5 0 5 10 15 20 25 30 35 40

Fahrenheit -4.0 5.0 14.0 23.0 32.0 41.0 50.0 59.0 68.0 77.0 86.0 95.0 104.0

El siguiente código resuelve el ejemplo: print("------------------") # Encabezado de la tabla C = -20

# valor inicial de C

dC = 5

# incremento de C en el ciclo

while C 0 es el periodo de la función. f (t) puede escribirse como ( ) ( )] ∞ [ a0 ∑ 2nπt 2nπt f (t) = + an cos + bn sin 2 T T n=1

1.5. Estructuras de control

15

Programacion con Python, Release 0.4

Para encontrar los coeficientes a0 , an y bn , se utilizan las siguientes fórmulas: a0 = an =

2 T 2 T

2 bn = T

∫

a+T

f (t) dt ∫

a

(

a+T

f (t) cos ∫

a

(

a+T

f (t) sin a

2nπt T 2nπt T

) dt ) dt

El uso más frecuente de la serie de Fourier es para representar funciones periódicas. Lo que se hace es definir la función para un intervalo (a, a + T ) con una serie de Fourier. El resultado es una función definida en todo el eje real. a0 siempre da el valor promedio de la función. El valor que da la serie geométrica de Fourier para El valor 2 los valores en que la función presenta discontinuidades es el valor promedio de la función. Un desarrollo en serie de Fourier para f (t) definida de acuerdo con la figura import matplotlib.pyplot as plt import numpy as np x1 = np.array([0.0, 1.0]) y1 = np.array([2.0, 2.0]) x2 = np.array([1.0, 2.0]) y2 = np.array([-1.0, -1.0]) plt.plot(x1, y1, 'b', x2, y2, 'b') plt.axis([0.0, 2.0, -2.0, 3.0]) plt.grid() plt.show()

16

Capítulo 1. Programación básica

Programacion con Python, Release 0.4

3 2 1 0 1 20.0

0.5

1.0

1.5

2.0

está dado por las siguientes expresiones: a0 = 1 an = 0 3 n bn = [1 − (−1) ] nπ T =2 es decir, la solución es: ∞

f (t) =

1 ∑ 3 n + [1 − (−1) ] sin(nπt) 2 n=1 nπ ∞

6 1 ∑ sin [(2n − 1)πt] = + 2 n=1 (2n − 1)π =

∞ 1 6 ∑ sin [(2n − 1)πt] + 2 π n=1 2n − 1

La gráfica siguiente se obtiene utilizando los primeros 3 términos de la sumatoria de la última expresión para f (t): import matplotlib.pyplot as plt import numpy as np from numpy import pi

1.5. Estructuras de control

17

Programacion con Python, Release 0.4

num_puntos = 301 t = np.linspace(-2.0, 4.0, num_puntos) f = np.zeros((3, num_puntos)) for n in range(3): impar = 2*n + 1 theta = impar*pi*t f[n] = np.sin(theta)/impar plt.plot(t, f[n]) y = 0.5 + 6.0/pi*f.sum(0) plt.plot(t, y) plt.grid() plt.show()

2.5 2.0 1.5 1.0 0.5 0.0 0.5 1.0 1.5 2

1

0

1

2

3

4

Modifique el programa anterior para visualizar la señal periódica que se obtiene al sumar los primeros diez términos de la serie. Después visualize la señal correspondiente a 100 términos. Ejercicios 1. Escriba un programa que genere todos los números impares desde 1 hasta n. Fije un valor para n al inicio del programa y use un ciclo while para calcular los números. Asegúrese de que si n es un número par, el número impar más grande generado es n - 1. Nombre del programa: impar.py 2. Modifique el programa del ejercicio 1 y almacene los números impares generados en una lista. Empiece con una lista vacía y use un ciclo while donde en cada iteración del ciclo agregue un nuevo elemento a la 18

Capítulo 1. Programación básica

Programacion con Python, Release 0.4

lista. Finalmente, imprima la lista de elementos en la pantalla. Nombre del programa: lista_impares1.py 3. Resuelva el ejercicio 2 utilizando una list comprehension (con lista_impares2.py

for

y

range).

Nombre del programa:

4. (Calculo de una sumatoria) Se supone que el siguiente código calcule la suma s =

M ∑ 1 : k

k=1

s = 0 k = 1 M = 100 while k < M: s += 1/k print(s)

Este programa no funciona correctamente. ¿Cuáles son los tres errores? Si usted intenta ejecutar el programa, nada sucederá en la pantalla. Teclee Ctrl-C, es decir, presione (y deje presionada) la tecla Control (Ctrl) y entonces teclee la letra c, para detener el programa. Escriba un programa correcto. Nombre del programa: suma_while.py. Hay dos formas básicas de encontrar errores en el programa: (a) Lea el programa cuidadosamente y piense en las consecuencias de cada instrucción (b) Imprima resultados intermedios y compare con los cálculos manuales Primero, intente el método (a) y encuentre tantos errores como pueda. Entonces, intente el método (b) para M = 3 y compare la evolución de s con sus propios cálculos manuales. Nota: Una forma pythonica de resolver este problema es mediante: print(sum(1.0/k for k in xrange(1, M+1, 1)))

5. El siguiente programa grafica la función trigonométrica sin(x) import matplotlib.pyplot as plt import numpy as np from numpy import pi x = np.linspace(-2*pi, 2*pi, 101) y = np.sin(x) plt.plot(x, y) plt.grid() plt.show()

1.5. Estructuras de control

19

Programacion con Python, Release 0.4

1.0

0.5

0.0

0.5

1.0 8

4

6

2

0

2

4

6

8

Suponga que desea calcular la función seno inversa: arcsin(x). El módulo numpy tiene una función para calcular dicha función, ¿pero cuál es el nombre de esta función? Grafique esta función en el dominio adecuado. 6. Sea: q = [['a', 'b', 'c'], ['d', 'e', 'f'], ['g', 'h']]

Indexe esta lista para extraer (a) La letra

a

(b) La lista

['d', 'e', 'f']

(c) El último elemento (d) El elemento

h

d

(e) Explique por qué Nombre del programa:

q[-1][-2]

tiene el valor g.

index_nested_list.py

7. Construya un doble loop sobre una lista anidada. Considere la lista del ejercicio 6. Es posible visitar todos los elementos de q usando el siguiente for loop anidado: for i in q: for j in range(len(i)): print i[j]

¿Cuál es el tipo de las variables

20

i

y j?

Capítulo 1. Programación básica

Programacion con Python, Release 0.4

8. Reescriba la generación de la lista anidada q: q = [r**2 for r in [10**i for i in range(5)]]

utilizando ciclos estándar

for

en lugar de list comprehensions.

9. Grafique la siguiente función utilizando

matplotlib

 0    x N (x) =  2−x    0

si si si si

x>> a.dtype

Además del tipo de los elementos del arreglo, también es importante conocer su “forma”: >>> a >>> array([0, 1, 2, 3, 4]) >>> a.shape >>> (5,)

El atributo shape (forma) del arreglo es una tupla, en este caso, una tupla de 1 elemento, que contiene el tamaño del arreglo en cada dimensión. Creando un arreglo multidimensional Muestre la forma y el tipo de los elementos del arreglos siguientes: >>> m1 = array([arange(2), arange(2)]) >>> m2 = array([ [1, 2], [3, 4] ])

Ahora cree una matriz de tamaño 3 × 3 con los siguientes elementos   0 1 2 3 4 5 6 7 8 Compruebe que también puede crearlo de la siguiente manera: >>> m3 = arange(9).reshape(3, 3)

Los elementos de m3 pueden accederse mediante la siguiente sintáxis: >>> m3[0, 0]

¿Cuál es la salida de

m3[2, 1]?

3.2 Álgebra Lineal básica con NumPy El paquete numpy.linalg contiene funciones de álgebra lineal. Con este módulo, usted puede invertir matrices, calcular eigen–valores, resolver sistemas de ecuaciones lineales, calcular determinantes, entre otras cosas.

36

Capítulo 3. Álgebra lineal

Programacion con Python, Release 0.4

3.2.1 La inversa de una matriz La inversa de una matriz cuadrada A de n × n se denota por A−1 , que cuando se multiplica por la matriz original, es igual a la matriz identidad. Esto puede escribirse como se indica a continuación: AA−1 = In

A−1 A = In

La función inv del paquete numpy.linalg puede invertir una matriz por nosotros. Vamos a resolver un ejemplo: 1. Crea la matriz A con la función

mat:

>>> A = numpy.mat("0 1 2; 1 0 3; 4 -3 8")

2. Invierte la matriz A con la función

inv:

>>> A_inv = numpy.linalg.inv(A)

Si la matriz es singular, o no es cuadrada, resulta un error

LinAlgError.

3. Verifica utilizando multiplicación matricial: >>> print("Verifica\n {}".format(A * A_inv))

3.2.2 Resolviendo un sistema de ecuaciones lineales La función solve de numpy.linalg resuelve sistemas de ecuaciones de la forma Ax = b, donde A es una matriz, b puede ser un arreglo unidimensional o bidimensional, y x es la variable desconocida. 1. Crea la matriz A y el vector b: >>> A = numpy.mat("1 -2 1; 0 2 -8; -4 5 9") >>> print("A\n {}".format(A)) >>> b = numpy.array([0, -8, 9]) >>> print("b\n {}".format(b))

2. Invoque la función

solve

para resolver el sistema de ecuaciones lineales:

>>> x = numpy.linalg.solve(A, b) >>> print("Solucion: {}".format(x))

3. Verifique que la solución es correcta (dot realiza una multiplicación matricial): >>> print("Verifica\n {}".format( numpy.dot(A, x) ))

3.3 Ejercicios Instrucciones: Utilice Numpy para resolver los siguientes ejercicios.     3.2 1.5 1. Resuelva las operaciones indicadas, suponiendo que ⃗u = −0.6, ⃗v =  4.1  −1.4 −0.2 (a) ⃗u · ⃗v (b) ∥⃗u∥ (c) 3⃗v − ⃗u 3.3. Ejercicios

37

Programacion con Python, Release 0.4

(d) ⃗u × ⃗v 2. Encuentre la solución del sistema de ecuaciones lineales dado por las matrices aumentadas   0 2 3 8 1 5 (a) 2 3 1 −1 −2 −5   1 1 1 1 4 1 2 3 4 10  (b)  1 3 6 10 20 1 4 10 20 35 3. Utilice la función numpy.linalg.eig para encontrar los eigenvalores y eigenvectores de la matriz cua  1 0 2 drada A = 3 −1 3 2 0 1

38

Capítulo 3. Álgebra lineal

CAPÍTULO

4

Polinomios e interpolación

4.1 Introducción a la Interpolación La interpolación se usa para estimar valores intermedios entre datos precisos. El método más común que se usa para este propósito es la interpolación polinomial. Para n + 1 puntos, hay uno y sólo un polinomio de orden n que pasa a través de todos los puntos. Por ejemplo, hay sólo una línea recta (un polinomio de primer orden) que conecta dos puntos. De manera similar únicamente una parábola conecta un conjunto de tres puntos. La interpolación polinomial consiste en determinar el único polinomio de n–ésimo orden que ajuste n + 1 puntos. Este polinomio entonces proporciona una fórmula para calcular valores intermedios.

4.1.1 Interpolación de polinomios de Lagrange con SciPy 4.1.2 Ejemplo Estime log10 5 mediante interpolación lineal. • Interpole entre log 4 = 0.60206 y log 6 = 0.7781513: import numpy as np from scipy.interpolate import lagrange x = np.array([4.0, 6.0]) y = np.array([0.60206, 0.7781513]) p = lagrange(x, y) print(p) print(p(4.0)) print(p(6.0)) print(p(5.0))

• Interpole entre log 4.5 = 0.6532125 y log 5.5 = 0.7403627. • Ajuste con un polinomio de interpolación de Lagrange de tercer orden para estimar log 5 utilizando los datos anteriores.

39

Programacion con Python, Release 0.4

4.1.3 El fenómeno de Runge Ejemplo Consideremos el polinomio interpolador de Lagrange de la función f (x) = 1/(1 + 12x2 ) en el intervalo [−1, 1] para nodos equiespaciados, ¿tiende a cero el error EN (x) = F (x) − PN (x) cuando N crece? No: el máximo del término del error EN (x) crece cuando N → ∞. Esta falta de convergencia se conoce como fenómeno de Runge import numpy as np import matplotlib.pyplot as plt from scipy.interpolate import lagrange def f(x): return 1.0/(1.0 + 12.0 * x**2) x = np.linspace(-1.0, 1.0, 200) plt.plot(x, f(x), '--') x2 = np.linspace(-1.0, 1.0, 11) p = lagrange(x2, f(x2)) y2 = p(x2) plt.plot(x2, y2, 'o') plt.plot(x, p(x), linewidth=2.0) plt.grid() plt.show()

Una mejor solución utilizando polinomios de Chebyshev import numpy as np import matplotlib.pyplot as plt from numpy.polynomial import Chebyshev as T def f(x): return 1.0/(1.0 + 12.0 * x**2) N = 11 x = np.cos([(2*k + 1)*np.pi/2.0/N for k in range(N)]) y = f(x) p = T.fit(x, y, 10) xx, yy = p.linspace() plt.plot(xx, f(xx), '--') plt.plot(x, y, 'o') plt.plot(xx, yy, lw=2) plt.grid() plt.show()

40

Capítulo 4. Polinomios e interpolación

Programacion con Python, Release 0.4

1.0 0.8 0.6 0.4 0.2 0.01.0

0.5

0.0

0.5

1.0

4.2 Interpolación segmentaria (splines)

4.2. Interpolación segmentaria (splines)

41

Programacion con Python, Release 0.4

42

Capítulo 4. Polinomios e interpolación

CAPÍTULO

5

Ajuste de funciones

Veremos ajuste por mínimos cuadrados.

5.1 Regresión por mínimos cuadrados 5.1.1 Regresión lineal Es el ajuste de un conjunto de pares de observaciones: (x1 , y1 ), (x2 , y2 ), ldots, (xn , yn ) a una línea recta. La expresión matemática para esta última es: y = a0 + a1 x + e

donde a0 y a1 son coeficientes que representan el intercepto y la pendiente, respectivamente, y e es el error, o residuo, entre el modelo y las observaciones, el cual se puede representar al reordenar la expresión anterior como: e = y − a0 − a1 x

(5.1)

Así, el error o residuo es la discrepancia entre el valor real de y y el valor aproximado, a0 + a1 x, predicho por la ecuación lineal. En la regresión lineal por mínimos cuadrados, la estrategia es (como su nombre lo indica), minimizar la suma de los cuadrados de los residuos entre la y medida y la y calculada con el modelo lineal: Sr =

n ∑

e2i =

i=1

n ∑

(yi − a0 − a1 xi )

2

(5.2)

i=1

donde n es el número de pares de observaciones. Ecuaciones normales (regresión lineal) Para determinar los valores de a0 y a1 que minimizan la Ec. (5.2), esta última es diferenciada parcialmente con respecto a cada coeficiente e igualando cada expresión resultante a cero. Como resultado se tienen las

43

Programacion con Python, Release 0.4

ecuaciones normales:

(∑ ) ∑ na0 + x i a1 = yi (∑ ) (∑ ) ∑ xi a0 + x2i a1 = xi yi

(5.3) (5.4)

que pueden resolverse en forma simultánea.

5.1.2 Ejemplo 1: Lo básico Ajustemos una línea recta y = mx + b, a través de unos puntos con un poco de ruido: import numpy as np x = np.array([ 0,

1,

2,

3])

y = np.array([-1, 0.2, 0.9, 2.1]) # Es posible escribir la ecuación lineal como y = Ap, donde # A = [[x 1]] A = np.vstack([x, np.ones(len(x))]).T print(A) # y # p = [[m], [b]] # Use np.linalg.lstsq para calcular p: m, b = np.linalg.lstsq(A, y)[0] print(m, b) # Comprueba que # m = 1.0 y b = -0.95

Ahora vamos a graficar los datos junto con la línea ajustada: import matplotlib.pyplot as plt import numpy as np x = np.array([ 0,

1,

2,

3])

y = np.array([-1, 0.2, 0.9, 2.1]) m, b = 1.0, -0.95 plt.plot(x, y, 'o', label='Datos originales', markersize=10) plt.plot(x, m*x + b, 'r', label='Datos ajustados') plt.legend(loc='best') plt.show()

44

Capítulo 5. Ajuste de funciones

Programacion con Python, Release 0.4

2.5 2.0

Datos originales Datos ajustados

1.5 1.0 0.5 0.0 0.5 1.00.0

0.5

1.0

1.5

2.0

2.5

3.0

Resuelva el ejercicio manualmente y compruebe sus resultados.

5.1.3 Forma general lineal por mínimos cuadrados La regresión lineal, polinomial y lineal múltiple pertenecen al siguiente modelo general de mínimos cuadrados lineales: y = a0 z0 + a1 z1 + a2 z2 + · · · + am zm + e

(5.5)

donde z0 , z1 , . . ., zm son las m + 1 funciones diferentes. La regresión lineal simple y múltiple encajan dentro de este modelo: z0 = 1, z1 = x1 , . . ., zm = xm . Además la regresión de polinomios se incluye también si las z son monomios simples como en z0 = x0 = 1, z1 = x, . . ., zm = xm . Observe que la terminología lineal se refiere sólo a la dependencia del modelo sobre sus parámetros (es decir, las a). Como en el caso de regresión de polinomios, las mismas funciones pueden ser altamente no lineales. La ecuación (5.5) puede se puede expresar en notación matricial como: Y = ZA + E

(5.6)

donde Z es una matriz d los valores calculados de las funciones z en los valores medidos de las variables

5.1. Regresión por mínimos cuadrados

45

Programacion con Python, Release 0.4

independientes 

z01  z02  Z= .  ..

z11 z12 .. .

··· ··· .. .

 zm1 zm2   ..  . 

z0n

z1n

···

zmn

donde m es el número de variables del modelo y n es el número de datos. Como n ≥ m + 1, la mayoría de las veces Z no es una matriz cuadrada. El vector columna Y contiene los valores observados de la variable dependiente:   y1  y2    Y = .   ..  yn

El vector columna A contiene los coeficientes desconocidos   a0  a1    A= .   ..  am

y el vector columna E contiene los residuos 

 e1  e2    E= .   ..  an

La suma de los cuadrados de los residuos para este modelo se puede definir como

Sr =

n ∑ i=1

e2i =

n ∑

 yi −

i=1

m ∑

2 aj zji 

j=0

Esta cantidad se puede minimizar al tomar su derivada parcial con respecto a cada uno de los coeficientes y fijar los resultados de la ecuación igual a cero. La salida de este proceso son las ecuaciones normales que se pueden expresar brevemente en forma de matriz como Z T ZA = Z T Y

(5.7)

Ejemplo 2: Regresión polinomial (de segundo orden) Ajustar a una parábola los datos siguientes:

46

Capítulo 5. Ajuste de funciones

Programacion con Python, Release 0.4

xi 0 1 2 3 4 5

yi 2.1 7.7 13.6 27.2 40.9 61.1

La matriz Z debe tomar la forma  1 1  Z = .  ..

x1 x2 .. .

 x21 x22   ..  . 

(5.8)

1 xn

x2n

que es un caso particular de la matriz de Vandermonde  1 x1 x21 1 x2 x22  Z = . .. ..  .. . .

··· ···

 xm 1  xm 2  ..  . 

x2n

···

xm n

1

xn

Podemos definir una función en python que calcule la matriz de Vandermonde de la siguiente manera: import numpy as np def vandermonde(x, N): return np.column_stack([x**(i) for i in range(N + 1)])

Compruebe que llamar esta función con [[ 1

0

0]

[ 1

1

1]

[ 1

2

4]

[ 1

3

9]

[ 1

4 16]

[ 1

5 25]]

x = np.arange(6)

y

N = 2

genera el siguiente arreglo de numpy:

Por lo tanto, el problema planteado puede resolverse planteando las ecuaciones normales (5.7): import numpy as np def vandermonde(x, N): return np.column_stack([x**(i) for i in range(N + 1)]) Z = vandermonde(np.arange(6), 2) Y = np.array([2.1, 7.7, 13.6, 27.2, 40.9, 61.1]) A = np.linalg.solve( np.dot(Z.T, Z), np.dot(Z.T,Y) ) print(A)

Utilizando

numpy.vander

y

numpy.linalg.lstsq,

el problema puede resolverse de la siguiente manera:

import numpy as np Z = np.vander(np.arange(6), 3) Y = np.array([2.1, 7.7, 13.6, 27.2, 40.9, 61.1]) A = np.linalg.lstsq(Z, Y)[0] print(A)

5.1. Regresión por mínimos cuadrados

47

Programacion con Python, Release 0.4

Interprete los resultados anteriores y grafique los datos junto con el resultado obtenido de la regresión polinomial

48

Capítulo 5. Ajuste de funciones

CAPÍTULO

6

Raíces de ecuaciones no lineales

6.1 Introducción: Raíces de ecuaciones no lineales 6.1.1 Método gráfico Aunque el método gráfico no es el más exacto para hallar las raíces de las ecuaciones no lineales, resulta de ayuda. Ejemplo 1 Encuentre la raíz de [

] 1 − exp(−0.146843c) f (c) = 667.38 − 40 c

en el intervalo de [0.5, 18.0] import numpy as np import matplotlib.pyplot as plt c = np.linspace(0.5, 18.0) f = 667.38*(1 - np.exp(-0.146843*c))/c - 40 plt.plot(c, f) plt.grid() plt.show()

49

Programacion con Python, Release 0.4

60 50 40 30 20 10 0 100

2

4

6

8

10

12

14

16

18

Puede notarse que la solución está entre 14 y 16. Si grafica en ese intervalo, puede notarse un poco mejor la raíz. import numpy as np import matplotlib.pyplot as plt c = np.linspace(14.0, 16.0) f = 667.38*(1 - np.exp(-0.146843*c))/c - 40 plt.plot(c, f) plt.grid() plt.show()

50

Capítulo 6. Raíces de ecuaciones no lineales

Programacion con Python, Release 0.4

2.0 1.5 1.0 0.5 0.0 0.5 1.0 1.5 2.0 2.5 14.0

14.5

15.0

15.5

16.0

Se puede utilizar Scipy para encontrar las raíces de manera numérica, con mucha simplicidad. A continuación se muestra cómo se resuelve el problema anterior con el método de bisección, uno de los más simples para encontrar raíces de ecuaciones no lineales. import numpy as np from scipy.optimize import bisect def f(x): return 667.38*(1 - np.exp(-0.146843*x))/x - 40 x0 = bisect(f, 0.5, 18.0) print(x0)

6.1.2 Ejercicios 1. Determine la raíz real más grande de h(x) = 2x3 − 11.7x2 + 17.7x − 5 (a) En forma gráfica. import numpy as np import matplotlib.pyplot as plt x = np.linspace(0.0, 4.0) h = ((2*x - 11.7)*x

+ 17.7)*x - 5.0

plt.plot(x, h)

6.1. Introducción: Raíces de ecuaciones no lineales

51

Programacion con Python, Release 0.4

plt.grid() plt.show()

8 6 4 2 0 2 4 60.0

0.5

1.0

1.5

2.0

2.5

3.0

3.5

4.0

import numpy as np import matplotlib.pyplot as plt x = np.linspace(3.5, 3.6) h = ((2*x - 11.7)*x

+ 17.7)*x - 5.0

plt.plot(x, h) plt.grid() plt.show()

52

Capítulo 6. Raíces de ecuaciones no lineales

Programacion con Python, Release 0.4

0.6 0.4 0.2 0.0 0.2 0.4 0.6 0.8 3.50

3.52

3.54

3.56

3.58

3.60

Puede notarse que la raíz más grande es aproximadamente 3.56. (b) Con el método de Newton-Raphson (tres iteraciones, x0 = 3, δ = 0.001). Este método es una aproximación a una raíz de f (x) = 0 a partir del valor inicial p0 mediante la iteración: pk = pk−1 −

f (pk−1 ) f ′ (pk−1 )

Un algoritmo para el método de Newton-Raphson se obtiene fácilmente al sustituir la ecuación anterior por la fórmula predictiva x = g(x) en el algoritmo de punto fijo. Observe, sin embargo, que el programa también debe modificarse para calcular la primera derivada. Considerando varios problemas potenciales con el método de Newton Raphson, el programa se podría mejorar incorporando algunas consideraciones adicionales (ver pág. 153, 154 del libro Métodos numéricos para ingenieros, de Chapra y Canale). def newton_raphson(g, dg, p0, delta, epsilon, num_max_iteraciones): # DATOS DE ENTRADA # g es la funcion de iteración. # dg es la derivada de g # p0 es la aproximacion inicial a un cero de f # delta es la tolerancia para la diferencia entre un valor de la #

iteracion y el valor de la iteracion anterior

# epsilon es la tolerancia para los valores de la funcion # numMaxIteraciones es el numero maximo de iteraciones # RESULTADOS: # err es la diferencia entre dos términos consecutivos

6.1. Introducción: Raíces de ecuaciones no lineales

53

Programacion con Python, Release 0.4

# P es la sucesión {pn} completa # k es el número de iteraciones realizadas # P[k] es la aproximación al punto fijo P = zeros(numMaxIteraciones + 1) P[0] = p0 for k in range(1, numMaxIteraciones + 1): gk = P[k - 1] - g(P[k - 1])/dg(P[k - 1]) P[k] = gk err = abs(P[k] - P[k - 1]) errorRelativo = 2*err/(abs(P[k]) + delta) y = g(P[k]) if (err < delta or errorRelativo < delta or abs(y) < epsilon): break if (k == numMaxIteraciones): print('Se ha excedido el numero de iteraciones') return P, k, errorRelativo def h(x): return ((2*x - 11.7)*x

+ 17.7)*x - 5.0

def dh(x): return (6.0*x - 23.4)*x + 17.7 P, k, error_relativo = newton_raphson(h, dh, 3.0, 0.001, 0.0001, 10)

(c) Utilizando el método

newton

del módulo

scipy.optimize.

def h(x): return ((2*x - 11.7)*x

+ 17.7)*x - 5.0

def dh(x): return (6.0*x - 23.4)*x + 17.7 from scipy.optimize import newton zero = newton(h, 3.0, dh) print(zero)

6.2 Un vistazo a la librería GNU Scientific Library (GSL) Aprovecharemos que la librería GSL tiene rutinas para encontrar las raíces de funciones unidimensionales arbitrarias que proporcionan acceso a los pasos intermedios de la solución, para estudiar cómo funcionan los métodos de intervalos. La librería GSL está escrita en C. No te preocupes si no comprendes todo el código del ejemplo siguiente (basado en el manual de referencia de la GNU GSL). Esta parte del curso es más ilustrativa y de ninguna manera exhaustiva. Lee con atención los comentarios que se dan en el código, pues te servirán de mucho para comprender la intención del mismo. /* * ===================================================================================== * * *

Nombre del archivo:

demo_fn.h

* *

54

Descripción:

Definición de

Capítulo 6. Raíces de ecuaciones no lineales

Programacion con Python, Release 0.4

*

+ Estructura parametros_cuadratica. Esta estructura almacena

*

los coeficientes {a, b, c} de un polinomio cuadrático

*

ax^2 + bx + c como números de tipo double.

*

Declaración de prototipos de funciones para

*

+ Evaluar un polinomio cuadrático, dado por los parámetros

*

apuntados por params, en un valor x (cuadratica)

*

+ Evaluación de la derivada del polinomio cuadrático

*

(deriv_cuadratica)

*

+ Evaluación simultánea del polinomio y de su derivada

*

(fdf_cuadratica)

* * ===================================================================================== */ struct parametros_cuadratica { double a, b, c; }; double cuadratica(double x, void *params); double deriv_cuadratica(double x, void *params); void fdf_cuadratica(double x, void *params, double *y, double *dy);

A continuación se muestra la implementación de los prototipos declarados en demo_fn.h /* * ===================================================================================== * *

Nombre del

*

archivo:

demo_fn.c

* *

Descripción:

Definición (implementación) de los prototipos en demo_fn.h

* * ===================================================================================== */ #include "demo_fn.h" double cuadratica(double x,

void *params)

{ struct parametros_cuadratica *p = (struct parametros_cuadratica *) params; double a = p->a; double b = p->b; double c = p->c; return (a * x + b) * x + c; } double deriv_cuadratica(double x, void *params) { struct parametros_cuadratica *p = (struct parametros_cuadratica *) params; double a = p->a; double b = p->b; double c = p->c; return 2.0 * a * x + b; }

6.2. Un vistazo a la librería GNU Scientific Library (GSL)

55

Programacion con Python, Release 0.4

void fdf_cuadratica(double x, void *params, double *y, double *dy) { struct parametros_cuadratica *p = (struct parametros_cuadratica *) params; double a = p->a; double b = p->b; double c = p->c; *y = (a * x + b) * x + c; *dy = 2.0 * a * x + b; }

Finalmente, el siguiente programa utiliza los dos archivos de código fuente anteriores para encontrar la raíz del polinomio cuadrático x2 − 5 = 0: /* * ===================================================================================== * *

Nombre del

*

archivo:

demo_root.c

* *

Descripción:

*

Ejemplo del uso de rutinas de GSL para encontrar las raíces de funciones unidimensionales arbitrarias.

*

Compilar con:

icc demo_root.c -lgsl -o demo_root

* ===================================================================================== */ #include #include #include #include #include "demo_fn.c" int main () { int estado; /* Indica el estado de una iteración del método */ int iter = 0; /* Contador de iteraciones */ const int MAX_ITER = 100; /* Número máximo de iteraciones */ const gsl_root_fsolver_type *T; /* Apuntador al _tipo_ de solucionador */ gsl_root_fsolver *s; /* Apuntador a la instancia del solucionador */ double r = 0.0; /* Valor inicial de la raíz */ double raiz_esperada = sqrt (5.0); /* Valores iniciales de los limites del intervalo de búsqueda */ double x_inferior = 0.0,

x_superior = 5.0;

/* * Estructura capaz de apuntar a *

+ la función cuya raíz desea encontrarse

*

+ los parámetros que necesita la función

*/ gsl_function F; struct parametros_cuadratica params = {1.0,

0.0,

-5.0};

F.function = &cuadratica; F.params = ¶ms; /* Se utiliza el método de bisección.

56

Capítulo 6. Raíces de ecuaciones no lineales

Programacion con Python, Release 0.4

* NOTA: Se pueden seleccionar entre diferentes tipos de solucionadores. *

Verifique el manual de la GSL

*/ T = gsl_root_fsolver_bisection; s = gsl_root_fsolver_alloc (T); /* Se crea la nueva instancia del solucionador */ /* * Se inicializa el solucionador s para que use la función F y * el intervalo inicial de búsqueda [x_inferior,

x_superior]

*/ gsl_root_fsolver_set (s,

&F,

x_inferior,

x_superior);

/* Imprimir el nombre del método usado */ printf ("Usando el método %s\n",

gsl_root_fsolver_name (s));

/* Imprimir el encabezado de resultados intermedios */ printf ("%5s [%9s, "iter",

%9s] %9s %10s %9s\n",

"inferior",

"superior",

"raiz",

"err",

"err(est)");

/* * Iterar mientras que * el estado sea GSL_CONTINUE y * no se haya alcanzado el número máximo de iteraciones */ do { iter++; /* Incrementar en uno el contador de iteraciones */ /* * Lleva a cabo una iteración del solucionador s. * Si la iteración encuentra un error inesperado, la variable estado * contiene un código de error */ estado = gsl_root_fsolver_iterate (s); /* Actualiza la mejor estimación de la raíz que se tiene al momento */ r = gsl_root_fsolver_root (s); /* Actualiza los límites del mejor intervalo que encierra una posible raíz */ x_inferior = gsl_root_fsolver_x_lower (s); x_superior = gsl_root_fsolver_x_upper (s); /* * Verificar la convergencia del método en el intervalo con * error absoluto = 0.0 * error relativo = 0.001 * El test devuelve GSL_SUCCESS si el método ha convergido */ estado = gsl_root_test_interval (x_inferior,

x_superior, 0, 0.001);

if (estado == GSL_SUCCESS) printf ("El método converge:\n"); /* Mostrar información sobre la iteración actual */ printf ("%5d [%.7f,

%.7f] %.7f %+.7f %.7f\n",

iter, x_inferior, x_superior, r, r - raiz_esperada, x_superior - x_inferior); } while (estado == GSL_CONTINUE && iter < MAX_ITER); /* Liberar la memoria asociada al solucionador s */ gsl_root_fsolver_free (s);

6.2. Un vistazo a la librería GNU Scientific Library (GSL)

57

Programacion con Python, Release 0.4

/* Terminar el programa devolviendo el valor de la variable estado */ return estado; } /* ----------

Fin de la función main

---------- */

El resultado de ejecutar el programa anterior es: Usando el método bisection iter [ inferior,

superior]

raiz

err

err(est)

1 [0.0000000,

2.5000000] 1.2500000 -0.9860680 2.5000000

2 [1.2500000,

2.5000000] 1.8750000 -0.3610680 1.2500000

3 [1.8750000,

2.5000000] 2.1875000 -0.0485680 0.6250000

4 [2.1875000,

2.5000000] 2.3437500 +0.1076820 0.3125000

5 [2.1875000,

2.3437500] 2.2656250 +0.0295570 0.1562500

6 [2.1875000,

2.2656250] 2.2265625 -0.0095055 0.0781250

7 [2.2265625,

2.2656250] 2.2460938 +0.0100258 0.0390625

8 [2.2265625,

2.2460938] 2.2363281 +0.0002601 0.0195312

9 [2.2265625,

2.2363281] 2.2314453 -0.0046227 0.0097656

10 [2.2314453,

2.2363281] 2.2338867 -0.0021813 0.0048828

11 [2.2338867,

2.2363281] 2.2351074 -0.0009606 0.0024414

El método converge: 12 [2.2351074,

2.2363281] 2.2357178 -0.0003502 0.0012207

Verifica manualmente los resultados de la tabla anterior.

58

Capítulo 6. Raíces de ecuaciones no lineales

CAPÍTULO

7

Diferenciación e integración numérica

7.1 Diferenciación numérica básica • Aproximación a la primera derivada con diferencia hacia adelante f ′ (xi ) =

f (xi+1 ) − f (xi ) + O(xi+1 − xi ) xi+1 − xi

• Aproximación a la primera derivada con diferencia hacia atrás f ′ (xi ) =

f (xi ) − f (xi−1 ) + O(xi − xi−1 ) xi − xi−1

• Aproximación a la primera derivada con diferencias centradas f ′ (xi ) =

( ) f (xi+1 ) − f (xi−1 ) + O h2 2(xi+1 − xi−1 )

El siguiente bloque de código resuelve el problema: def derivadasConDifDivididas(f, xi, h): difHaciaAdelante = (f(xi + h) - f(xi))/h print("Con diferencia hacia adelante: {}".format(difHaciaAdelante)) difHaciaAtras = (f(xi) - f(xi - h))/h print("Con diferencia hacia atras: {}".format(difHaciaAtras)) difCentradas =

(f(xi + h) - f(xi - h))/(2*h)

print("Con diferencia centrada: {}".format(difCentradas)) def f(x): return (((-0.1*x - 0.15)*x - 0.5)*x - 0.25)*x + 1.2 derivadasConDifDivididas(f, 0.5, 0.5) derivadasConDifDivididas(f, 0.5, 0.25)

59

Programacion con Python, Release 0.4

Utilizando Scipy: from scipy.misc import derivative def f(x): return (((-0.1*x - 0.15)*x - 0.5)*x - 0.25)*x + 1.2 df = derivative(f, 0.5, 0.25)

7.2 Integración numérica Algunas funciones pueden ser integradas de manera analítica, mientras que muchas otras no, y entonces se debe recurrir a una aproximación numérica a la integral. Por ejemplo, en el campo de la termodinámica estadística, el modelo de Debye para calcular la capacidad calórica de un sólido considera la siguiente función ∫ x t3 Φ(x) = dt t 0 e −1 Puesto que ho hay una expresión analítica para Φ(x), debemos usar algún método de integración numérica para calcular sus valores.

7.2.1 Ejemplo 1 ∫

4

x2 dx utilizando

Calcule

scipy.integrate

y compare con el resultado analítico:

0

from scipy import integrate x2 = lambda x: x**2 integ = integrate.quad(x2, 0.0, 4.0) print integ

Ahora calcule de manera numérica: 1. Φ(5) ∫ ∞ 2. e−x dx. Pista: utilice

inf

para denotar +∞.

0

La regla trapezoidal La regla trapezoidal es la más sencilla de todas las fórmulas de integración cerrada de Newton-Cotes. Geométricamente, la regla trapezoidal es equivalente a aproximar el área bajo la curva con la del trapecio bajo la línea recta que conecta f (a) y f (b) en la Figura Regla trapezoidal. Recuerde que la fórmula para calcular el área de un trapecio es la altura por el promedio de las bases. En nuestro caso, el concepto es el mismo, pero el trapecio está sobre su lado. La regla compuesta del trapecio Un método intuitivo para hallar el área limitada por una curva y = f (x) en un intervalo [a, b] es dar una aproximación a dicha área sumando las áreas de una serie de trapecios construidos sobre los intervalos {[xk , xk+1 ]} de una partición de [a, b]. 60

Capítulo 7. Diferenciación e integración numérica

Programacion con Python, Release 0.4

Figura 7.1: Regla trapezoidal

from pylab import * from matplotlib.patches import Polygon def func(x): return (x-3) * (x-5) * (x-7) + 85 ax = subplot(111) a, b = 2, 9 # Area de la integral x = arange(0, 10, 0.01) y = func(x) plot(x, y, linewidth=1) # make the shaded region ix = linspace(a, b, 4) iy = func(ix) verts = [(a,0)] + zip(ix,iy) + [(b,0)] poly = Polygon(verts, facecolor='0.8', edgecolor='k') ax.add_patch(poly) markerline, stemlines, baseline = stem(ix, iy, 'k-.') setp(markerline, 'markerfacecolor', 'b') axis([0,10, 0, 180]) figtext(0.9, 0.05, 'x') figtext(0.1, 0.9, 'y') ax.set_xticks((a,b)) ax.set_xticklabels(('a','b'))

7.2. Integración numérica

61

Programacion con Python, Release 0.4

ax.set_yticks([]) show()

Supongamos que se divide el intervalo [a, b] en M subintervalos [xk , xk+1 ] de anchura común h = (b − a)/M mediante una partición cuyos nodos xk = a + kh, para k = 0, 1, . . . , M , están equiespaciados. La regla compuesta del trapecio con M subintervalos se puede expresar de las siguientes formas equivalentes: h∑ (f (xk−1 ) + f (xk )) 2

(7.1)

h (f0 + 2f1 + 2f2 + 2f3 + · · · + 2fM −2 + 2fM −1 + fM ) 2

(7.2)

M −1 ∑ h (f (a) − f (b)) + h f (xk ) 2

(7.3)

M

T (f, h) =

k=1

o bien T (f, h) =

o bien T (f, h) =

k=1

Este valor es una aproximación a la integral de f (x) en [a, b], lo que se escribe como ∫

b

f (x) dx ≈ T (f, h) a

7.2.2 Ejemplo 2 √ Consideremos f (x) = 2 + sin(2 x). Vamos a usar la regla compuesta del trapecio con 11 nodos para calcular una aproximación a la integral de f (x) en el intervalo [1, 6]. Para los once nodos, tomamos M = 10, con lo que h = (6−1)/10 = 1/2. Usando la fórmula (7.3), los cálculos son: ( ) 1 1/2 T f, = (f (1) + f (6)) + 2 2 ( ( ) ( ) ( ) ( ) ( )) 1 3 5 7 9 11 f + f (2) + f + f (3) + f + f (4) + f + f (5) + f 2 2 2 2 2 2 1 = (2.9093 + 1.01735) + 4 1 (2.63815 + 2.30807 + 1.97931 + 1.68305 + 1.43530 + 1.24319 + 1.10831 + 1.02872 + 1.00024) 2 1 1 = (3.92665) + (14.42438) 4 2 = 0.98166 + 7.21219 = 8.19385 Compruebe el resultado anterior manualmente y después utilizando el siguiente código:

62

Capítulo 7. Diferenciación e integración numérica

Programacion con Python, Release 0.4

from numpy import linspace, sin, sqrt from scipy.integrate import trapz x = linspace(1, 6, 11) y = 2 + sin( 2*sqrt(x) ) T = trapz(y, x) print(T)

7.2. Integración numérica

63

Programacion con Python, Release 0.4

64

Capítulo 7. Diferenciación e integración numérica

CAPÍTULO

8

Ecuaciones Diferenciales Ordinarias

8.1 Ecuaciones Diferenciales Ordinarias Considere la ecuación: dy = 1 − e−t dt Esta ecuación es separable, por lo que puede resolverse directamente mediante integración: y(t) = t + e−t + C

donde C es la constante de integración. En esta solución existe un grado de libertad, dado por la elección de la constante de integración C. import matplotlib.pyplot as plt import numpy as np num_puntos = 81 t = np.linspace(-2.0, 5.0, num_puntos) num_curvas = 5 y = np.zeros((num_curvas, num_puntos)) for i, C in enumerate(range(-2, num_curvas - 2)): y[i] =

t + np.exp(-t) + C

plt.plot(t, y[i]) plt.grid() plt.show()

65

Programacion con Python, Release 0.4

8 7 6 5 4 3 2 1 0 12

1

0

1

2

3

4

5

A diferencia de la ecuación diferencial dada, no todas pueden resolverse de manera exacta (ya sea mediante integración directa o por alguno de los métodos “tradicionales”). A menudo se da el caso de que no hay una soluciñon analítica conocida, por lo que necesitamos aproximaciones numéricas. A modo de ejemplo, consideremos, en el contexto de la dinámica de poblaciones, un sistema no lineal que es una modificación de las ecuaciones de Lotka-Volterra 1 2 x 10 1 y˙ = g(t, x, y) = xy − y − y 2 20

x˙ = f (t, x, y) = x − xy −

para 0 ≤ t ≤ 30 con la condición inicial x(0) = 2 e y(0) = 1. La solución numérica se consigue con el siguiente código: import matplotlib.pyplot as plt import numpy as np from scipy.integrate import odeint x0 = [2, 1] def f(xVec, t): x, y = xVec f0 = x - x*y - x**2 / 10.0 f1 = x*y - y - y**2 / 20.0 return [f0, f1]

66

Capítulo 8. Ecuaciones Diferenciales Ordinarias

Programacion con Python, Release 0.4

t = np.linspace(0.0, 30.0, 3001) xVec = odeint(f, x0, t) x = xVec[:, 0] y = xVec[:, 1] plt.plot(x, y) plt.show()

8.1.1 Mezclas de soluciones salinas La mezcla de dos soluciones salinas de distintas concentraciones da lugar a una ecuación diferencial de primer orden, que define la cantidad de sal que contiene la mezcla. Supongamos que un tanque está parcialmente lleno con 100 galones de salmuera (sal disuelta en agua), con 1 lb de sal por galón a razón de 6 gal/min. El contenido del 10 lb de sal disuelta. Le entra salmuera con 2 tanque está bien mezclado y de él sale a razón de 6 gal/min. El contenido del tanque está bien mezclado y de él sale a razón de 4 gal/min de solución. Calcule la cantidad de libras de sal que hay en el tanque a los 30 minutos. ( )( ) 1 lb gal lb Ri = 6 =3 2 gal min min )( ) ( lb gal 2A lb A 4 = Ro = 100 + 2t gal min 50 + t min dA 2A = Ri − Ro = 3 − dt 50 + t La ecuación diferencial anterior se puede expresar en la forma estándar para las ecuaciones lineales: dA 2 + A=3 dt 50 + t

(8.1)

El factor integrante es: ∫

e

2 50+t

dt

= e2 ln(50+t) = (50 + t)2

(8.2)

Multiplicando la Ec. (8.1) por el factor integrante (8.2) ( ) dA 2 2 (50 + t) + A = 3(50 + t)2 dt 50 + t

identificamos el lado izquierdo como

] d [ (50 + t)2 A , por lo que tenemos la expresión más simple: dt ] d [ (50 + t)2 A = 3(50 + t)2 dt

que hay que resolver por integración de ambos miembros: ∫ (50 + t)2 A = 3 (50 + t)2 dt = (50 + t)3 + C

Sustituyendo la condición inicial A(0) = 10, en la expresión anterior, obtenemos C = −100000. La solución de la Ec. (8.1) es: A(t) = 50 + t −

8.1. Ecuaciones Diferenciales Ordinarias

100000 (50 + t)2

67

Programacion con Python, Release 0.4

Por lo tanto, a los 30 min., la cantidad de sal disuelta en el tanque es A(30) = 64.375 lb. Resuelva de manera numérica el problema. Grafique la solución exacta junto con la solución numérica, hasta los 35 min.

8.1.2 Un resorte de Hooke con estímulo senoidal El movimiento vertical de un peso sujeto a un resorte que obedece la Ley de Hooke, y que tiene un estímulo senoidal externo se puede modelar con la ecuación diferencial: y¨ + 2cy˙ + ω02 y = F0 sin(ωt)

(8.3)

donde c, ω0 y ω son números reales positivos. Considéremos el caso en que c < ω0 , que corresponde a un sistema sub-amortiguado. Además suponga que ω < ω0 . Con estas suposiciones, la solución general de (8.3) se puede expresar como y = yu + yd , donde (√ ) (√ ) (8.4) yu = C1 e−ct cos t ω02 − c2 + C2 e−ct sin t ω02 − c2 yd = F0 M (ω) sin[ωt + ϕ(ω)]

(8.5)

C1 , C2 dependen de las condiciones iniciales y 1 M (ω) = √ 2 (ω0 − ω 2 )2 + 4c2 ω 2 ( ) −2cω ϕ(ω) = arctan − π < ϕ(ω) < 0 ω02 − ω 2

(8.6) (8.7)

Ejemplo El sistema y¨ + y˙ + 4y = 12 sin t

con c = 1/2, ω0 = 2, ω = 1, sujeto a las condiciones iniciales y(0) = 0, y(0) ˙ = 0 tiene la solución [ (√ ) ( √ )] 15 2√ 15 6√ −t/2 6 y(t) = e cos t − 15 sin t + 10 sin(t − 0.32175) 5 2 5 2 5

import matplotlib.pyplot as plt import numpy as np tfinal = 20.0 num_puntos = tfinal * 100 + 1 t = np.linspace(0.0, tfinal, num_puntos) theta = np.sqrt(15)*t/2.0 yu = np.exp(-0.5*t)*(1.2*np.cos(theta) - 1.54919*np.sin(theta)) yd = 1.2*np.sqrt(10)*np.sin(t - 0.32175) y = yu + yd plt.plot(t, y) plt.grid() plt.show()

68

Capítulo 8. Ecuaciones Diferenciales Ordinarias

Programacion con Python, Release 0.4

4 3 2 1 0 1 2 3 40

5

10

15

20

Una solución numérica de este problema, utilizando un método de Runge-Kutta de orden (4)5 debido a Dormand y Prince (con control del tamaño de paso y salida densa) está dada por el siguiente código: import matplotlib.pyplot as plt import numpy as np from scipy.integrate import ode from math import sin y0 = [0, 0] t0 = 0 def f(t, y, c, omega0, F0, omega): x, v = y xDot = v vDot = -omega0*omega0 *x - 2*c*v + F0*sin(omega*t) return [xDot, vDot] r = ode(f).set_integrator('dopri5') r.set_initial_value(y0, t0) #(c = 0.5, omega0 = 2.0, F0 = 12.0, omega = 1.0) r.set_f_params(0.5, 2.0, 12.0, 1.0) t1 = 20.0 t = np.linspace(0.0, t1, t1*100 + 1)

8.1. Ecuaciones Diferenciales Ordinarias

69

Programacion con Python, Release 0.4

dt = t[1] y = np.zeros(len(t)) y[0] = y0[0] i = 1 while r.successful() and r.t