8.1 Regresión vs interpolación

Referencia: Chapra 17.1 p 466. Burden 8.1 p498

Cuando los datos de un experimento presentan variaciones o errores sustanciales respecto al modelo matemático, la interpolación polinomial presentada en la Unidad 4 es inapropiada para predecir valores intermedios.

Minimos Cuadrados Lagrange 01

En el ejemplo de Chapra 17.1 p470, se presentan los datos de un experimento mostrados en la imagen y la siguiente tabla:

xi = [1,   2,   3,  4,  5,   6, 7]
yi = [0.5, 2.5, 2., 4., 3.5, 6, 5.5]

Un polinomio de interpolación, por ejemplo de Lagrange de grado 6 pasará por todos los puntos, pero oscilando.

Una función de aproximación que se ajuste a la tendencia general, que no necesariamente pasa por los puntos de muestra puede ser una mejor respuesta. Se busca una "curva" que minimice las diferencias entre los puntos y la curva, llamada regresión por mínimos cuadrados.


Descripción con la media yi

Considere una aproximación para la relación entre los puntos xi, yi como un polinomio, grado 0 que sería la media de yi. Para este caso, los errores se presentan en la gráfica:

Minimos Cuadrados 03

Otra forma sería aproximar el comportamiento de los datos es usar un polinomio de grado 1. En la gráfica se pueden observar que para los mismos puntos el error disminuye considerablemente.

El polinomio de grado 1 recibe el nombre de regresión por mínimos cuadrados, que se desarrolla en la siguiente sección.

Instrucciones Python

# representación con la media
import numpy as np
import matplotlib.pyplot as plt

# INGRESO
xi = [1,   2,   3,  4,  5,   6, 7]
yi = [0.5, 2.5, 2., 4., 3.5, 6, 5.5]

# PROCEDIMIENTO
xi = np.array(xi,dtype=float)
yi = np.array(yi,dtype=float)
n = len(xi)

xm = np.mean(xi)
ym = np.mean(yi)

# SALIDA
print('ymedia = ', ym)

# grafica
plt.plot(xi,yi,'o',label='(xi,yi)')
plt.stem(xi,yi,bottom=ym, linefmt = '--')
plt.xlabel('xi')
plt.ylabel('yi')
plt.legend()
plt.show()

Instrucciones Python - compara interpolación y regresión.

Para ilustrar el asunto y para comparar los resultados se usa Python, tanto para interpolación y mínimos cuadrados usando las librerías disponibles. Luego se desarrolla el algoritmo paso a paso.

# mínimos cuadrados
import numpy as np
import scipy.interpolate as sci
import matplotlib.pyplot as plt

# INGRESO
xi = [1,   2,   3,  4,  5,   6, 7]
yi = [0.5, 2.5, 2., 4., 3.5, 6, 5.5]

# PROCEDIMIENTO
xi = np.array(xi)
yi = np.array(yi)
n = len(xi)

# polinomio Lagrange
px = sci.lagrange(xi,yi)
xj = np.linspace(min(xi),max(xi),100)
pj = px(xj)

# mínimos cuadrados
A = np.vstack([xi, np.ones(n)]).T
[m0, b0] = np.linalg.lstsq(A, yi, rcond=None)[0]
fx = lambda x: m0*(x)+b0
fi = fx(xi)

# ajusta límites
ymin = np.min([np.min(pj),np.min(fi)])
ymax = np.max([np.max(pj),np.max(fi)])

# SALIDA
plt.subplot(121)
plt.plot(xi,yi,'o',label='(xi,yi)')
plt.plot(xj,pj,label='P(x) Lagrange')
plt.ylim(ymin,ymax)
plt.xlabel('xi')
plt.ylabel('yi')
plt.title('Interpolación Lagrange')

plt.subplot(122)
plt.plot(xi,yi,'o',label='(xi,yi)')
plt.plot(xi,fi,label='f(x)')
plt.ylim(ymin,ymax)
plt.xlabel('xi')
plt.title('Mínimos cuadrados')

plt.show()

Medir el tiempo de ejecución en Python

Con varios métodos para resolver un problema, una media de eficiencia es el tiempo de ejecución de un algoritmo.

dibujo cronómetro y banderaPara determinar el tiempo de ejecución de un algoritmo, se toman las lecturas de tiempo antes y después del bloque, calculando el tiempo de ocupación de las diferencias entre tiempos.

La librería "time" permite obtener lectura del reloj del computador.

Como un algoritmo de prueba se usa la suma de los m primeros números enteros.

# tiempos de ejecuci n de algoritmos
import time

# INGRESO
m = 100

# PROCEDIMIENTO
t1 = time.clock()

# Ejecuta Operaciones del algoritmo
suma = 0
for i in range(0,m,1):
    suma = suma + i

t2 = time.clock()
# Tiempo para ejecutar operaciones
ocupado = t2-t1

# SALIDA
print('tiempos:')
print(t1, t2)
print('tiempo de ejecuci n:')
print(ocupado)

con lo que se obtienen los siguientes resultados:

tiempos:
8.210955878428587e-07 5.46028565915501e-05
tiempo de ejecución:
5.378176100370724e-05
>>> 

Sympy - Fórmulas y funciones simbólicas

[ expresiones ] [ operaciones ] [ evaluar f(x) ] [ f(x) a lambda ] [ f(x) matemáticas ] [ derivadas ] [ integrales ]

SymPy es una librería usada para manejar de forma algebraica las expresiones matemáticas. Se requiere definir el símbolo que representa la variable en la expresión, por ejemplo la letra 'x'. Si la librería Sympy no está disponible o muestra un error la puede revisar las instrucciones del enlace instalar con pip.

Una formula o función matemática f(x) descrita como fx se puede derivar, integrar, simplificar sus términos. También se puede construir una expresión matemática, por ejemplo desarrollar los términos para un polinomio como Taylor.

Referencia: https://www.sympy.org/es/

[ expresiones ] [ operaciones ] [ evaluar f(x) ] [ f(x) a lambda ] [ f(x) matemáticas ] [ derivadas ] [ integrales ]
..


1. Expresiones con Sympy

Una función f(x) como un polinomio

f(x) = (1-x)^5 + 5 x ((1-x)^4) - 0.4

Se pueden manejar de forma algebraica con Sympy al declarar la variable independiente 'x' como un símbolo. La expresión se asigna a fx

import sympy as sym
x = sym.Symbol('x')
fx = (1-x)**5 + 5*x*((1-x)**4) - 0.4

Se simplifican o se expande los términos del polinomio con solo una instrucción,

fx = fx.expand()
print(fx)

mostrando el siguiente resultado:

4*x**5 - 15*x**4 + 20*x**3 - 10*x**2 + 0.6

o una presentación diferente con sym.pprint():

>>> sym.pprint(fx)
           4          5      
5*x*(1 - x)  + (1 - x)  - 0.4
>>> 

Referencia: https://docs.sympy.org/latest/tutorials/intro-tutorial/printing.html#setting-up-pretty-printing

[ expresiones ] [ operaciones ] [ evaluar f(x) ] [ f(x) a lambda ] [ f(x) matemáticas ] [ derivadas ] [ integrales ]
..


2. Operaciones, combinar otras expresiones

A una función simbólica fx se pueden añadir más términos con el mismo símbolo

>>> import sympy as sym
>>> sym.Symbol('x')
>>> fx = sym.cos(x)
>>> gx = fx + x
>>> gx
x + cos(x)
>>> gx = gx - 2
>>> gx
x + cos(x) - 2

Por lo que las funciones simbólicas son útiles cuando se construyen expresiones, como en el caso de series de Taylor descritas en la Unidad01

[ expresiones ] [ operaciones ] [ evaluar f(x) ] [ f(x) a lambda ] [ f(x) matemáticas ] [ derivadas ] [ integrales ]
..


3. Evaluar una expresión f(x) con Sympy

Las expresiones simbólicas se pueden evaluar en un punto, ejemplo x0=0.1 usando la instrucción fx.subs(x,x0), que sustituye el símbolo de la variable x con el valor definido para x0.

>>> import sympy as sym
>>> sym.Symbol('x')
>>> fx = sym.cos(x)
>>> fx.subs(x,0.1)
0.995004165278026

Referencia: https://docs.sympy.org/latest/tutorials/intro-tutorial/basic_operations.html#substitution

[ expresiones ] [ operaciones ] [ evaluar f(x) ] [ f(x) a lambda ] [ f(x) matemáticas ] [ derivadas ] [ integrales ]
..


4. Conversión de forma simbólica a forma numérica Lambda

Para evaluar varios puntos en la expresión  en forma numérica para usar 'x' con valores en un vector o matriz como un arreglo, se convierte a la forma numérica Lambda con la instrucción sym.lambdify(x,fx)

>>> >>> import sympy as sym
>>> sym.Symbol('x')
>>> fx = sym.cos(x)
&>>> f_x = sym.lambdify(x,fx)
>>> f_x(0.1)
0.9950041652780258
>>> f_x([0.1,0.3,0.5])
array([0.99500417, 0.95533649, 0.87758256])
>>> 

Referencia: https://docs.sympy.org/latest/tutorials/intro-tutorial/basic_operations.html#lambdify

Ejemplo: Polinomio de Taylor – Ejemplos con Sympy-Python

[ expresiones ] [ operaciones ] [ evaluar f(x) ] [ f(x) a lambda ] [ f(x) matemáticas ] [ derivadas ] [ integrales ]
..


5. expresiones de Funciones matemáticas

Las Funciones matemáticas usadas en otras librerías como Numpy tienen también su representación simbólica en Sympy. Ejemplo cos(x):

f(x) = \cos (x)
>>> import sympy as sym
>>> x = sym.Symbol('x')
>>> fx = sym.cos(x)

Realice pruebas con otras funciones: sin(x), exp(x), log(x), log10(x), Heaviside(x)

[ expresiones ] [ operaciones ] [ evaluar f(x) ] [ f(x) a lambda ] [ f(x) matemáticas ] [ derivadas ] [ integrales ]
..


6. Derivadas con Sympy

Las expresiones de la derivada se obtienen con la expresión fx.diff(x,k), indicando la k-ésima  derivada de la expresión.

>>> import sympy as sym
>>> x = sym.Symbol('x')
>>> fx = sym.cos(x)
>>> x
x
>>> fx
cos(x)
>>> fx.diff(x,1)
-sin(x)
>>> fx.diff(x,2)
-cos(x)

Referencia: https://docs.sympy.org/latest/tutorials/intro-tutorial/calculus.html#derivatives

Ejemplos: Sistemas LTI CT – Respuesta a entrada cero con Sympy-Python

Derivadas como expresión de Sympy sin evaluar

Cuando se requiere expresar tan solo la operación de derivadas, que será luego usada o reemplazada con otra expresión, se usa la derivada sin evaluar. Ejemplo:

y = \frac{d}{dx}f(x)

Para más adelante definir f(x).

En Sympy, la expresión de y se realiza indicando que f será una variable

x = sym.Symbol('x', real=True)
f = sym.Symbol('f', real=True)
y = sym.diff(f,x, evaluate=False) # derivada sin evaluar
g = sym.cos(x) + x**2
yg = y.subs(f,g).doit() # sustituye f con g y evalua .doit()
>>> y
Derivative(f, x)
>>> g
x**2 + cos(x)
>>> yg
2*x - sin(x)

Ejemplos:

Polinomio de Taylor – Ejemplos con Sympy-Python

EDO lineal – solución complementaria y particular con Sympy-Python

EDO lineal – ecuaciones auxiliar, general y complementaria con Sympy-Python

EDP Parabólica – analítico con Sympy-Python

EDP Elípticas – analítico con Sympy-Python

Evaluación de las expresiones Sympy

Se define la expresión 'derivada' y se la usa con la instrucción fx.subs(x,valor)

>>> fx.subs(x,0)
1
>>> fx.subs(x,1/3)
0.944956946314738
>>> derivada = fx.diff(x,1)
>>> derivada
-sin(x)
>>> derivada.subs(x,1/3)
-0.327194696796152
>>> 

[ expresiones ] [ operaciones ] [ evaluar f(x) ] [ f(x) a lambda ] [ f(x) matemáticas ] [ derivadas ] [ integrales ]
..


7. Integrales con Sympy

Sympy incorpora la operación del integral en sus expresiones, que pueden ser integrales definidas en un intervalo o expresiones sin evaluar.

>>> import sympy as sym
>>> t = sym.Symbol('t',real=True)
>>> x = 10*sym.exp(-3*t)
>>> x
10*exp(-3*t)
>>> y = sym.integrate(x,(t,0,10))
>>> y
10/3 - 10*exp(-30)/3
>>> y = sym.integrate(x,(t,0,sym.oo))
>>> y
10/3

Referencia: https://docs.sympy.org/latest/tutorials/intro-tutorial/calculus.html#integrals

Ejemplos:
Transformada de Laplace – Integral con Sympy-Python

Series de Fourier de señales periódicas con Sympy-Python

Integral de Convolución entre x(t) y h(t) causal con Sympy-Python

[ expresiones ] [ operaciones ] [ evaluar f(x) ] [ f(x) a lambda ] [ f(x) matemáticas ] [ derivadas ] [ integrales ]

Numpy - Vectores y matrices

Resumen de algunas funciones para vectores y matrices como arreglos en la librería Numpy de Python, usadas en el curso para facilitan el cálculo numérico. El orden de las instrucciones es el que aparece en las entradas del blog.

Lo primero es hacer el llamado a las librerías con el alias 'np'

import numpy as np

Vectores

puntos espaciados entre [a,b] - np.linspace()

para obtener puntos en el rango [a,b] para una cantidad de tramos. Se obtienen tramos+1 puntos .

a = 1
b = 2
tramos = 4
x = np.linspace(a,b,tramos+1)
>>> x
array([ 1.  ,  1.25,  1.5 ,  1.75,  2.  ])

rango de valores en vector - np.arange(a,b,dt)

crea un vector con valores en el rango [a,b) y espaciados dt.

t=np.arange(0,10,2)
>>> t
array([0, 2, 4, 6, 8])

transponer vector - np.transpose()

>>> fila = np.array([1,2,3])
>>> columna = np.transpose([fila])
>>> columna
array([[1],
       [2],
       [3]])

indices para ordenar - np.argsort()

para ordenar por una columna específica:
- se obtiene un vector con la columna que se requiere ordenar: tabla[:,0]
- con el vector se determinan los índices para ordenar la tabla por filas: np.argsort()
- se aplica los índices a una copia de la matriz: tabla[orden]

tabla = np.array([[5,3],
                  [4,2],
                  [3,1],
                  [2,4],
                  [1,5]])

# ordenar por primera columna
referencia = tabla[:,0]
orden = np.argsort(referencia)
ordenada = tabla[orden]

print(orden)
print(ordenada)

resultado:

[4 3 2 1 0]
[[1 5]
 [2 4]
 [3 1]
 [4 2]
 [5 3]]
# ordenar por segunda columna
referencia = tabla[:,1]
orden = np.argsort(referencia)
ordenada = tabla[orden]
print(orden)
print(ordenada)

resultado:

[2 1 0 3 4]
[[3 1]
 [4 2]
 [5 3]
 [2 4]
 [1 5]]

Matrices en Numpy

concatenar- np.concatenate((A,b), axis=1)

concatenar usa el parámetro axis: 0 para filas, 1 para columnas.

import numpy as np
cantidad = np.array([[4,2,5],
                     [2,5,8],
                     [5,4,3]])

pagado = np.array([[60.70],
                   [92.90],
                   [56.30]])

matriz = np.concatenate((cantidad,pagado),axis=1)
>>> matriz
array([[  4. ,   2. ,   5. ,  60.7],
       [  2. ,   5. ,   8. ,  92.9],
       [  5. ,   4. ,   3. ,  56.3]])

resolver sistema de ecuaciones - np.linalg.solve(A,B)

Resuelve el sistema de ecuaciones dado por una matriz A y un vector B. siendo, por ejemplo:

 0 =  c1 +  c2
-5 = -c1 - 2c2

c1 = -5 
c2 = 5
>>> A = [[ 1, 1],
	 [-1,-2]]
>>> B = [0,-5]
>>> np.linalg.solve(A,B)
array([-5.,  5.])
>>>

Otras instrucciones

np.pi constante con valor π

>>> np.pi
3.141592653589793
np.sin(t)
np.cos(t)
función trigonométrica en radianes. La variable t puede ser un escalar o un arreglo.

>>> t=0.65
>>> np.sin(0.65)
0.60518640573603955
>>> t=[0, 0.3, 0.6]
>>> np.sin(t)
array([ 0. , 0.29552021, 0.56464247])
np.abs() obtiene el valor absoluto de un número. En el caso de un número complejo obtiene la parte real.
np.real(complejo)
np.imag(complejo)
obtiene la parte real de los números complejos en un vector. Se aplica lo mismo para la parte imaginaria del número complejo.
complex(a,b) crea el número complejo a partir de los valores de a y b.
a=2
b=3
el resultado es: 2+3j
np.piecewise(t, t>=donde, [1,0]) función que crea a partir de t, los valores de la condición t>=donde, ubicando los valores de 1, para otro caso es 0. Usada en la función escalón.
np.roots([a,b,c]) obtiene las raíces del polinomio:
ax2+bx+c
x2 + 3 x + 2 = (x+1)(x+2)

>>> np.roots([1,3,2])
array([-2., -1.])

Funciones lambda vs def-return

Las dos formas de escritura de funciones matemáticas básicamente hacen lo mismo. Sin embargo, cuando la función se define por tramos, la forma def-return se convierte en la más usada.

Se adjunta un ejemplo:

Funciones Lambda

Permite describir funciones sencillas de una sola línea, que no está conformada por intervalos.

f(x) = x \cos(x) - 2 x^2 + 3 x -1
import numpy as np
fx = lambda x: x*np.cos(x) - 2*(x**2) + 3*x -1
>>> fx(0.2)
-0.28398668443175157

Funciones def-return

Permiten describir en detalle lo que sucede con una función matemática por intervalos. Tiene la ventaja de permitir definir valores por intervalos, puntos discontínuos, etc.

f(x)= \begin{cases} x \cos(x) - 2 x^2 + 3 x -1, & x>0 \\1, & x \leq 0 \end{cases}
import numpy as np
def fx(x):
    if x>0:
        fi = x*np.cos(x) - 2*(x**2) + 3*x -1
    if x<=0:
        fi = 1 
    return(fi)
>>> fx(0.2)
 -0.28398668443175157

7.3 EDP hiperbólicas con Python

[ concepto ] [ analítico ] [ algoritmo ]
..


1. Concepto y ejercicio

Referencia:  Chapra PT8.1 p860,  Rodríguez 10.4 p435

Las Ecuaciones Diferenciales Parciales tipo hiperbólicas semejantes a la mostrada, para un ejemplo en particular, representa la ecuación de onda de una dimensión u[x,t], que representa el desplazamiento vertical de una cuerda.

\frac{\partial ^2 u}{\partial t^2}=c^2\frac{\partial ^2 u}{\partial x^2}

EDP Cuerda 01

Los extremos de la cuerda de longitud unitaria están sujetos y referenciados a una posición 0 a la izquierda y 1 a la derecha.

u[x,t] , 0<x<1, t≥0
u(0,t) = 0 , t≥0
u(1,t) = 0 , t≥0

Al inicio, la cuerda está estirada por su punto central:

u(x,0) = \begin{cases} -0.5x &, 0\lt x\leq 0.5 \\ 0.5(x-1) &, 0.5\lt x \lt 1 \end{cases}

Se suelta la cuerda, con velocidad cero desde la posición inicial:

\frac{\partial u(x,0)}{\partial t} = 0

[ concepto ] [ analítico ] [ algoritmo ]
..


2. Desarrollo analítico

La solución se realiza de forma semejante al procedimiento para EDP parabólicas y elípticas. Se cambia a la forma discreta  usando diferencias finitas divididas:

\frac{u_{i,j+1}-2u_{i,j}+u_{i,j-1}}{(\Delta t)^2} =c^2 \frac{u_{i+1,j}-2u_{i,j}+u_{i-1,j}}{(\Delta x)^2}

El error es del orden O(\Delta x)^2 + O(\Delta t)^2.
se reagrupa de la forma:

u_{i,j+1}-2u_{i,j}+u_{i,j-1} = \frac{c^2 (\Delta t)^2}{(\Delta x)^2} \big( u_{i+1,j}-2u_{i,j}+u_{i-1,j} \big)

El término constante, por facilidad se simplifica con el valor de 1

\lambda = \frac{c^2 (\Delta t)^2}{(\Delta x)^2} =1

si c = 2 y Δx = 0.2, se deduce que Δt = 0.1

que al sustituir en la ecuación, se simplifica anulando el término u[i,j]:

u_{i,j+1}+u_{i,j-1} = u_{i+1,j}+u_{i-1,j}

EDP Cuerda 02

en los que intervienen solo los puntos extremos. Despejando el punto superior del rombo:

u_{i,j+1} = u_{i+1,j}-u_{i,j-1}+u_{i-1,j}

Convergencia:

\lambda = \frac{c^2 (\Delta t)^2}{(\Delta x)^2} \leq 1

para simplificar aún más el problema, se usa la condición velocidad inicial de la cuerda igual a cero

\frac{\delta u_{i,0}}{\delta t}=\frac{u_{i,1}-u_{i,-1}}{2\Delta t} = 0 u_{i,-1}=u_{i,1}

que se usa para cuando el tiempo es cero, permite calcular los puntos para t[1]:

j=0

u_{i,1} = u_{i+1,0}-u_{i,-1}+u_{i-1,0} u_{i,1} = u_{i+1,0}-u_{i,1}+u_{i-1,0} 2 u_{i,1} = u_{i+1,0}+u_{i-1,0} u_{i,1} = \frac{u_{i+1,0}+u_{i-1,0}}{2}

Aplicando solo cuando j = 0

EDP Cuerda 03

que al ponerlos en la malla se logra un sistema explícito para cada u[i,j], con lo que se puede resolver con un algoritmo:

EDP hiperbolica 01

[ concepto ] [ analítico ] [ algoritmo ]
..


3. Algoritmo con Python

# Ecuaciones Diferenciales Parciales
# Hiperbólica. Método explicito
import numpy as np

def cuerdainicio(xi):
    n = len(xi)
    y = np.zeros(n,dtype=float)
    for i in range(0,n,1):
        if (xi[i]<=0.5):
            y[i]=-0.5*xi[i]
        else:
            y[i]=0.5*(xi[i]-1)
    return(y)

# INGRESO
x0 = 0 # Longitud de cuerda
xn = 1
y0 = 0 # Puntos de amarre
yn = 0
t0 = 0 # tiempo inicial
c = 2  # constante EDP
# discretiza
tramosx = 16
tramost = 32
dx = (xn-x0)/tramosx 
dt = dx/c

# PROCEDIMIENTO
xi = np.arange(x0,xn+dx/2,dx)
tj = np.arange(0,tramost*dt+dt/2,dt)
n = len(xi)
m = len(tj)

u = np.zeros(shape=(n,m),dtype=float)
u[:,0] = cuerdainicio(xi)
u[0,:] = y0
u[n-1,:] = yn
# Aplicando condición inicial
j = 0
for i in range(1,n-1,1):
    u[i,j+1] = (u[i+1,j]+u[i-1,j])/2
# Para los otros puntos
for j in range(1,m-1,1):
    for i in range(1,n-1,1):
        u[i,j+1] = u[i+1,j]-u[i,j-1]+u[i-1,j]

# SALIDA
np.set_printoptions(precision=2)
print('xi =')
print(xi)
print('tj =')
print(tj)
print('matriz u =')
print(u)

con lo que se obtienen los resultados numéricos, que para mejor interpretación se presentan en una gráfica estática y otra animada.

# GRAFICA
import matplotlib.pyplot as plt

for j in range(0,m,1):
    y = u[:,j]
    plt.plot(xi,y)

plt.title('EDP hiperbólica')
plt.xlabel('x')
plt.ylabel('y')
plt.show()

# **** GRÁFICO CON ANIMACION ***********
import matplotlib.animation as animation

# Inicializa parametros de trama/foto
retardo = 70   # milisegundos entre tramas
tramas  = m
maximoy = np.max(np.abs(u))
figura, ejes = plt.subplots()
plt.xlim([x0,xn])
plt.ylim([-maximoy,maximoy])

# lineas de función y polinomio en gráfico
linea_poli, = ejes.plot(xi,u[:,0], '-')
plt.axhline(0, color='k')  # Eje en x=0

plt.title('EDP hiperbólica')
# plt.legend()
# txt_x = (x0+xn)/2
txt_y = maximoy*(1-0.1)
texto = ejes.text(x0,txt_y,'tiempo:',
                  horizontalalignment='left')
plt.xlabel('x')
plt.ylabel('y')
plt.grid()

# Nueva Trama
def unatrama(i,xi,u):
    # actualiza cada linea
    linea_poli.set_ydata(u[:,i])
    linea_poli.set_xdata(xi)
    linea_poli.set_label('tiempo linea: '+str(i))
    texto.set_text('tiempo['+str(i)+']')
    # color de la línea
    if (i<=9):
        lineacolor = 'C'+str(i)
    else:
        numcolor = i%10
        lineacolor = 'C'+str(numcolor)
    linea_poli.set_color(lineacolor)
    return linea_poli, texto

# Limpia Trama anterior
def limpiatrama():
    linea_poli.set_ydata(np.ma.array(xi, mask=True))
    linea_poli.set_label('')
    texto.set_text('')
    return linea_poli, texto

# Trama contador
i = np.arange(0,tramas,1)
ani = animation.FuncAnimation(figura,
                              unatrama,
                              i ,
                              fargs=(xi,u),
                              init_func=limpiatrama,
                              interval=retardo,
                              blit=True)
# Graba Archivo video y GIFAnimado :
# ani.save('EDP_hiperbólica.mp4')
ani.save('EDP_hiperbolica.gif', writer='imagemagick')
plt.draw()
plt.show()

Una vez realizado el algoritmo, se pueden cambiar las condiciones iniciales de la cuerda y se observan los resultados.

Se recomienda realizar otros ejercicios de la sección de evaluaciones para EDP Hiperbólicas y observar los resultados con el algoritmo modificado.

[ concepto ] [ analítico ] [ algoritmo ]

7.2.2 EDP Elípticas método implícito con Python

EDP Elípticas [ concepto ] Método implícito: [ Analítico ] [ Algoritmo ]
..


1. EDP Elípticas: Método Implícito – Desarrollo Analítico

con el resultado desarrollado en EDP elípticas para:

\frac{\partial ^2 u}{\partial x^2} + \frac{\partial ^2 u}{ \partial y^2} = 0

y con el supuesto que: \lambda = \frac{(\Delta y)^2}{(\Delta x)^2} = 1

se puede plantear que:

u_{i+1,j}-4u_{i,j}+u_{i-1,j} + u_{i,j+1} +u_{i,j-1} = 0

EDP Elipticas Iterativo

con lo que para el método implícito, se plantea un sistema de ecuaciones para determinar los valores en cada punto desconocido.

j=1, i =1

u_{2,1}-4u_{1,1}+u_{0,1} + u_{1,2} +u_{1,0} = 0 u_{2,1}-4u_{1,1}+Ta + u_{1,2} +Tc= 0 -4u_{1,1}+u_{2,1}+u_{1,2} = -(Tc+Ta)

j=1, i =2

u_{3,1}-4u_{2,1}+u_{1,1} + u_{2,2} +u_{2,0} = 0 u_{3,1}-4u_{2,1}+u_{1,1} + u_{2,2} +Tc = 0 u_{1,1}-4u_{2,1}+u_{3,1}+ u_{2,2}= -Tc

j=1, i=3

u_{4,1}-4u_{3,1}+u_{2,1} + u_{3,2} +u_{3,0} = 0 Tb-4u_{3,1}+u_{2,1} + u_{3,2} +Tc = 0 u_{2,1} -4u_{3,1} + u_{3,2} = -(Tc+Tb)

j=2, i=1

u_{2,2}-4u_{1,2}+u_{0,2} + u_{1,3} +u_{1,1} = 0 u_{2,2}-4u_{1,2}+Ta + u_{1,3} +u_{1,1} = 0 -4u_{1,2}+u_{2,2}+u_{1,1}+u_{1,3} = -Ta

j = 2, i = 2

u_{1,2}-4u_{2,2}+u_{3,2} + u_{2,3} +u_{2,1} = 0

j = 2, i = 3

u_{4,2}-4u_{3,2}+u_{2,2} + u_{3,3} +u_{3,1} = 0 Tb-4u_{3,2}+u_{2,2} + u_{3,3} +u_{3,1} = 0 u_{2,2} -4u_{3,2}+ u_{3,3} +u_{3,1} = -Tb

j=3, i = 1

u_{2,3}-4u_{1,3}+u_{0,3} + u_{1,4} +u_{1,2} = 0 u_{2,3}-4u_{1,3}+Ta + Td +u_{1,2} = 0 -4u_{1,3}+u_{2,3}+u_{1,2} = -(Td+Ta)

j=3, i = 2

u_{3,3}-4u_{2,3}+u_{1,3} + u_{2,4} +u_{2,2} = 0 u_{3,3}-4u_{2,3}+u_{1,3} + Td +u_{2,2} = 0 +u_{1,3} -4u_{2,3}+u_{3,3} +u_{2,2} = -Td

j=3, i=3

u_{4,3}-4u_{3,3}+u_{2,3} + u_{3,4} +u_{3,2} = 0 Tb-4u_{3,3}+u_{2,3} + Td +u_{3,2} = 0 u_{2,3}-4u_{3,3}+u_{3,2} = -(Td+Tb)

con las ecuaciones se arma una matriz:

A = np.array([
    [-4, 1, 0, 1, 0, 0, 0, 0, 0],
    [ 1,-4, 1, 0, 1, 0, 0, 0, 0],
    [ 0, 1,-4, 0, 0, 1, 0, 0, 0],
    [ 1, 0, 0,-4, 1, 0, 1, 0, 0],
    [ 0, 1, 0, 1,-4, 1, 0, 1, 0],
    [ 0, 0, 1, 0, 1,-4, 0, 0, 1],
    [ 0, 0, 0, 1, 0, 0,-4, 1, 0],
    [ 0, 0, 0, 0, 1, 0, 1,-4, 1],
    [ 0, 0, 0, 0, 0, 1, 0, 1,-4],
    ])
B = np.array([-(Tc+Ta),-Tc,-(Tc+Tb),
                  -Ta,   0,    -Tb,
              -(Td+Ta),-Td,-(Td+Tb)])

que al resolver el sistema de ecuaciones se obtiene:

>>> Xu
array([ 56.43,  55.71,  56.43,  60.  ,  60.  ,  60.  ,  63.57,  64.29,
        63.57])

ingresando los resultados a la matriz u:

xi=
[ 0.   0.5  1.   1.5  2. ]
yj=
[ 0.    0.38  0.75  1.12  1.5 ]
matriz u
[[ 60.    60.    60.    60.    60.  ]
 [ 50.    56.43  60.    63.57  70.  ]
 [ 50.    55.71  60.    64.29  70.  ]
 [ 50.    56.43  60.    63.57  70.  ]
 [ 60.    60.    60.    60.    60.  ]]
>>>

EDP Elípticas [ concepto ] Método implícito: [ Analítico ] [ Algoritmo ]

..


2. EDP Elípticas: Método Implícito – Desarrollo con algoritmo

Instrucciones en Python

# Ecuaciones Diferenciales Parciales
# Elipticas. Método implícito
import numpy as np

# INGRESO
# Condiciones iniciales en los bordes
Ta = 60
Tb = 60
Tc = 50
Td = 70
# dimensiones de la placa
x0 = 0
xn = 2
y0 = 0
yn = 1.5
# discretiza, supone dx=dy
tramosx = 4
tramosy = 4
dx = (xn-x0)/tramosx 
dy = (yn-y0)/tramosy 
maxitera = 100
tolera = 0.0001

A = np.array([
    [-4, 1, 0, 1, 0, 0, 0, 0, 0],
    [ 1,-4, 1, 0, 1, 0, 0, 0, 0],
    [ 0, 1,-4, 0, 0, 1, 0, 0, 0],
    [ 1, 0, 0,-4, 1, 0, 1, 0, 0],
    [ 0, 1, 0, 1,-4, 1, 0, 1, 0],
    [ 0, 0, 1, 0, 1,-4, 0, 0, 1],
    [ 0, 0, 0, 1, 0, 0,-4, 1, 0],
    [ 0, 0, 0, 0, 1, 0, 1,-4, 1],
    [ 0, 0, 0, 0, 0, 1, 0, 1,-4],
    ])
B = np.array([-(Tc+Ta),-Tc,-(Tc+Tb),
              -Ta,0,-Tb,
              -(Td+Ta),-Td,-(Td+Tb)])


# PROCEDIMIENTO
# Resuelve sistema ecuaciones
Xu = np.linalg.solve(A,B)
[nx,mx] = np.shape(A)

xi = np.arange(x0,xn+dx/2,dx)
yj = np.arange(y0,yn+dy/2,dy)
n = len(xi)
m = len(yj)

u = np.zeros(shape=(n,m),dtype=float)
u[:,0]   = Tc
u[:,m-1] = Td
u[0,:]   = Ta
u[n-1,:] = Tb
u[1:1+3,1] = Xu[0:0+3]
u[1:1+3,2] = Xu[3:3+3]
u[1:1+3,3] = Xu[6:6+3]

# SALIDA
np.set_printoptions(precision=2)
print('xi=')
print(xi)
print('yj=')
print(yj)
print('matriz u')
print(u)

La gráfica de resultados se obtiene de forma semejante al ejercicio con método iterativo.

Se podría estandarizar un poco más el proceso para que sea realizado por el algoritmo y sea más sencillo generar la matriz con más puntos. Tarea.

7.2.1 EDP Elípticas método iterativo con Python

EDP Elípticas [ concepto ] [ ejercicio ]
Método iterativo: [ Analítico ] [ Algoritmo ]

..


1. Ejercicio

Referencia: Chapra 29.1 p866, Rodríguez 10.3 p425, Burden 12.1 p694

Siguiendo el tema anterior en EDP Elípticas, para resolver la parte numérica asuma como:

Valores de frontera: Ta = 60 , Tb = 25 , Tc = 50, Td = 70
Longitud en:  x0 = 0, xn = 2 , y0 = 0, yn = 1.5
Tamaño de paso dx = 0.25 dy = dx
iteramax = 100 tolera = 0.0001

Para la ecuación diferencial parcial Elíptica

\frac{\partial ^2 u}{\partial x^2} + \frac{\partial ^2 u}{ \partial y^2} = 0


EDP Elípticas [ concepto ] [ ejercicio ]
Método iterativo: [ Analítico ] [ Algoritmo ]
..


2. EDP Elípticas: Método iterativo - Desarrollo Analítico

A partir del resultado desarrollado en EDP elípticas:

\frac{(\Delta y)^2}{(\Delta x)^2} \Big( u_{i+1,j}-2u_{i,j} +u_{i-1,j} \Big)+ u_{i,j+1}-2u_{i,j}+u_{i,j-1}=0 \lambda \Big( u_{i+1,j}-2u_{i,j} +u_{i-1,j} \Big)+ u_{i,j+1}-2u_{i,j}+u_{i,j-1}=0 \lambda u_{i+1,j} +(-2\lambda-2) u_{i,j} +\lambda u_{i-1,j} + u_{i,j+1} +u_{i,j-1}=0

y con el supuesto que: \lambda = \frac{(\Delta y)^2}{(\Delta x)^2} = 1

se puede plantear que:

u_{i+1,j}-4u_{i,j}+u_{i-1,j} + u_{i,j+1} +u_{i,j-1} = 0

que reordenando para un punto central desconocido se convierte a:

u_{i,j} = \frac{1}{4} \big[ u_{i+1,j}+u_{i-1,j} + u_{i,j+1} +u_{i,j-1} \big]

con lo que se interpreta que cada punto central es el resultado del promedio de los puntos alrededor del rombo formado en la malla.

El cálculo numérico se puede realizar de forma iterativa haciendo varias pasadas en la matriz, promediando cada punto. Para revisar las iteraciones se controla la convergencia junto con un máximo de iteraciones.

EDP Elípticas [ concepto ] [ ejercicio ]
Método iterativo: [ Analítico ] [ Algoritmo ]

..


3. Algoritmo en Python

Para el ejercicio dado, se presenta el resultado para tres iteraciones y el resultado final con un gráfico en 3D:

u inicial:
[[50.   50.   50.   50.   50.   50.   50.   50.   50.  ]
 [60.   51.25 51.25 51.25 51.25 51.25 51.25 51.25 25.  ]
 [60.   51.25 51.25 51.25 51.25 51.25 51.25 51.25 25.  ]
 [60.   51.25 51.25 51.25 51.25 51.25 51.25 51.25 25.  ]
 [60.   51.25 51.25 51.25 51.25 51.25 51.25 51.25 25.  ]
 [60.   51.25 51.25 51.25 51.25 51.25 51.25 51.25 25.  ]
 [70.   70.   70.   70.   70.   70.   70.   70.   70.  ]]
itera: 0   ;  u:
[[50.   50.   50.   50.   50.   50.   50.   50.   50.  ]
 [60.   53.12 51.41 50.98 50.87 50.84 50.84 44.27 25.  ]
 [60.   53.91 51.95 51.36 51.18 51.13 51.12 42.91 25.  ]
 [60.   54.1  52.14 51.5  51.3  51.23 51.21 42.59 25.  ]
 [60.   54.15 52.2  51.55 51.34 51.27 51.24 42.52 25.  ]
 [60.   58.85 58.07 57.72 57.58 57.52 57.5  48.76 25.  ]
 [70.   70.   70.   70.   70.   70.   70.   70.   70.  ]]
errado_u:  8.728094100952148
itera: 1   ;  u:
[[50.   50.   50.   50.   50.   50.   50.   50.   50.  ]
 [60.   53.83 51.69 50.98 50.75 50.68 49.02 41.73 25.  ]
 [60.   54.97 52.54 51.55 51.18 51.05 48.55 39.47 25.  ]
 [60.   55.31 52.89 51.82 51.39 51.23 48.4  38.85 25.  ]
 [60.   56.59 54.78 53.91 53.54 53.38 50.45 40.76 25.  ]
 [60.   61.17 60.91 60.6  60.42 60.33 57.38 48.29 25.  ]
 [70.   70.   70.   70.   70.   70.   70.   70.   70.  ]]
errado_u:  3.7443900108337402
itera: 2   ;  u:
[[50.   50.   50.   50.   50.   50.   50.   50.   50.  ]
 [60.   54.17 51.92 51.06 50.73 50.2  47.62 40.52 25.  ]
 [60.   55.5  52.97 51.76 51.23 50.3  46.45 37.7  25.  ]
 [60.   56.25 53.95 52.75 52.19 51.07 46.71 37.54 25.  ]
 [60.   58.05 56.71 55.9  55.47 54.33 49.8  40.16 25.  ]
 [60.   62.24 62.39 62.18 61.99 60.93 57.25 48.1  25.  ]
 [70.   70.   70.   70.   70.   70.   70.   70.   70.  ]]
errado_u:  2.0990657806396484
... continua
Método Iterativo EDP Elíptica
iteraciones: 41  converge =  True
xi: [0.   0.25 0.5  0.75 1.   1.25 1.5  1.75 2.  ]
yj: [0.   0.25 0.5  0.75 1.   1.25 1.5 ]
Tabla de resultados en malla EDP Elíptica
j, U[i,j]
6 [70. 70. 70. 70. 70. 70. 70. 70. 70.]
5 [60.   64.02 64.97 64.71 63.62 61.44 57.16 47.96 25.  ]
4 [60.   61.1  61.14 60.25 58.35 54.98 49.23 39.67 25.  ]
3 [60.   59.23 58.25 56.8  54.53 50.89 45.13 36.48 25.  ]
2 [60.   57.56 55.82 54.19 52.09 48.92 43.91 36.14 25.  ]
1 [60.   55.21 53.27 52.05 50.73 48.78 45.46 39.15 25.  ]
0 [50. 50. 50. 50. 50. 50. 50. 50. 50.]
>>>

cuyos valores se interpretan mejor en una gráfica, en este caso 3D:

EDP Elipticas Iterativo

Instrucciones en Python

# Ecuaciones Diferenciales Parciales
# Elipticas. Método iterativo
# d2u/dx2  + du/dt = f(x,y)
import numpy as np

# INGRESO
fxy = lambda x,y: 0*x+0*y # f(x,y) = 0 , ecuacion de Poisson
# Condiciones iniciales en los bordes, valores de frontera
Ta = 60     # izquierda de la placa
Tb = 25     # derecha de la placa
#Tc = 50    # inferior 
fxc = lambda x: 50+0*x #  f(x) inicial inferior
# Td = 70   # superior
fxd = lambda x: 70+0*x #  f(x) inicial superior

# dimensiones de la placa
x0 = 0    # longitud en x
xn = 2
y0 = 0    # longitud en y
yn = 1.5
# discretiza, supone dx=dy
dx = 0.25  # Tamaño de paso
dy = dx

iteramax = 100  # maximo de iteraciones
tolera = 0.0001
verdigitos = 2   # decimales a mostrar en tabla de resultados
vertabla = True  # ver iteraciones

# coeficientes de U[x,t]. factores P,Q,R
lamb = (dy**2)/(dx**2)
P = lamb       # izquierda P*U[i-1,j]
Q = -2*lamb-2  # centro Q*U[i,j]
R = lamb       # derecha R*U[i+1,j]

# PROCEDIMIENTO
# iteraciones en xi,yj ancho,profundidad
xi = np.arange(x0,xn+dx/2,dx)
yj = np.arange(y0,yn+dy/2,dy)
n = len(xi)
m = len(yj)

# Matriz u[xi,yj], tabla de resultados
u = np.zeros(shape=(n,m),dtype=float)

# llena u con valores en fronteras
u[0,:]   = Ta   # izquierda
u[n-1,:] = Tb   # derecha
fic = fxc(xi)   # Tc inferior
u[:,0]   = fic
fid = fxd(xi)   # Td superior
u[:,m-1] = fid  
# estado inicial dentro de la placa
# promedio = (Ta+Tb+Tc+Td)/4
promedio = (Ta+Tb+np.max(fic)+np.max(fid))*(1/-Q)
u[1:n-1,1:m-1] = promedio

np.set_printoptions(precision=verdigitos)
if vertabla == True:
    print('u inicial:')
    print(np.transpose(u))

# iterar puntos interiores
U_3D = [np.copy(u)]
U_error = ['--']
itera = 0
converge = False
while (itera<=iteramax and converge==False):
    itera = itera +1
    Uantes = np.copy(u)
    for i in range(1,n-1):
        for j in range(1,m-1):
            u[i,j] = P*u[i-1,j]+R*u[i+1,j]+u[i,j-1]+u[i,j+1]
            u[i,j] = u[i,j] - fxy(xi[i],yj[j])*(dy**2)
            u[i,j] = (1/-Q)*u[i,j]
    diferencia = u - Uantes
    errado_u = np.max(np.abs(diferencia))
    if (errado_u<tolera):
        converge=True
    U_3D.append(np.copy(u))
    U_error.append(np.around(errado_u,verdigitos))
    
    if (vertabla==True) and (itera<=3):
        np.set_printoptions(precision=verdigitos)
        print('itera:',itera-1,'  ; ','u:')
        print(np.transpose(u))
        print('errado_u: ', errado_u)
    if (vertabla==True) and (itera==4):
        print('... continua')
        
# SALIDA
print('Método Iterativo EDP Elíptica')
print('iteraciones:',itera,' converge = ', converge)
print('errado_u: ', errado_u)
print('xi:', xi)
print('yj:', yj)
print('Tabla de resultados en malla EDP Elíptica')
print('j, U[i,j]')
for j in range(m-1,-1,-1):
    print(j,u[:,j])

La gráfica de resultados requiere ajuste de ejes, pues el índice de filas es el eje x, y las columnas es el eje y. La matriz y sus datos en la gráfica se obtiene como la transpuesta de u

# GRAFICA 3D ------
import matplotlib.pyplot as plt
# Malla para cada eje X,Y
Xi, Yi = np.meshgrid(xi, yj)
U = np.transpose(u) # ajuste de índices fila es x

fig_3D = plt.figure()
graf_3D = fig_3D.add_subplot(111, projection='3d')
graf_3D.plot_wireframe(Xi,Yi,U,
                       color ='blue')
graf_3D.plot(Xi[1,0],Yi[1,0],U[1,0],'o',color ='orange')
graf_3D.plot(Xi[1,1],Yi[1,1],U[1,1],'o',color ='red',
             label='U[i,j]')
graf_3D.plot(Xi[1,2],Yi[1,2],U[1,2],'o',color ='salmon')
graf_3D.plot(Xi[0,1],Yi[0,1],U[0,1],'o',color ='green')
graf_3D.plot(Xi[2,1],Yi[2,1],U[2,1],'o',color ='salmon')

graf_3D.set_title('EDP elíptica')
graf_3D.set_xlabel('x')
graf_3D.set_ylabel('y')
graf_3D.set_zlabel('U')
graf_3D.legend()
graf_3D.view_init(35, -45)
plt.show()

EDP Elípticas [ concepto ] [ ejercicio ]
Método iterativo: [ Analítico ] [ Algoritmo ]

 

7.2 EDP Elípticas

EDP Elípticas [ concepto ] Método [ iterativo ] [ implícito ]
..


1. EDP Elípticas

Referencia: Chapra 29.1 p866, Rodríguez 10.3 p425, Burden 12.1 p694

Las Ecuaciones Diferenciales Parciales tipo elípticas semejantes a la mostrada:

\frac{\partial ^2 u}{\partial x^2} + \frac{\partial ^2 u}{ \partial y^2} = 0

(ecuación de Laplace, Ecuación de Poisson con f(x,y)=0)

se interpreta como una placa metálica de dimensiones Lx y Ly, delgada con aislante que recubren las caras de la placa, y sometidas a condiciones en las fronteras:

Lx = dimensión x de placa metálica
Ly = dimensión y de placa metálica
u[0,y]  = Ta
u[Lx,y] = Tb
u[x,0]  = Tc
u[x,Ly] = Td

Placa Metalica 01

Para el planteamiento se usa una malla en la que cada nodo corresponden a los valores u[xi,yj]. Para simplificar la nomenclatura se usan los subíndices i para el eje de las x,  j para el eje t, quedando u[i,j].

EDP Elipticas Iterativo

vista superior, plano xy:

Placa Metalica 02

La ecuación se cambia a la forma discreta, usando diferencias divididas que se sustituyen en la ecuación, por ejemplo:

\frac{u_{i+1,j}-2u_{i,j}+u_{i-1,j}}{(\Delta x)^2} + \frac{u_{i,j+1}-2u_{i,j}+u_{i,j-1}}{(\Delta y)^2}=0

Se agrupan los términos de los diferenciales:

\frac{(\Delta y)^2}{(\Delta x)^2} \Big( u_{i+1,j}-2u_{i,j} +u_{i-1,j} \Big)+ u_{i,j+1}-2u_{i,j}+u_{i,j-1}=0 \lambda \Big( u_{i+1,j}-2u_{i,j} +u_{i-1,j} \Big)+ u_{i,j+1}-2u_{i,j}+u_{i,j-1}=0 \lambda u_{i+1,j} - 2\lambda u_{i,j} + \lambda u_{i-1,j} + u_{i,j+1} - 2 u_{i,j} + u_{i,j-1} = 0 \lambda u_{i+1,j} + (-2\lambda-2) u_{i,j} +\lambda u_{i-1,j} + u_{i,j+1} +u_{i,j-1}=0

Obteniendo así la solución numérica conceptual. La forma de resolver el problema determina el nombre del método a seguir.

Una forma de simplificar la expresión, es hacer λ=1, es decir  \lambda = \frac{(\Delta y)^2}{(\Delta x)^2} = 1, se determina que los tamaño de paso Δx, Δy son iguales.

u_{i+1,j}-4u_{i,j}+u_{i-1,j} + u_{i,j+1} +u_{i,j-1} = 0

EDP Elípticas [ concepto ] Método [ iterativo ] [ implícito ]