8.3.1 Regresión polinomial de grado m – Ejercicio Temperatura para un día

De una estación meteorológica se obtiene un archivo.csv con los datos de los sensores disponibles durante una semana.

2021OctubreEstMetorologica.csv


Para analizar el comportamiento de la variable de temperatura, se requiere disponer de un modelo polinomial que describa la temperatura a lo largo del día, para un solo día.

Como valores de variable independiente utilice un equivalente numérico decimal de la hora, minuto y segundo.

a. Realice la lectura del archivo, puede usar las instrucciones descritas en el  enlace: Archivos.csv con Python – Ejercicio con gráfica de temperatura y Humedad (CCPG1001: Fundamentos de programación)

b. Seleccione los datos del día 1 del mes para realizar la gráfica, y convierta tiempo en equivalente decimal.

c. Mediante pruebas, determine el grado del polinomio más apropiado para minimizar los errores.

Desarrollo

literal a

Se empieza con las instrucciones del enlace dadas añadiendo los parpametros de entrada como undia = 0 y grado del polinomio como gradom = 2 como el ejercicio anterior.

literal b

En el modelo polinomial, el eje x es un número decimal, por lo que los valores de hora:minuto:segundo se convierte a un valor decimal. Para el valor decimal de xi se usa la unidad de horas y las fracciones proporcionales de cada minuto y segundo.

literal c

Se inicia con el valor de gradom = 2, observando el resultado se puede ir incrementando el grado del polinomio observando los parámetros de error y coeficiente de determinación hasta cumplir con la tolerancia esperada segun la aplicación del ejercicio.

Resultados obtenidos:

columnas en tabla: 
Index(['No', 'Date', 'Time', 'ColdJunc0', 'PowerVolt', 'PowerKind', 'WS(ave)',
       'WD(ave)', 'WS(max)', 'WD(most)', 'WS(inst_m)', 'WD(inst_m)',
       'Max_time', 'Solar_rad', 'TEMP', 'Humidity', 'Rainfall', 'Bar_press.',
       'Soil_temp', 'fecha', 'horadec'],
      dtype='object')
ymedia =  25.036805555555553
 f = -6.57404678141012e-6*x**6 + 0.00052869201494877*x**5 - 0.0152875582352169*x**4 + 0.184200388015364*x**3 - 0.761164009032398*x**2 + 0.393278389794015*x + 22.1142936414255
coef_determinacion r2 =  0.9860621424061621
98.61% de los datos
     se describe con el modelo

Instrucciones en Python

# lecturas archivo.csv de estación meteorológica
import numpy as np
import sympy as sym
import pandas as pd
import matplotlib.pyplot as plt
from matplotlib.dates import DateFormatter, DayLocator

# INGRESO
narchivo = "2021OctubreEstMetorologica.csv"
sensor = 'TEMP'
undia  = 2 # dia a revisar
gradom = 6 # grado del polinomio

# PROCEDIMIENTO
tabla = pd.read_csv(narchivo, sep=';',decimal=',')
n_tabla = len(tabla)

# fechas concatenando columnas de texto
tabla['fecha'] = tabla['Date']+' '+tabla['Time']

# convierte a datetime
fechaformato = "%d/%m/%Y %H:%M:%S"
tabla['fecha'] = pd.to_datetime(tabla['fecha'],
                                format=fechaformato)
# serie por dia, busca indices
diaIndice = [0] # indice inicial
for i in range(1,n_tabla-1,1):
    i0_fecha = tabla['fecha'][i-1]
    i1_fecha = tabla['fecha'][i]
    if i0_fecha.day != i1_fecha.day:
        diaIndice.append(i)
diaIndice.append(len(tabla)-1) # indice final
m = len(diaIndice)

# horas decimales en un día
horadia = tabla['fecha'].dt.hour
horamin = tabla['fecha'].dt.minute
horaseg = tabla['fecha'].dt.second
tabla['horadec'] = horadia + horamin/60 + horaseg/3600

# PROCEDIMIENTO Regresión Polinomial grado m
m = gradom
# Selecciona dia
i0 = diaIndice[undia]
i1 = diaIndice[undia+1]
# valores a usar en regresión
xi = tabla['horadec'][i0:i1]
yi = tabla[sensor][i0:i1]
n  = len(xi)

# llenar matriz a y vector B
k = m + 1
A = np.zeros(shape=(k,k),dtype=float)
B = np.zeros(k,dtype=float)
for i in range(0,k,1):
    for j in range(0,i+1,1):
        coeficiente = np.sum(xi**(i+j))
        A[i,j] = coeficiente
        A[j,i] = coeficiente
    B[i] = np.sum(yi*(xi**i))

# coeficientes, resuelve sistema de ecuaciones
C = np.linalg.solve(A,B)

# polinomio
x = sym.Symbol('x')
f = 0
fetiq = 0
for i in range(0,k,1):
    f = f + C[i]*(x**i)
    fetiq = fetiq + np.around(C[i],4)*(x**i)

fx = sym.lambdify(x,f)
fi = fx(xi)

# errores
ym = np.mean(yi)
xm = np.mean(xi)
errado = fi - yi

sr = np.sum((yi-fi)**2)
syx = np.sqrt(sr/(n-(m+1)))
st = np.sum((yi-ym)**2)

# coeficiente de determinacion
r2 = (st-sr)/st
r2_porcentaje = np.around(r2*100,2)

# SALIDA
print(' columnas en tabla: ')
print(tabla.keys())
print('ymedia = ',ym)
print(' f =',f)
print('coef_determinacion r2 = ',r2)
print(str(r2_porcentaje)+'% de los datos')
print('     se describe con el modelo')

# grafica
plt.plot(xi,yi,'o',label='(xi,yi)')
plt.plot(xi,fi, color='orange', label=fetiq)

# lineas de error
for i in range(i0,i1,1):
    y0 = np.min([yi[i],fi[i]])
    y1 = np.max([yi[i],fi[i]])
    plt.vlines(xi[i],y0,y1, color='red',
               linestyle = 'dotted')

plt.xlabel('xi - hora en decimal')
plt.ylabel('yi - '+ sensor)
plt.legend()
etiq_titulo = sensor+ ' dia '+str(undia)
plt.title(etiq_titulo+': Regresion polinomial grado '+str(m))
plt.show()

Tarea

Determinar el polinomio de regresión para los días 3 y 5, y repetir el proceso para el sensor de Humedad (‘Humidity’)

Referencia: Archivos.csv con Python – Ejercicio con gráfica de temperatura y Humedad en Fundamentos de Programación

8.3 Regresión polinomial de grado m

Referencia: Chapra 17.2 p482, Burden 8.1 p501

Una alternativa a usar transformaciones ente los ejes para ajustar las curvas es usar regresión polinomial extendiendo el procedimiento de los mínimos cuadrados.

Suponga que la curva se ajusta a un polinomio de segundo grado o cuadrático

y = a_0 + a_1 x + a_2 x^2 +e

usando nuevamente la suma de los cuadrados de los residuos, se tiene,

S_r = \sum_{i=1}^n (y_i- (a_0 + a_1 x_i + a_2 x_i^2))^2

para minimizar los errores se deriva respecto a cada coeficiente: a0, a1, a2

\frac{\partial S_r}{\partial a_0} = -2\sum (y_i- a_0 - a_1 x_i - a_2 x_i^2) \frac{\partial S_r}{\partial a_1} = -2\sum x_i (y_i- a_0 - a_1 x_i - a_2 x_i^2) \frac{\partial S_r}{\partial a_2} = -2\sum x_i^2 (y_i- a_0 - a_1 x_i - a_2 x_i^2)

Se busca el mínimo de cada sumatoria, igualando a cero la derivada y reordenando, se tiene el conjunto de ecuaciones:

a_0 (n) + a_1 \sum x_i + a_2 \sum x_i^2 = \sum y_i a_0 \sum x_i + a_1 \sum x_i^2 + a_2 \sum x_i^3 =\sum x_i y_i a_0 \sum x_i^2 + a_1 \sum x_i^3 + a_2 \sum x_i^4 =\sum x_i^2 y_i

con incógnitas a0, a1 y a2, cuyos coeficientes se pueden evaluar con los puntos observados. Se puede usar un médoto directo de la unidad 3 para resolver el sistema de ecuaciones Ax=B.

\begin{bmatrix} n & \sum x_i & \sum x_i^2 \\ \sum x_i & \sum x_i^2 & \sum x_i^3 \\ \sum x_i^2 & \sum x_i^3 & \sum x_i^4 \end{bmatrix} \begin{bmatrix} a_0 \\ a_1 \\ a_2 \end{bmatrix} = \begin{bmatrix} \sum y_i \\ \sum x_i y_i \\ \sum x_i^2 y_i \end{bmatrix}
C = np.linalg.solve(A,B)
y = C_0 + C_1 x + C_2 x^2

El error estándar se obtiene como:

S_{y/x} = \sqrt{\frac{S_r}{n-(m+1)}}

siendo m el grado del polinomio usado, para el caso presentado m = 2

S_r = \sum{(yi-fi)^2}

Sr es la suma de los cuadrados de los residuos alrededor de la línea de regresión.

xi yi (yi-ymedia)2 (yi-fi)2
∑yi St Sr
r^2 = \frac{S_t-S_r}{S_t} S_t = \sum{(yi-ym)^2}

siendo St la suma total de los cuadrados alrededor de la media para la variable dependiente y. Este valor es la magnitud del error residual asociado con la variable dependiente antes de la regresión.

El algoritmo se puede desarrollar de forma semejante a la presentada en la sección anterior,

Ejercicio Chapra 17.5 p484

Los datos de ejemplo para la referencia son:

xi = [0,   1,    2,    3,    4,   5]
yi = [2.1, 7.7, 13.6, 27.2, 40.9, 61.1]

resultado de algoritmo:

fx  =  1.86071428571429*x**2 + 2.35928571428571*x + 2.47857142857143
Syx =  1.117522770621316
coef_determinacion r2 =  0.9985093572984048
99.85% de los datos se describe con el modelo
>>> 

Instrucciones en Python

# regresion con polinomio grado m=2
import numpy as np
import sympy as sym
import matplotlib.pyplot as plt

# INGRESO
xi = [0,   1,    2,    3,    4,   5]
yi = [2.1, 7.7, 13.6, 27.2, 40.9, 61.1]
m  = 2

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

# llenar matriz a y vector B
k = m + 1
A = np.zeros(shape=(k,k),dtype=float)
B = np.zeros(k,dtype=float)
for i in range(0,k,1):
    for j in range(0,i+1,1):
        coeficiente = np.sum(xi**(i+j))
        A[i,j] = coeficiente
        A[j,i] = coeficiente
    B[i] = np.sum(yi*(xi**i))

# coeficientes, resuelve sistema de ecuaciones
C = np.linalg.solve(A,B)

# polinomio
x = sym.Symbol('x')
f = 0
fetiq = 0
for i in range(0,k,1):
    f = f + C[i]*(x**i)
    fetiq = fetiq + np.around(C[i],4)*(x**i)

fx = sym.lambdify(x,f)
fi = fx(xi)

# errores
ym = np.mean(yi)
xm = np.mean(xi)
errado = fi - yi

sr = np.sum((yi-fi)**2)
syx = np.sqrt(sr/(n-(m+1)))
st = np.sum((yi-ym)**2)

# coeficiente de determinacion
r2 = (st-sr)/st
r2_porcentaje = np.around(r2*100,2)

# SALIDA
print('ymedia = ',ym)
print(' f =',f)
print('coef_determinacion r2 = ',r2)
print(str(r2_porcentaje)+'% de los datos')
print('     se describe con el modelo')

# grafica
plt.plot(xi,yi,'o',label='(xi,yi)')
# plt.stem(xi,yi, bottom=ym)
plt.plot(xi,fi, color='orange', label=fetiq)

# lineas de error
for i in range(0,n,1):
    y0 = np.min([yi[i],fi[i]])
    y1 = np.max([yi[i],fi[i]])
    plt.vlines(xi[i],y0,y1, color='red',
               linestyle = 'dotted')

plt.xlabel('xi')
plt.ylabel('yi')
plt.legend()
plt.title('Regresion polinomial grado '+str(m))
plt.show()

8.2 Regresión por Mínimos cuadrados con Python

Referencia: Chapra 17.1.2 p 469. Burden 8.1 p499

Criterio de ajuste a una línea recta por mínimos cuadrados

Una aproximación de la relación entre los puntos xi, yi por medio de un polinomio de grado 1, busca minimizar la suma de los errores residuales de los datos.

y_{i,modelo} = p_1(x) = a_0 + a_1 x_i

Se busca el valor mínimo para los cuadrados de las diferencias entre los valores de yi con la línea recta.

S_r = \sum_{i=1}^{n} \Big( y_{i,medida} - y_{i,modelo} \Big)^2 S_r= \sum_{i=1}^{n} \Big( y_{i} - (a_0 + a_1 x_i) \Big)^2

para que el error acumulado sea mínimo, se deriva respecto a los coeficientes de la recta a0 y a1 y se iguala a cero,

\frac{\partial S_r}{\partial a_0}= (-1)2 \sum_{i=1}^{n} \Big( y_{i} - a_0 - a_1 x_i \Big) \frac{\partial S_r}{\partial a_1}= (-1)2 \sum_{i=1}^{n} \Big( y_{i} - a_0 - a_1 x_i \Big)x_i 0 = \sum_{i=1}^{n} y_i - \sum_{i=1}^{n} a_0 - \sum_{i=1}^{n} a_1 x_i 0= \sum_{i=1}^{n} y_i x_i - \sum_{i=1}^{n} a_0x_i - \sum_{i=1}^{n} a_1 x_i^2

simplificando,

\sum_{i=1}^{n} a_0 + a_1 \sum_{i=1}^{n} x_i = \sum_{i=1}^{n} y_{i} a_0 \sum_{i=1}^{n} x_i + a_1 \sum_{i=1}^{n} x_i^2 = \sum_{i=1}^{n} y_i x_i

que es un conjunto de dos ecuaciones lineales simultaneas con dos incógnitas a0 y a1, los coeficientes del sistema de ecuaciones son las sumatorias que se obtienen completando la siguiente tabla:

Tabla de datos y columnas para ∑ en la fórmula
xi yi xiyi xi2 yi2
x0 y0 x0y0 x02 y02
∑ xi ∑ yi ∑ xiyi ∑ xi2 ∑ yi2
\begin{bmatrix} n & \sum x_i \\ \sum x_i & \sum x_i^2 \end{bmatrix} \begin{bmatrix} a_0 \\ a_1 \end{bmatrix} = \begin{bmatrix} \sum y_i \\ \sum x_i y_i \end{bmatrix}

cuya solución es:

a_1 = \frac{n \sum x_i y_i - \sum x_i \sum y_i}{n \sum x_i^2 - \big( \sum x_i \big) ^2 } a_0 = \overline{y} - a_1 \overline{x}

usando la media de los valores en cada eje para encontrar a0

Coeficiente de correlación

El coeficiente de correlación se puede obtener con las sumatorias anteriores usando la siguiente expresión:

r= \frac{n \sum x_i y_i - \big( \sum x_i \big) \big( \sum y_i\big)} {\sqrt{n \sum x_i^2 -\big(\sum x_i \big) ^2 }\sqrt{n \sum y_i^2 - \big( \sum y_i \big)^2}}

En un ajuste perfecto, Sr = 0 y r = r2 = 1, significa que la línea explica
el 100% de la variabilidad de los datos.

Aunque el coeficiente de correlación ofrece una manera fácil de medir la bondad del ajuste, se deberá tener cuidado de no darle más significado del que ya tiene.

El solo hecho de que r sea “cercana” a 1 no necesariamente significa que el ajuste sea “bueno”. Por ejemplo, es posible obtener un valor relativamente alto de r cuando la relación entre y y x no es lineal.

Los resultados indicarán que el modelo lineal explicó r2 % de la incertidumbre original


Algoritmo en Python

Siguiendo con los datos propuestos del ejemplo en Chapra 17.1 p470:

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

Aplicando las ecuaciones para a0 y a1 se tiene el siguiente resultado para los datos de prueba:

 f =  0.839285714285714*x + 0.0714285714285712
coef_correlación   r  =  0.9318356132188194
coef_determinación r2 =  0.8683176100628931
86.83% de los datos está descrito en el modelo lineal
>>>

con las instrucciones:

# mínimos cuadrados, regresión con polinomio grado 1
import numpy as np
import sympy as sym
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)

# sumatorias y medias
xm  = np.mean(xi)
ym  = np.mean(yi)
sx  = np.sum(xi)
sy  = np.sum(yi)
sxy = np.sum(xi*yi)
sx2 = np.sum(xi**2)
sy2 = np.sum(yi**2)

# coeficientes a0 y a1
a1 = (n*sxy-sx*sy)/(n*sx2-sx**2)
a0 = ym - a1*xm

# polinomio grado 1
x = sym.Symbol('x')
f = a0 + a1*x

fx = sym.lambdify(x,f)
fi = fx(xi)

# coeficiente de correlación
numerador = n*sxy - sx*sy
raiz1 = np.sqrt(n*sx2-sx**2)
raiz2 = np.sqrt(n*sy2-sy**2)
r = numerador/(raiz1*raiz2)

# coeficiente de determinacion
r2 = r**2
r2_porcentaje = np.around(r2*100,2)

# SALIDA
# print('ymedia =',ym)
print(' f = ',f)
print('coef_correlación   r  = ', r)
print('coef_determinación r2 = ', r2)
print(str(r2_porcentaje)+'% de los datos')
print('     está descrito en el modelo lineal')

# grafica
plt.plot(xi,yi,'o',label='(xi,yi)')
# plt.stem(xi,yi,bottom=ym,linefmt ='--')
plt.plot(xi,fi, color='orange',  label=f)

# lineas de error
for i in range(0,n,1):
    y0 = np.min([yi[i],fi[i]])
    y1 = np.max([yi[i],fi[i]])
    plt.vlines(xi[i],y0,y1, color='red',
               linestyle ='dotted')
plt.legend()
plt.xlabel('xi')
plt.title('minimos cuadrados')
plt.show()

Coeficiente de correlación con Numpy

Tambien es posible usar la librería numpy para obtener el resultado anterior,

>>> coeficientes = np.corrcoef(xi,yi)
>>> coeficientes
array([[1.        , 0.93183561],
       [0.93183561, 1.        ]])
>>> r = coeficientes[0,1]
>>> r
0.9318356132188195

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.

En el ejemplo de Chapra 17.1 p470, se presentan los datos de un experimento mostados 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:

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

2.6 Sistemas de Ecuaciones no lineales – Newton-Raphson

Referencia: Chapra 6.5 p162, Chapra Ejercicio 6.11 p.166/pdf190

Con el método de Newton-Raphson para múltiples ecuaciones, determine las raíces para:

x^2+xy =10 y + 3xy^2 = 57

Observe que un par correcto de raíces es x=2 y y=3.
Use como valore iniciales x=1.5, y=3.5

Planteamiento

Las ecuaciones se expresan de la forma f(x,y) = 0

x^2+xy -10 = 0 y + 3xy^2 -57 = 0

Se puede usar extensines de los métodos abiertos para resolver ecuacioens simples, por ejemplo Newton-Raphson.

u_{i+1} = u_i + (x_{i+1}-x_i)\frac{\partial u_i}{\partial x} + (y_{i+1}-y_i) \frac{\partial u_i}{\partial y} v_{i+1} = v_i + (x_{i+1}-x_i)\frac{\partial v_i}{\partial x} + (y_{i+1}-y_i) \frac{\partial v_i}{\partial y}

ecuaciones que se pueden reordenar y encontrar la solución a partir de la matriz Jacobiano.

Instrucciones en Python

Usando un algoritmo para resolver el Jacobiano y estimar los puntos luego de cada iteración se obtienen:

iteración:  1
Jacobiano con puntos iniciales: 
Matrix([[6.50000000000000, 1.50000000000000], [36.7500000000000, 32.5000000000000]])
determinante:  156.12499999999994
puntos xi,yi: 2.03602882305845 2.84387510008006
error: 0.656124899919936
iteración:  2
Jacobiano con puntos iniciales: 
Matrix([[6.91593274619696, 2.03602882305845], [24.2628767545662, 35.7412700376474]])
determinante:  197.78430344142245
puntos xi,yi: 1.99870060905582 3.00228856292451
error: 0.158413462844444
iteración:  3
Jacobiano con puntos iniciales: 
Matrix([[6.99968978103616, 1.99870060905582], [27.0412098452019, 37.0040558756713]])
determinante:  204.96962918261596
puntos xi,yi: 1.99999998387626 2.99999941338891
error: 0.00228914953559523
iteración:  4
Jacobiano con puntos iniciales: 
Matrix([[6.99999938114143, 1.99999998387626], [26.9999894410015, 36.9999926704397]])
determinante:  204.9999473486533
puntos xi,yi: 1.99999999999998 3.00000000000008
error: 5.86611161867978e-7
Resultado: 
1.99999999999998 3.00000000000008
>>> 

Algoritmo presentado para dos ecuaciones y dos incógnitas, en la unidad 3 se puede ampliar la propuesta. Revisar el método de Gauss-Seidel y Jacobi.

# Ejercicio Chapra Ej:6.11
# Sistemas de ecuaciones no lineales
# con método de Newton Raphson para xy

import numpy as np
import sympy as sym

def matrizJacobiano(variables, funciones):
    n = len(funciones)
    m = len(variables)
    # matriz Jacobiano inicia con ceros
    Jcb = sym.zeros(n,m)
    for i in range(0,n,1):
        unafi = sym.sympify(funciones[i])
        for j in range(0,m,1):
            unavariable = variables[j]
            Jcb[i,j] = sym.diff(unafi, unavariable)
    return Jcb

# PROGRAMA ----------
# INGRESO
x = sym.Symbol('x')
y = sym.Symbol('y')

f1 = x**2 + x*y - 10
f2 = y + 3*x*(y**2)-57

x0 = 1.5
y0 = 3.5

tolera = 0.0001

# PROCEDIMIENTO
funciones = [f1,f2]
variables = [x,y]
n = len(funciones)
m = len(variables)

Jxy = matrizJacobiano(variables, funciones)

# valores iniciales
xi = x0
yi = y0

# tramo inicial, mayor que tolerancia
itera = 0
tramo = tolera*2

while (tramo>tolera):
    J = Jxy.subs([(x,xi),(y,yi)])

    # determinante de J
    Jn = np.array(J,dtype=float)
    determinante =  np.linalg.det(Jn)

    # iteraciones
    f1i = f1.subs([(x,xi),(y,yi)])
    f2i = f2.subs([(x,xi),(y,yi)])

    numerador1 = f1i*Jn[n-1,m-1]-f2i*Jn[0,m-1]
    xi1 = xi - numerador1/determinante
    numerador2 = f2i*Jn[0,0]-f1i*Jn[n-1,0]
    yi1 = yi -numerador2/determinante
    
    tramo = np.max(np.abs([xi1-xi,yi1-yi]))
    xi = xi1
    yi = yi1

    itera = itera +1
    print('iteración: ',itera)
    print('Jacobiano con puntos iniciales: ')
    print(J)
    print('determinante: ', determinante)
    print('puntos xi,yi:',xi,yi)
    print('error:',tramo)
    
# SALIDA
print('Resultado: ')
print(xi,yi)

Unidades – Temas de estudio

Interpolar – Mascota 50 años

Ejercicio de interpolación

  • Desarrollar por grupos
  • integrar soluciones en un mismo gráfico
  • mostrar polinomios y rangos de existencia

# caparazon superior
xi = [113., 117, 134, 153, 169, 184, 194, 203]
fi = [127., 141, 161, 166, 160, 155, 140, 132]
puntos = [[xi,fi]]

# caparazon inferior 1
xi = [113., 123, 130, 149, 182, 197, 208, 211]
fi = [127., 116, 112, 107, 110, 114, 112, 108]
puntos.append([xi,fi])

# caparazon inferior 2
xi = [107., 114, 120, 143, 170, 192, 210]
fi = [124., 114, 108, 99, 99, 106, 107]
puntos.append([xi,fi])

# Patas 01
xi = [110., 116, 120, 121, 134, 143, 153, 168, 173, 177, 181, 188, 194, 201, 205, 210, 214, 218]
fi = [92., 86, 89, 86, 87, 95, 91, 93, 85, 85, 83, 86, 88, 86, 91, 87, 92, 92]
puntos.append([xi,fi])

# Patas 02
xi = [109., 115, 117, 120, 125, 130, 134, 138, 142, 145, 150, 152, 154]
fi = [91., 87, 82, 82, 78, 78, 81, 78, 76, 79, 77, 80, 86]
puntos.append([xi,fi])

# Patas 03
xi = [172., 175, 182, 186, 191, 195, 201, 205, 210, 211, 217]
fi = [86., 79, 79, 77, 82, 78, 81, 79, 82, 85, 90]
puntos.append([xi,fi])

# Patas 04
xi = [113., 118, 122, 127, 132, 136, 140, 144, 149, 152]
fi = [88., 90, 88, 86, 85, 83, 87, 84, 86, 85]
puntos.append([xi,fi])

# Rabito 01
xi = [97., 102, 108, 111, 117, 121]
fi = [120., 113, 108, 105, 102, 102]
puntos.append([xi,fi])

# cabeza 01
xi = [194., 196, 203, 210, 211, 209]
fi = [177., 167, 157, 149, 138, 135]
puntos.append([xi,fi])

# cabeza 02
xi = [195., 199, 207, 214, 219, 224, 229, 234, 239, 242, 244]
fi = [177., 185, 190, 190, 188, 193, 192, 185, 175, 172, 168]
puntos.append([xi,fi])

# cabeza 03
xi = [220., 226, 234, 239, 241]
fi = [164., 162, 163, 163, 163]
puntos.append([xi,fi])

# cabeza 04
xi = [203., 211, 214, 219, 223, 224, 225, 230, 236, 241]
fi = [115., 119, 125, 124, 127, 137, 146, 149, 154, 162]
puntos.append([xi,fi])

# cabeza 05
xi = [208., 212, 215, 219, 221, 225, 228]
fi = [174., 177, 178, 179, 182, 183, 178]
puntos.append([xi,fi])

# cabeza 06
xi = [206., 210, 214, 218, 220, 224, 229, 233, 232]
fi = [179., 182, 182, 181, 174, 171, 170, 175, 180]
puntos.append([xi,fi])

7.3 EDP hiperbólicas

Referencia:  Chapra PT8.1 p.860 pdf.884,  Rodriguez 10.4 p.435

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}

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

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}

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

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:

Algoritmo en 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
xn = 1 # Longitud de cuerda
y0 = 0
yn = 0 # Puntos de amarre
t0 = 0
c = 2
# discretiza
tramosx = 16
tramost = 32
dx = (xn-x0)/tramosx 
dt = dx/c

# PROCEDIMIENTO
xi = np.arange(x0,xn+dx,dx)
tj = np.arange(0,tramost*dt+dt,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.

7.2.2 EDP Elípticas método implícito

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

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.  ]]
>>>

Algoritmo usado para resolver el problema:

# 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,dx)
yj = np.arange(y0,yn+dy,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.

# Gráfica
import matplotlib.pyplot as plt
from matplotlib import cm
from mpl_toolkits.mplot3d import Axes3D

X, Y = np.meshgrid(xi, yj)
U = np.transpose(u) # ajuste de índices fila es x

figura = plt.figure()
ax = Axes3D(figura)
ax.plot_surface(X, Y, U, rstride=1, cstride=1, cmap=cm.Reds)

plt.title('EDP elíptica')
plt.xlabel('x')
plt.ylabel('y')
plt.show()

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.