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

EDP Elípticas [ contínua a discreta ] || [ sympy_iterativo ] [ sympy_implícito ]
..


1. Ejercicio

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

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

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

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

EDP Elípticas [ contínua a discreta ] || [ sympy_iterativo ] [ sympy_implícito ]

..


2. EDP Elípticas contínua a discreta

El desarrollo analítico inicia convirtiendo la ecuación diferencial parcial elíptica de la forma contínua a su representación discreta usando las expresiones de derivadas en  diferencias divididas. El cambio de la expresión se realiza usando  Sympy.

La EDP se escribe en formato Sympy en el bloque de ingreso.  La ecuación EDP se compone del lado izquierdo (LHS) y lado derecho (RHS), indicando u como una función de las variables x,y.

# INGRESO
fxy = 0*x+0*y # f(x,y) = 0 , ecuacion de Poisson
# ecuacion: LHS=RHS
LHS = sym.diff(u,x,2) + sym.diff(u,y,2)
RHS = fxy
EDP = sym.Eq(LHS-RHS,0)

las diferencias divididas se ingresan como un diccionario:

dif_dividida ={sym.diff(u,x,2): (U.subs(i,i-1)-2*U+U.subs(i,i+1))/(Dx**2),
               sym.diff(u,y,2): (U.subs(j,j-1)-2*U+U.subs(j,j+1))/(Dy**2),
               sym.diff(u,y,1): (U - U.subs(j,j-1))/Dy}

las condiciones iniciales en los bordes, pueden ser funciones matemáticas, por lo que se usa el formato lambda. En el ejercicio básico presentado en la parte teórica, los bordes tienen temperaturas constantes, por lo que en ésta sección se escriben las condiciones como la constante mas cero por la variable independiente, facilitando la evaluación de un vector xi por ejemplo.

# Condiciones iniciales en los bordes
fya = lambda y: 60 +0*y  # izquierda
fyb = lambda y: 25 +0*y  # derecha
fxc = lambda x: 50 +0*x  # inferior, función inicial 
fxd = lambda x: 70 +0*x  # superior, función inicial 

los demás parámetros del ejercicio se ingresan más adelante de forma semejante al ejercicio con Numpy.

Todas las expresiones se escriben al lado izquierdo (LHS) de la igualdad, la parte del lado derecho(RHS) se la prefiere mantener en cero. La expresión se organiza manteniendo el coeficiente en 1 para el termino de Dx de mayor orden. Se sustituyen las derivadas por diferencias divididas, para luego simplificar la expresión usando λ como referencia.

El resultado de las operaciones en la expresión usando el algoritmo es:

EDP Elíptica contínua a discreta
 ordenDx : 2
 ordenDy : 2

 edp=0 :
  2              2             
 d              d              
---(u(x, y)) + ---(u(x, y)) = 0
  2              2             
dx             dy              
 K_ : 1.0

 discreta=0 :
-2*U(i, j) + U(i, j - 1) + U(i, j + 1)   -2*U(i, j) + U(i - 1, j) + U(i + 1, j)
-------------------------------------- + --------------------------------------
                   2                                        2                 
                 Dy                                       Dx                  
 Lambda_L :
  2
Dy 
---
  2
Dx 
 Lambda L_k : 1

 (discreta=0)*Dy**ordeny :
                                             2             2                 2
                                         2*Dy *U(i, j)   Dy *U(i - 1, j)   Dy *U(i + 1, j)
-2*U(i, j) + U(i, j - 1) + U(i, j + 1) - ------------- + --------------- + ---------------
                                                2                2                 2
                                              Dx               Dx                Dx  
 discreta_L=0 :
L*U(i - 1, j) + L*U(i + 1, j) + (-2*L - 2)*U(i, j) + U(i, j - 1) + U(i, j + 1) = 0

 discreta :
-4*U(i, j) + U(i, j - 1) + U(i, j + 1) + U(i - 1, j) + U(i + 1, j) = 0

 discreta_iterativa :
          U(i, j - 1) + U(i, j + 1) + U(i - 1, j) + U(i + 1, j)
U(i, j) = -----------------------------------------------------
                                    4                          

Instrucciones en Python

Los resultados al aplicar una operación a la expresión se guardan en un diccionario, para revisar cada paso intermedio al final

# Ecuaciones Diferenciales Parciales Elipticas
# EDP Elípticas contínua a discreta con Sympy
import numpy as np
import sympy as sym

# funciones continuas y variables simbólicas usadas
y = sym.Symbol('y',real=True)
x = sym.Symbol('x',real=True)
u = sym.Function('u')(x,y)
f = sym.Function('f')(x,y) # funcion complemento
# funciones discretas y variables simbólicas usadas
i  = sym.Symbol('i',integer=True,positive=True)
j  = sym.Symbol('j',integer=True,positive=True)
Dx = sym.Symbol('Dx',real=True,positive=True)
Dy = sym.Symbol('Dy',real=True,positive=True)
L  = sym.Symbol('L',real=True)
U  = sym.Function('U')(i,j)

# INGRESO
fxy = 0*x+0*y  # f(x,y) = 0 , ecuacion de Poisson
# ecuacion edp : LHS=RHS
LHS = sym.diff(u,x,2) + sym.diff(u,y,2)
RHS = fxy
edp = sym.Eq(LHS-RHS,0)

dif_dividida ={sym.diff(u,x,2): (U.subs(i,i-1)-2*U+U.subs(i,i+1))/(Dx**2),
               sym.diff(u,y,2): (U.subs(j,j-1)-2*U+U.subs(j,j+1))/(Dy**2),
               sym.diff(u,y,1): (U - U.subs(j,j-1))/Dy}

# Condiciones iniciales en los bordes
fya = lambda y: 60 +0*y  # izquierda
fyb = lambda y: 25 +0*y  # derecha
fxc = lambda x: 50 +0*x  # inferior 
fxd = lambda x: 70 +0*x  # 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
tolera = 0.0001

verdigitos = 2 # para mostrar en tabla de resultados

# PROCEDIMIENTO
def edp_discreta(edp,dif_dividida,x,y,u):
    ''' EDP contínua a discreta, usa diferencias divididas
        proporcionadas en parámetros, indica las variables x,y
        con función u de (x,y)
    '''
    resultado={}
    # expresión todo a la izquierda LHS (Left Hand side)
    LHS = edp.lhs
    RHS = edp.rhs
    if not(edp.rhs==0):
        LHS = LHS-RHS
        RHS = 0
        edp = sym.Eq(LHS,RHS)
    # para orden de derivada por x, y
    edp_x = edp.subs(x,0)
    edp_y = edp.subs(y,0)
    ordenDx = sym.ode_order(edp_x,u)
    ordenDy = sym.ode_order(edp_y,u)
    resultado['ordenDx'] = ordenDx
    resultado['ordenDy'] = ordenDy
    # convierte coeficiente d2u/dx2 a 1
    coeff_x = edp_coef_Dx(edp,x,ordenDx)
    if not(coeff_x==1):
        LHS = LHS/coeff_x
        RHS = RHS/coeff_x
    edp = sym.Eq(LHS,RHS)
    K_ = edp_coef_Dx(edp,y,ordenDy)
    resultado['edp=0']  = edp
    resultado['K_']  = K_
    # edp discreta
    discreta = edp.lhs 
    for derivada in dif_dividida:
        discreta = discreta.subs(derivada,dif_dividida[derivada])
    resultado['discreta=0'] = discreta
    return (resultado)

def edp_coef_Dx(edp,x,ordenx):
    ''' Extrae el coeficiente de la derivada Dx de ordenx
    '''
    coeff_x = 1.0
    # separa cada término de suma
    term_suma = sym.Add.make_args(edp.lhs)
    for term_k in term_suma:
        if term_k.is_Mul: # término de mas de un factor
            factor_Mul = sym.Mul.make_args(term_k)
            coeff_temp = 1; coeff_usar=0
            # separa cada factor de término 
            for factor_k in factor_Mul:
                if not(factor_k.is_Derivative):
                    coeff_temp = coeff_temp*factor_k
                else: # factor con derivada de ordenx
                    partes = factor_k.args
                    if partes[1]==(x,ordenx):
                        coeff_usar = 1
            if coeff_usar==1:
                coeff_x = coeff_x*coeff_temp
    return(coeff_x)

def edp_simplificaLamba(resultado,dx,dy):
    '''simplifica ecuacion con valores de lambda, dx y dy
    entregando la edp discreta simplificada
    '''
    discreta = resultado['discreta=0']
    ordenDy = resultado['ordenDy']
    ordenDx = resultado['ordenDx']
    K_ = resultado['K_']
    lamb = (Dy**ordenDy)/(Dx**ordenDx)
    if ordenDy==1 and ordenDx==2:
        lamb = lamb/K_
    resultado['Lambda_L'] = lamb
    # valor de Lambda en ecuacion edp
    L_k = lamb.subs([(Dx,dx),(Dy,dy),(1.0,1)])
    resultado['Lambda L_k'] = L_k
    # simplifica con lambda L
    discreta_L = sym.expand(discreta*(Dy**ordenDy),mul=True)
    resultado['(discreta=0)*Dy**ordeny'] = discreta_L
    discreta_L = edp_sustituye_L(resultado)
    discreta_L = discreta_L.subs(lamb,L)
    discreta_L = sym.collect(discreta_L,U)
    discreta_L = sym.Eq(discreta_L,0)
    resultado['discreta_L = 0'] = discreta_L
    # sustituye constantes en ecuación a iterar
    discreta_L = discreta_L.subs([(Dx,dx),(Dy,dy),
                                  (L,L_k),(1.0,1)])
    resultado['discreta'] = discreta_L
    return(resultado)

def edp_sustituye_L(resultado):
    ''' sustituye lambda con Dy**ordeny/Dx**x/K_
    por L, al simplificar Lambda
    '''
    discreta = resultado['(discreta=0)*Dy**ordeny']
    ordenDy = resultado['ordenDy']
    ordenDx = resultado['ordenDx']
    discreta_L = 0
    # separa cada término de suma
    term_suma = sym.Add.make_args(discreta)
    for term_k in term_suma:
        # busca partes de L y cambia por valor L
        cambiar = 0 # por orden de derivada
        if term_k.has(Dx) and term_k.has(Dy):
            partes = term_k.args
            ordeny=1
            ordenx=1
            for unaparte in partes:    
                if unaparte.has(Dy):
                    if unaparte.is_Pow:
                        partey = unaparte.args
                        ordeny = partey[1]
                if unaparte.has(Dx):
                    if unaparte.is_Pow:
                        partey = unaparte.args
                        ordenx = partey[1]
            if (ordeny<=ordenDy and ordenx<=-ordenDx):
                cambiar=1
        if cambiar:
            term_k = term_k*L/resultado['Lambda_L']
        discreta_L = discreta_L + term_k
        # simplifica unos con decimal a entero 
        discreta_L = discreta_L.subs(1.0,1)
    return(discreta_L)

# PROCEDIMIENTO
# transforma edp continua a discreta
resultado = edp_discreta(edp,dif_dividida,x,y,u)
resultado = edp_simplificaLamba(resultado,dx,dy)
discreta = resultado['discreta']

# Método iterativo para EDP Elíptica
# separa término U[j,j] del centro,con valor no conocido
buscar = U.subs(j,j)
discreta = sym.solve(discreta,buscar)
discreta = sym.factor_terms(discreta[0])
discreta = sym.Eq(buscar,discreta)
resultado['discreta_iterativa'] = discreta

# SALIDA
np.set_printoptions(precision=verdigitos)
algun_numero = [int,float,str,'Lambda L_k']
print('EDP Elíptica contínua a discreta')
for entrada in resultado:
    tipo = type(resultado[entrada])
    if tipo in algun_numero or entrada in algun_numero:
        print('',entrada,':',resultado[entrada])
    else:
        print('\n',entrada,':')
        sym.pprint(resultado[entrada])

EDP Elípticas [ contínua a discreta ] || [ sympy_iterativo ] [ sympy_implícito ]
..


3. EDP Elípticas – Método iterativo con Sympy-Python

El desarrollo del método iterativo para una ecuación diferencial parcial elíptica usando Sympy, sigue a continuación del anterior.

Se añaden las instrucciones para iterar la expresión discreta para cada i,j dentro de la matriz U creada para el ejercicio. Se aplican las condiciones iniciales a los bordes, para luego proceder a iterar.

Como en las actividades del curso se requiere realizar al menos 3 iteraciones para las expresiones del algoritmo con papel y lápiz, solo se presentarán los resultados para la primera fila en j=0

Los resultados del algoritmo  luego del resultado del algoritmo anterior se presentan como:

 discreta_iterativa :
          U(i, j - 1) + U(i, j + 1) + U(i - 1, j) + U(i + 1, j)
U(i, j) = -----------------------------------------------------
                                    4                          

 iterar i,j:
j:  1  ;  i:  1
          U(0, 1)   U(1, 0)   U(1, 2)   U(2, 1)
U(1, 1) = ------- + ------- + ------- + -------
             4         4         4         4   
U(i, j) = 53.125
j:  1  ;  i:  2
          U(1, 1)   U(2, 0)   U(2, 2)   U(3, 1)
U(2, 1) = ------- + ------- + ------- + -------
             4         4         4         4   
U(i, j) = 51.40625
j:  1  ;  i:  3
          U(2, 1)   U(3, 0)   U(3, 2)   U(4, 1)
U(3, 1) = ------- + ------- + ------- + -------
             4         4         4         4   
U(i, j) = 50.9765625
j:  1  ;  i:  4
          U(3, 1)   U(4, 0)   U(4, 2)   U(5, 1)
U(4, 1) = ------- + ------- + ------- + -------
             4         4         4         4   
U(i, j) = 50.869140625
j:  1  ;  i:  5
          U(4, 1)   U(5, 0)   U(5, 2)   U(6, 1)
U(5, 1) = ------- + ------- + ------- + -------
             4         4         4         4   
U(i, j) = 50.84228515625
j:  1  ;  i:  6
          U(5, 1)   U(6, 0)   U(6, 2)   U(7, 1)
U(6, 1) = ------- + ------- + ------- + -------
             4         4         4         4   
U(i, j) = 50.8355712890625
j:  1  ;  i:  7
          U(6, 1)   U(7, 0)   U(7, 2)   U(8, 1)
U(7, 1) = ------- + ------- + ------- + -------
             4         4         4         4   
U(i, j) = 44.271392822265625

Método iterativo EDP Elíptica
continuar el desarrollo con: 
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 ]
 j, U[i,j]
2 [60.   51.25 51.25 51.25 51.25 51.25 51.25 51.25 25.  ]
1 [60.   53.12 51.41 50.98 50.87 50.84 50.84 44.27 25.  ]
0 [50. 50. 50. 50. 50. 50. 50. 50. 50.]
>>> 

Se pueden sustituir los valores de las expresiones, evaluando cada u[i,j] y completando la matriz y para la gráfica. Esta parte queda como tarea.

Instrucciones en Python

Las instrucciones adicionales al algoritmo de EDP elíptica contínua a discreta empiezan con la creación de la matriz u_xy para los valores, los vectores xi,yj y una matriz u_mask que indica si se dispone de  un valor calculado o conocido en ese punto (x,y) . Para el método iterativo, la matriz se rellena con el promedio de los valores máximos de las funciones dadas para los bordes de la placa.

Tarea: El algoritmo desarrolla el cálculo para j_k=1, solo la primera fila. Se deja como tarea desarrollar para toda la matriz.

# ITERAR para cada i,j dentro de U ------------
# x[i] , y[j]  valor en posición en cada eje
xi = np.arange(x0,xn+dx,dx)
yj = np.arange(y0,yn+dy,dy)
n = len(xi)
m = len(yj)

# Matriz U
u_xy = np.zeros(shape=(n,m),dtype = float)
u_xy = u_xy*np.nan # valor inicial dentro de u
# llena u con valores en fronteras
u_xy[0,:]   = fya(yj)  # izquierda Ta
u_xy[n-1,:] = fyb(yj)  # derecha   Tb
u_xy[:,0]   = fxc(xi)  # inferior  Tc
u_xy[:,m-1] = fxd(xi)  # superior  Td
u0 = np.copy(u_xy)     # matriz u inicial

# u_mask[i,j] con valores iniciales o calculados:True
u_mask = np.zeros(shape=(n,m),dtype=bool)
u_mask[0,:]  = True # izquierda
u_mask[-1,:] = True # derecha
u_mask[:,0]  = True # inferior
u_mask[:,-1] = True # superior

# valor inicial de iteración dentro de u
# promedio = (Ta+Tb+Tc+Td)/4
promedio = (np.max(fya(yj))+np.max(fyb(yj))+\
            np.max(fxc(xi))+np.max(fxd(xi)))/4
u_xy[1:n-1,1:m-1]   = promedio
u_mask[1:n-1,1:m-1] = True
u0 = np.copy(u_xy) # copia para revisión

def edp_sustituyeValorU(discreta,xi,yj,u_xy,u_mask):
    '''Sustituye en edp discreta los valores conocidos de U[i,j]
    tomados desde u_xy, marcados con u_mask
    u_mask indica si el valor se ha calculado con edp.
    '''
    LHS = discreta.lhs # lado izquierdo de ecuacion
    RHS = discreta.rhs # lado derecho
    # sustituye U[i,j] con valores conocidos
    A_diagonal = [] # lista de i,j para matriz de coeficientes A
    # Separa términos suma
    term_suma = sym.Add.make_args(LHS)
    for term_k in term_suma:
        # busca U[i,j] y cambia por valor uxt[i,j]
        cambiar = 0 ; cambiar_valor = 0 ; cambiar_factor = 0
        # separa cada factor de término
        factor_Mul = sym.Mul.make_args(term_k)
        for factor_k in factor_Mul:
            # busca U[i,j] en matriz uxt[i,j]
            if factor_k.is_Function:
                [i_k,j_k] = factor_k.args
                if not(np.isnan(u_xy[i_k,j_k])):
                    cambiar = u_mask[i_k,j_k]
                    cambiar_factor = factor_k
                    cambiar_valor = u_xy[i_k,j_k]
                else:
                    A_diagonal.append([i_k,j_k,term_k/factor_k])
        # encontró valor U[i,j],term_k va a la derecha de ecuación
        if cambiar:
            LHS = LHS - term_k
            term_ki = term_k.subs(cambiar_factor,cambiar_valor)
            RHS = RHS - term_ki
    discreta = sym.Eq(LHS,RHS)
    B_diagonal = RHS
    resultado = [discreta,A_diagonal,B_diagonal]
    return (resultado)

# Desarrollo de iteraciones en cada nodo i,j, para la fila j_k
# genera el valor en cada punto
# Nota: Solo la primera fila, Tarea desarrollar para la matriz
print('\n iterar i,j:')
j_k = 1
for i_k in range(1,n-1,1):
    print('j: ',j_k,' ; ','i: ',i_k)
    discreta_ij = discreta.subs({i:i_k,j:j_k,
                                 x:xi[i_k],y:yj[j_k]})
    sym.pprint(discreta_ij)
    RHS = discreta_ij.rhs
    discreta_ij = sym.Eq(RHS,0)
    # usa valores de frontera segun u_mask con True
    discreta_k = edp_sustituyeValorU(discreta_ij,
                                    xi,yj,u_xy,u_mask)
    u_xy[i_k,j_k] = sym.Float(-discreta_k[2])
    print(buscar,'=',u_xy[i_k,j_k])
    
# SALIDA
np.set_printoptions(precision=verdigitos)
print('\nMétodo iterativo EDP Elíptica')
print('continuar el desarrollo con: ')
print('xi:',xi)
print('yj:',yj)
print(' j, U[i,j]')
for j_k in range(2,-1,-1):
    print(j_k, (u_xy[:,j_k]))

EDP Elípticas [ contínua a discreta ] || [ sympy_iterativo ] [ sympy_implícito ]
..


4. EDP Elípticas – Método Implícito con Sympy-Python

Desarrollo analítico del método implicito para una ecuación diferencial parcial elíptica usando Sympy. El algoritmo reutiliza el algoritmo para la EDP Elíptica de contínua a discreta, la creación de la matriz de valores u_xy, y la función
edp_sustituyeValorU() para buscar los valores conocidos de la u(x,y).

El algoritmo usa la ecuación discreta para en cada iteración i,j reemplazar los valores de U conocidos. Los valores de U se escriben en una matriz u_xy, para diferenciar si el valor existe se usa la matriz de estados u_mask.

Con las ecuaciones de cada iteración se llena la matriz A de coeficientes y el vector B de las constantes. Al resolver el sistema de ecuaciones se obtienen todos los valores de la matriz U, completando el ejercicio. Se desarrola la solución al sistema de ecuaciones usando Sympy como alternativa a usar numpy con np.linalg.solve(A.B) que se encuentra como comentario entre las instrucciones.

 discreta :
-4*U(i, j) + U(i, j - 1) + U(i, j + 1) + U(i - 1, j) + U(i + 1, j) = 0

Método Implícito - EDP Elíptica
j: 1 ; i: 1 ; ecuacion: 0
U(0, 1) + U(1, 0) - 4*U(1, 1) + U(1, 2) + U(2, 1) = 0
-4*U(1, 1) + U(1, 2) + U(2, 1) = -110.0
j: 1 ; i: 2 ; ecuacion: 1
U(1, 1) + U(2, 0) - 4*U(2, 1) + U(2, 2) + U(3, 1) = 0
U(1, 1) - 4*U(2, 1) + U(2, 2) + U(3, 1) = -50.0
j: 1 ; i: 3 ; ecuacion: 2
U(2, 1) + U(3, 0) - 4*U(3, 1) + U(3, 2) + U(4, 1) = 0
U(2, 1) - 4*U(3, 1) + U(3, 2) + U(4, 1) = -50.0
j: 1 ; i: 4 ; ecuacion: 3
U(3, 1) + U(4, 0) - 4*U(4, 1) + U(4, 2) + U(5, 1) = 0
U(3, 1) - 4*U(4, 1) + U(4, 2) + U(5, 1) = -50.0
j: 1 ; i: 5 ; ecuacion: 4
U(4, 1) + U(5, 0) - 4*U(5, 1) + U(5, 2) + U(6, 1) = 0
U(4, 1) - 4*U(5, 1) + U(5, 2) + U(6, 1) = -50.0
j: 1 ; i: 6 ; ecuacion: 5
U(5, 1) + U(6, 0) - 4*U(6, 1) + U(6, 2) + U(7, 1) = 0
U(5, 1) - 4*U(6, 1) + U(6, 2) + U(7, 1) = -50.0
j: 1 ; i: 7 ; ecuacion: 6
U(6, 1) + U(7, 0) - 4*U(7, 1) + U(7, 2) + U(8, 1) = 0
U(6, 1) - 4*U(7, 1) + U(7, 2) = -75.0
j: 2 ; i: 1 ; ecuacion: 7
U(0, 2) + U(1, 1) - 4*U(1, 2) + U(1, 3) + U(2, 2) = 0
U(1, 1) - 4*U(1, 2) + U(1, 3) + U(2, 2) = -60.0
j: 2 ; i: 2 ; ecuacion: 8
U(1, 2) + U(2, 1) - 4*U(2, 2) + U(2, 3) + U(3, 2) = 0
j: 2 ; i: 3 ; ecuacion: 9
U(2, 2) + U(3, 1) - 4*U(3, 2) + U(3, 3) + U(4, 2) = 0
j: 2 ; i: 4 ; ecuacion: 10
U(3, 2) + U(4, 1) - 4*U(4, 2) + U(4, 3) + U(5, 2) = 0
j: 2 ; i: 5 ; ecuacion: 11
U(4, 2) + U(5, 1) - 4*U(5, 2) + U(5, 3) + U(6, 2) = 0
j: 2 ; i: 6 ; ecuacion: 12
U(5, 2) + U(6, 1) - 4*U(6, 2) + U(6, 3) + U(7, 2) = 0
j: 2 ; i: 7 ; ecuacion: 13
U(6, 2) + U(7, 1) - 4*U(7, 2) + U(7, 3) + U(8, 2) = 0
U(6, 2) + U(7, 1) - 4*U(7, 2) + U(7, 3) = -25.0
j: 3 ; i: 1 ; ecuacion: 14
U(0, 3) + U(1, 2) - 4*U(1, 3) + U(1, 4) + U(2, 3) = 0
U(1, 2) - 4*U(1, 3) + U(1, 4) + U(2, 3) = -60.0
j: 3 ; i: 2 ; ecuacion: 15
U(1, 3) + U(2, 2) - 4*U(2, 3) + U(2, 4) + U(3, 3) = 0
j: 3 ; i: 3 ; ecuacion: 16
U(2, 3) + U(3, 2) - 4*U(3, 3) + U(3, 4) + U(4, 3) = 0
j: 3 ; i: 4 ; ecuacion: 17
U(3, 3) + U(4, 2) - 4*U(4, 3) + U(4, 4) + U(5, 3) = 0
j: 3 ; i: 5 ; ecuacion: 18
U(4, 3) + U(5, 2) - 4*U(5, 3) + U(5, 4) + U(6, 3) = 0
j: 3 ; i: 6 ; ecuacion: 19
U(5, 3) + U(6, 2) - 4*U(6, 3) + U(6, 4) + U(7, 3) = 0
j: 3 ; i: 7 ; ecuacion: 20
U(6, 3) + U(7, 2) - 4*U(7, 3) + U(7, 4) + U(8, 3) = 0
U(6, 3) + U(7, 2) - 4*U(7, 3) + U(7, 4) = -25.0
j: 4 ; i: 1 ; ecuacion: 21
U(0, 4) + U(1, 3) - 4*U(1, 4) + U(1, 5) + U(2, 4) = 0
U(1, 3) - 4*U(1, 4) + U(1, 5) + U(2, 4) = -60.0
j: 4 ; i: 2 ; ecuacion: 22
U(1, 4) + U(2, 3) - 4*U(2, 4) + U(2, 5) + U(3, 4) = 0
j: 4 ; i: 3 ; ecuacion: 23
U(2, 4) + U(3, 3) - 4*U(3, 4) + U(3, 5) + U(4, 4) = 0
j: 4 ; i: 4 ; ecuacion: 24
U(3, 4) + U(4, 3) - 4*U(4, 4) + U(4, 5) + U(5, 4) = 0
j: 4 ; i: 5 ; ecuacion: 25
U(4, 4) + U(5, 3) - 4*U(5, 4) + U(5, 5) + U(6, 4) = 0
j: 4 ; i: 6 ; ecuacion: 26
U(5, 4) + U(6, 3) - 4*U(6, 4) + U(6, 5) + U(7, 4) = 0
j: 4 ; i: 7 ; ecuacion: 27
U(6, 4) + U(7, 3) - 4*U(7, 4) + U(7, 5) + U(8, 4) = 0
U(6, 4) + U(7, 3) - 4*U(7, 4) + U(7, 5) = -25.0
j: 5 ; i: 1 ; ecuacion: 28
U(0, 5) + U(1, 4) - 4*U(1, 5) + U(1, 6) + U(2, 5) = 0
U(1, 4) - 4*U(1, 5) + U(2, 5) = -130.0
j: 5 ; i: 2 ; ecuacion: 29
U(1, 5) + U(2, 4) - 4*U(2, 5) + U(2, 6) + U(3, 5) = 0
U(1, 5) + U(2, 4) - 4*U(2, 5) + U(3, 5) = -70.0
j: 5 ; i: 3 ; ecuacion: 30
U(2, 5) + U(3, 4) - 4*U(3, 5) + U(3, 6) + U(4, 5) = 0
U(2, 5) + U(3, 4) - 4*U(3, 5) + U(4, 5) = -70.0
j: 5 ; i: 4 ; ecuacion: 31
U(3, 5) + U(4, 4) - 4*U(4, 5) + U(4, 6) + U(5, 5) = 0
U(3, 5) + U(4, 4) - 4*U(4, 5) + U(5, 5) = -70.0
j: 5 ; i: 5 ; ecuacion: 32
U(4, 5) + U(5, 4) - 4*U(5, 5) + U(5, 6) + U(6, 5) = 0
U(4, 5) + U(5, 4) - 4*U(5, 5) + U(6, 5) = -70.0
j: 5 ; i: 6 ; ecuacion: 33
U(5, 5) + U(6, 4) - 4*U(6, 5) + U(6, 6) + U(7, 5) = 0
U(5, 5) + U(6, 4) - 4*U(6, 5) + U(7, 5) = -70.0
j: 5 ; i: 7 ; ecuacion: 34
U(6, 5) + U(7, 4) - 4*U(7, 5) + U(7, 6) + U(8, 5) = 0
U(6, 5) + U(7, 4) - 4*U(7, 5) = -95.0

 A : 
 [[-4.  1.  0. ...  0.  0.  0.]
 [ 1. -4.  1. ...  0.  0.  0.]
 [ 0.  1. -4. ...  0.  0.  0.]
 ...
 [ 0.  0.  0. ... -4.  1.  0.]
 [ 0.  0.  0. ...  1. -4.  1.]
 [ 0.  0.  0. ...  0.  1. -4.]]

 B : 
 [-110.  -50.  -50.  -50.  -50.  -50.  -75.  -60.    0.    0.    0.    0.
    0.  -25.  -60.    0.    0.    0.    0.    0.  -25.  -60.    0.    0.
    0.    0.    0.  -25. -130.  -70.  -70.  -70.  -70.  -70.  -95.]
Resultados para U(x,y)
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 ]
 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.81 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.]
>>> 

Instrucciones en Python

Las instrucciones completas con Sympy-Python son:

# Ecuaciones Diferenciales Parciales Elipticas
# EDP Elípticas contínua a discreta con Sympy
import numpy as np
import sympy as sym

# funciones continuas y variables simbólicas usadas
y = sym.Symbol('y',real=True)
x = sym.Symbol('x',real=True)
u = sym.Function('u')(x,y)
f = sym.Function('f')(x,y) # funcion complemento
# funciones discretas y variables simbólicas usadas
i  = sym.Symbol('i',integer=True,positive=True)
j  = sym.Symbol('j',integer=True,positive=True)
Dx = sym.Symbol('Dx',real=True,positive=True)
Dy = sym.Symbol('Dy',real=True,positive=True)
L  = sym.Symbol('L',real=True)
U  = sym.Function('U')(i,j)

# INGRESO
fxy = 0*x+0*y  # f(x,y) = 0 , ecuacion de Poisson
# ecuacion edp : LHS=RHS
LHS = sym.diff(u,x,2) + sym.diff(u,y,2)
RHS = fxy
edp = sym.Eq(LHS-RHS,0)

dif_dividida ={sym.diff(u,x,2): (U.subs(i,i-1)-2*U+U.subs(i,i+1))/(Dx**2),
               sym.diff(u,y,2): (U.subs(j,j-1)-2*U+U.subs(j,j+1))/(Dy**2),
               sym.diff(u,y,1): (U - U.subs(j,j-1))/Dy}

# Condiciones iniciales en los bordes
fya = lambda y: 60 +0*y  # izquierda
fyb = lambda y: 25 +0*y  # derecha
fxc = lambda x: 50 +0*x  # inferior, función inicial 
fxd = lambda x: 70 +0*x  # superior, función inicial 

# 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
tolera = 0.0001

verdigitos = 2 # para mostrar en tabla de resultados

# PROCEDIMIENTO
def edp_discreta(edp,dif_dividida,x,y,u):
    ''' EDP contínua a discreta, usa diferencias divididas
        proporcionadas en parámetros, indica las variables x,y
        con función u de (x,y)
    '''
    resultado={}
    # expresión todo a la izquierda LHS (Left Hand side)
    LHS = edp.lhs
    RHS = edp.rhs
    if not(edp.rhs==0):
        LHS = LHS-RHS
        RHS = 0
        edp = sym.Eq(LHS,RHS)
    # para orden de derivada por x, y
    edp_x = edp.subs(x,0)
    edp_y = edp.subs(y,0)
    ordenDx = sym.ode_order(edp_x,u)
    ordenDy = sym.ode_order(edp_y,u)
    resultado['ordenDx'] = ordenDx
    resultado['ordenDy'] = ordenDy
    # convierte coeficiente d2u/dx2 a 1
    coeff_x = edp_coef_Dx(edp,x,ordenDx)
    if not(coeff_x==1):
        LHS = LHS/coeff_x
        RHS = RHS/coeff_x
    edp = sym.Eq(LHS,RHS)
    K_ = edp_coef_Dx(edp,y,ordenDy)
    resultado['edp=0']  = edp
    resultado['K_']  = K_
    # edp discreta
    discreta = edp.lhs 
    for derivada in dif_dividida:
        discreta = discreta.subs(derivada,dif_dividida[derivada])
    resultado['discreta=0'] = discreta
    return (resultado)

def edp_coef_Dx(edp,x,ordenx):
    ''' Extrae el coeficiente de la derivada Dx de ordenx
    '''
    coeff_x = 1.0
    # separa cada término de suma
    term_suma = sym.Add.make_args(edp.lhs)
    for term_k in term_suma:
        if term_k.is_Mul: # término de mas de un factor
            factor_Mul = sym.Mul.make_args(term_k)
            coeff_temp = 1; coeff_usar=0
            # separa cada factor de término 
            for factor_k in factor_Mul:
                if not(factor_k.is_Derivative):
                    coeff_temp = coeff_temp*factor_k
                else: # factor con derivada de ordenx
                    partes = factor_k.args
                    if partes[1]==(x,ordenx):
                        coeff_usar = 1
            if coeff_usar==1:
                coeff_x = coeff_x*coeff_temp
    return(coeff_x)

def edp_simplificaLamba(resultado,dx,dy):
    '''simplifica ecuacion con valores de lambda, dx y dy
    entregando la edp discreta simplificada
    '''
    discreta = resultado['discreta=0']
    ordenDy = resultado['ordenDy']
    ordenDx = resultado['ordenDx']
    K_ = resultado['K_']
    lamb = (Dy**ordenDy)/(Dx**ordenDx)
    if ordenDy==1 and ordenDx==2:
        lamb = lamb/K_
    resultado['Lambda_L'] = lamb
    # valor de Lambda en ecuacion edp
    L_k = lamb.subs([(Dx,dx),(Dy,dy),(1.0,1)])
    resultado['Lambda L_k'] = L_k
    # simplifica con lambda L
    discreta_L = sym.expand(discreta*(Dy**ordenDy),mul=True)
    resultado['(discreta=0)*Dy**ordeny'] = discreta_L
    discreta_L = edp_sustituye_L(resultado)
    discreta_L = discreta_L.subs(lamb,L)
    discreta_L = sym.collect(discreta_L,U)
    discreta_L = sym.Eq(discreta_L,0)
    resultado['discreta_L=0'] = discreta_L
    # sustituye constantes en ecuación a iterar
    discreta_L = discreta_L.subs([(Dx,dx),(Dy,dy),
                                  (L,L_k),(1.0,1)])
    resultado['discreta'] = discreta_L
    return(resultado)

def edp_sustituye_L(resultado):
    ''' sustituye lambda con Dy**ordeny/Dx**x/K_
    por L, al simplificar Lambda
    '''
    discreta = resultado['(discreta=0)*Dy**ordeny']
    ordenDy = resultado['ordenDy']
    ordenDx = resultado['ordenDx']
    discreta_L = 0
    # separa cada término de suma
    term_suma = sym.Add.make_args(discreta)
    for term_k in term_suma:
        # busca partes de L y cambia por valor L
        cambiar = 0 # por orden de derivada
        if term_k.has(Dx) and term_k.has(Dy):
            partes = term_k.args
            ordeny=1
            ordenx=1
            for unaparte in partes:    
                if unaparte.has(Dy):
                    if unaparte.is_Pow:
                        partey = unaparte.args
                        ordeny = partey[1]
                if unaparte.has(Dx):
                    if unaparte.is_Pow:
                        partey = unaparte.args
                        ordenx = partey[1]
            if (ordeny<=ordenDy and ordenx<=-ordenDx):
                cambiar=1
        if cambiar:
            term_k = term_k*L/resultado['Lambda_L']
        discreta_L = discreta_L + term_k
        # simplifica unos con decimal a entero 
        discreta_L = discreta_L.subs(1.0,1)
    return(discreta_L)

def edp_sustituyeValorU(discreta,xi,yj,u_xy,u_mask):
    '''Sustituye en edp discreta los valores conocidos de U[i,j]
    tomados desde u_xy, marcados con u_mask
    u_mask indica si el valor se ha calculado con edp.
    '''
    LHS = discreta.lhs # lado izquierdo de ecuacion
    RHS = discreta.rhs # lado derecho
    # sustituye U[i,j] con valores conocidos
    A_diagonal = [] # lista de i,j para matriz de coeficientes A
    # Separa términos suma
    term_suma = sym.Add.make_args(LHS)
    for term_k in term_suma:
        # busca U[i,j] y cambia por valor uxt[i,j]
        cambiar = 0 ; cambiar_valor = 0 ; cambiar_factor = 0
        # separa cada factor de término
        factor_Mul = sym.Mul.make_args(term_k)
        for factor_k in factor_Mul:
            # busca U[i,j] en matriz uxt[i,j]
            if factor_k.is_Function:
                [i_k,j_k] = factor_k.args
                if not(np.isnan(u_xy[i_k,j_k])):
                    cambiar = u_mask[i_k,j_k]
                    cambiar_factor = factor_k
                    cambiar_valor = u_xy[i_k,j_k]
                else:
                    A_diagonal.append([i_k,j_k,term_k/factor_k])
        # encontró valor U[i,j],term_k va a la derecha de ecuación
        if cambiar:
            LHS = LHS - term_k
            term_ki = term_k.subs(cambiar_factor,cambiar_valor)
            RHS = RHS - term_ki
    discreta = sym.Eq(LHS,RHS)
    B_diagonal = RHS
    resultado = [discreta,A_diagonal,B_diagonal]
    return (resultado)


# PROCEDIMIENTO
# transforma edp continua a discreta
resultado = edp_discreta(edp,dif_dividida,x,y,u)
resultado = edp_simplificaLamba(resultado,dx,dy)
discreta = resultado['discreta']

# SALIDA
np.set_printoptions(precision=verdigitos)
algun_numero = [int,float,str,'Lambda L_k']
print('EDP Elíptica contínua a discreta')
for entrada in resultado:
    tipo = type(resultado[entrada])
    if tipo in algun_numero or entrada in algun_numero:
        print('',entrada,':',resultado[entrada])
    else:
        print('\n',entrada,':')
        sym.pprint(resultado[entrada])

# ITERAR para cada i,j dentro de U ------------
# x[i] , y[j]  valor en posición en cada eje
xi = np.arange(x0,xn+dx,dx)
yj = np.arange(y0,yn+dy,dy)
n = len(xi)
m = len(yj)

# Matriz U
u_xy = np.zeros(shape=(n,m),dtype = float)
u_xy = u_xy*np.nan # valor inicial dentro de u
# llena u con valores en fronteras
u_xy[0,:]   = fya(yj)  # izquierda Ta
u_xy[n-1,:] = fyb(yj)  # derecha   Tb
u_xy[:,0]   = fxc(xi)  # inferior  Tc
u_xy[:,m-1] = fxd(xi)  # superior  Td
u0 = np.copy(u_xy)     # matriz u inicial

# u_mask[i,j] con valores iniciales o calculados:True
u_mask = np.zeros(shape=(n,m),dtype=bool)
u_mask[0,:]  = True # izquierda
u_mask[-1,:] = True # derecha
u_mask[:,0]  = True # inferior
u_mask[:,-1] = True # superior

# Método implícito para EDP Elíptica
# ITERAR para plantear las ecuaciones en [i,j]
resultado = {}
eq_itera = [] ; tamano = (n-2)*(m-2)
A = np.zeros(shape=(tamano,tamano),dtype=float)
B = np.zeros(tamano,dtype=float)
for j_k in range(1,m-1,1): # no usar valores en bordes
    for i_k in range(1,n-1,1): 
        eq_conteo = (j_k - 1)*(n-2)+(i_k-1)
        discreta_ij = discreta.subs({i:i_k,j:j_k,
                                     x:xi[i_k],y:yj[j_k]})
        resultado[eq_conteo]= {'j':j_k, 'i':i_k,
                                 'discreta_ij': discreta_ij}
        # usa valores de frontera segun u_mask con True
        discreta_k = edp_sustituyeValorU(discreta_ij,
                                        xi,yj,u_xy,u_mask)
        discreta_ij = discreta_k[0]
        A_diagonal  = discreta_k[1] # lista de (i,j,coeficiente) 
        B_diagonal  = discreta_k[2]
        resultado[eq_conteo]['discreta_k'] = discreta_k[0]
        # añade ecuacion a resolver
        eq_itera.append(discreta_ij)
        # Aplica coeficientes de ecuacion en A y B:
        # A_diagonal tiene lista de (i,j,coeficiente) 
        for uncoeff in A_diagonal:
            columna = (uncoeff[1]-1)*(n-2)+(uncoeff[0]-1)
            fila = (j_k - 1)*(n-2)+(i_k-1)
            A[fila,columna] = uncoeff[2]
        B[eq_conteo] = float(B_diagonal) # discreta_ij.rhs
resultado['A'] = np.copy(A)
resultado['B'] = np.copy(B)

# resuelve el sistema de ecuaciones en eq_itera en Sympy
X_k = sym.solve(eq_itera)[0]
# actualiza uxt[i,j] , u_mask segun X_k en Sympy
for nodo_Uij in X_k: 
    [i_k,j_k] = nodo_Uij.args
    u_xy[i_k,j_k] = X_k[nodo_Uij]
    u_mask[i_k,j_k] = True
# resuelve el sistema A.X=B en Numpy
#X = np.linalg.solve(A,B)
# tarea: llevar valores X a u_xy

# SALIDA
np.set_printoptions(precision=verdigitos)
algun_numero = [int,float,str,'Lambda L_k']
print('\nMétodo Implícito - EDP Elíptica')
for entrada in resultado:
    tipo = type(resultado[entrada])
    if tipo in algun_numero or entrada in algun_numero:
        print('',entrada,':',resultado[entrada])
    elif (tipo==dict):
        print('j:',resultado[entrada]['j'],'; '
              'i:',resultado[entrada]['i'],'; '
              'ecuacion:',entrada)
        sym.pprint(resultado[entrada]['discreta_ij'])
        if resultado[entrada]['discreta_k'].rhs!=0:
            sym.pprint(resultado[entrada]['discreta_k'])
    else:
        print('\n',entrada,': \n',resultado[entrada])
print('Resultados para U(x,y)')
print('xi:',xi)
print('yj:',yj)
print(' j, U[i,j]')
for j_k in range(m-1,-1,-1):
    print(j_k, (u_xy[:,j_k]))

EDP Elípticas [ contínua a discreta ] || [ sympy_iterativo ] [ sympy_implícito ]

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

EDP Parabólica [ contínua a discreta ] || [ sympy_explícito ] [sympy_ implícito ]
..


1. Ejercicio

Referencia:  Chapra 30.2 p888 pdf912, Burden 9Ed 12.2 p725, Rodríguez 10.2 p406

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

Valores de frontera: Ta = 60, Tb = 40, T0 = 25
Longitud en x a = 0, b = 1,  Constante K= 4
Tamaño de paso dx = 0.1, dt = dx/10

EDP Parabólica [ contínua a discreta ] || [ sympy_explícito ] [sympy_ implícito ]
..


2. EDP Parabólica contínua a discreta

El desarrollo analítico comienza en escribir la ecuación diferencial parcial parabólica de la forma contínua a la representación discreta, siguiendo los criterios indicados en la parte teórica usando derivadas expresadas en diferencias divididas.

La EDP se escribe en formato Sympy en el bloque de ingreso en dos partes: lado izquierdo (LHS) de la ecuación y lado derecho (RHS), definiendo u como una función de variables x,y.

# INGRESO
fxy = 0*x+0*y  # f(x,y) = 0 , ecuacion complementaria
# ecuacion edp : LHS=RHS
LHS = sym.diff(u,x,2) 
RHS = 4*sym.diff(u,y,1) + fxy
edp = sym.Eq(LHS-RHS,0)

La ecuación f(x,y) se añade para ser usada con otros ejercicios que puede revisar en la sección de Evaluaciones.

Las expresiones de diferencias divididas seleccionadas se ingresan como un diccionario:

dif_dividida ={sym.diff(u,x,2): (U.subs(i,i-1)-2*U+U.subs(i,i+1))/(Dx**2),
               sym.diff(u,y,1): (U.subs(j,j+1)-U)/Dy}

Las condiciones iniciales en los extremos a y b, y la temperatura en la barra, pueden ser funciones matemáticas, por lo que se propone usar el formato lambda de las variables x o y. En el ejercicio básico presentado en la parte teórica, la barra tiene temperatura constante, por lo que escribe como una  constante mas cero por la variable independiente, esta forma facilita la evaluación de un vector xi por ejemplo, en caso de disponer una expresión diferente.

# Valores de frontera
fya = lambda y: 60 +0*y  # izquierda
fyb = lambda y: 40 +0*y  # derecha
fxc = lambda x: 25 +0*x  # inferior, función inicial

los demás parámetros del ejercicio se ingresan más adelante de forma semejante al ejercicio con Numpy.

Todas las expresiones se escriben al lado izquierdo (LHS) de la igualdad, la parte del lado derecho(RHS) se la prefiere mantener en cero. La expresión se organiza manteniendo el coeficiente en 1 para el termino de Dx de mayor orden. Se sustituyen las derivadas por diferencias divididas, para luego simplificar la expresión usando λ como referencia.

Al desarrollar el método explícito se requiere indicar el nodo u(i,j+1) a buscar

buscar = U.subs(j,j+1) # método explícito

El valor de K se puede determinar en la expresión al reordenar para que el término de Dx de mayor grado sea la unidad. El valor de K será el del término Dy con signo negativo, al encontrarse en el lado izquierdo (LHS) de la expresión, en la propuesta teórica se encontraba como positivo del lado derecho.

El término de mayor grado se puede obtener al revisar la expresión por cada término de la suma, buscando el factor que contiene Dx o Dy en cada caso. Por ser un poco ordenado se añade la función edp_coef_Dx().

El resultado del algoritmo para el ejercicio propuesto es:

EDP Parabólica contínua a discreta
 ordenDx : 2
 ordenDy : 1

 edp=0 :
  2                             
 d               d              
---(u(x, y)) - 4*--(u(x, y)) = 0
  2              dy             
dx                              
 K_ : 4.0

 Lambda L :
0.25*Dy
-------
    2  
  Dx   
 Lambda L_k : 0.24999999999999997

 (discreta=0)*Dy**ordeny :
                            2*Dy*U(i, j)   Dy*U(i - 1, j)   Dy*U(i + 1, j)
4*U(i, j) - 4*U(i, j + 1) - ------------ + -------------- + --------------
                                  2               2                2      
                                Dx              Dx               Dx       

 discreta_L=0 :
/    2*Dy\                           Dy*U(i - 1, j)   Dy*U(i + 1, j)    
|4 - ----|*U(i, j) - 4*U(i, j + 1) + -------------- + -------------- = 0
|      2 |                                  2                2          
\    Dx  /                                Dx               Dx           

 discreta :
2.0*U(i, j) - 4*U(i, j + 1) + 1.0*U(i - 1, j) + 1.0*U(i + 1, j) = 0

>>>

Instrucciones en Python

Para revisar los resultados de algunos pasos para obtener la expresión discreta de EDP, se usa un diccionario a ser usado en el bloque de salida. Para las expresiones de ecuación se usa sym.pprint().

# Ecuaciones Diferenciales Parciales Parabólicas
# EDP Parabólicas contínua a discreta con Sympy
import numpy as np
import sympy as sym

# funciones continuas y variables simbólicas usadas
t = sym.Symbol('t',real=True)
y = sym.Symbol('y',real=True)
x = sym.Symbol('x',real=True)
u = sym.Function('u')(x,y)
K = sym.Symbol('K',real=True)
f = sym.Function('f')(x,y) # funcion complemento
# funciones discretas y variables simbólicas usadas
i  = sym.Symbol('i',integer=True,positive=True)
j  = sym.Symbol('j',integer=True,positive=True)
Dt = sym.Symbol('Dt',real=True,positive=True)
Dx = sym.Symbol('Dx',real=True,positive=True)
Dy = sym.Symbol('Dy',real=True,positive=True)
L  = sym.Symbol('L',real=True)
U  = sym.Function('U')(i,j)

# INGRESO
fxy = 0*x+0*y  # f(x,y) = 0 , ecuacion complementaria
# ecuacion edp : LHS=RHS
LHS = sym.diff(u,x,2) 
RHS = 4*sym.diff(u,y,1) + fxy
edp = sym.Eq(LHS-RHS,0)

dif_dividida ={sym.diff(u,x,2): (U.subs(i,i-1)-2*U+U.subs(i,i+1))/(Dx**2),
               sym.diff(u,y,1): (U.subs(j,j+1)-U)/Dy}

buscar = U.subs(j,j+1) # método explícito

# Valores de frontera
fya = lambda y: 60 +0*y  # izquierda
fyb = lambda y: 40 +0*y  # derecha
fxc = lambda x: 25 +0*x  # inferior, función inicial

# dimensiones de la barra
a = 0  # longitud en x
b = 1
y0 = 0 # tiempo inicial, aumenta con dt en n iteraciones

# discretiza, supone dx=dy
x_tramos = 10
dx  = (b-a)/x_tramos
dy  = dx/10

verdigitos = 2 # para mostrar en tabla de resultados

n_iteraciones = 4 # iteraciones en tiempo

# PROCEDIMIENTO

def edp_coef_Dx(edp,x,ordenx):
    ''' Extrae el coeficiente de la derivada Dx de ordenx
    '''
    coeff_x = 1.0
    # separa cada término de suma
    term_suma = sym.Add.make_args(edp.lhs)
    for term_k in term_suma:
        if term_k.is_Mul: # mas de un factor
            factor_Mul = sym.Mul.make_args(term_k)
            coeff_temp = 1; coeff_usar=0
            # separa cada factor de término 
            for factor_k in factor_Mul:
                if not(factor_k.is_Derivative):
                    coeff_temp = coeff_temp*factor_k
                else: # factor con derivada de ordenx
                    partes = factor_k.args
                    if partes[1]==(x,ordenx):
                        coeff_usar = 1
            if coeff_usar==1:
                coeff_x = coeff_x*coeff_temp
    return(coeff_x)

resultado={}
# expresión todo a la izquierda LHS (Left Hand side)
if not(edp.rhs==0):
    LHS = EDP.lhs
    RHS = EDP.rhs
    edp = sym.Eq(LHS,RHS)
# para orden de derivada por x, y
edp_x = edp.subs(x,0)
edp_y = edp.subs(y,0)
ordenDx = sym.ode_order(edp_x,u)
ordenDy = sym.ode_order(edp_y,u)
resultado['ordenDx'] = ordenDx
resultado['ordenDy'] = ordenDy
# convierte coeficiente d2u/dx2 a 1
coeff_x = edp_coef_Dx(edp,x,ordenDx)
if not(coeff_x==1):
    LHS = LHS/coeff_x
    RHS = RHS/coeff_x
edp = sym.Eq(LHS-RHS,0)
K_ = -edp_coef_Dx(edp,y,ordenDy)
resultado['edp=0']  = edp
resultado['K_']  = float(K_)

lamb = (Dy**ordenDy)/(Dx**ordenDx)
if ordenDy==1 and ordenDx==2:
        lamb = lamb/K_
resultado['Lambda_L'] = lamb
# valor de Lambda en ecuacion edp
L_k = lamb.subs([(Dx,dx),(Dy,dy),(1.0,1)])
resultado['Lambda L_k'] = float(L_k)

discreta = edp.lhs # EDP discreta
for derivada in dif_dividida:
    discreta = discreta.subs(derivada,dif_dividida[derivada])

# simplifica con Lambda L
discreta_L = sym.expand(discreta*(Dy**ordenDy),mul=True)
resultado['(discreta=0)*Dy**ordeny'] = discreta_L

discreta_L = discreta_L.subs(lamb,L)
discreta_L = sym.collect(discreta_L,U)
discreta_L = sym.Eq(discreta_L,0)
resultado['discreta_L=0'] = discreta_L
        
# sustituye constantes en ecuación a iterar
discreta_L = discreta_L.subs([(Dx,dx),(Dy,dy),
                              (L,L_k),(1.0+0*x,1)])
discreta_L = discreta_L.subs(1.0+0*x,1)
resultado['discreta'] = discreta_L

# SALIDA
np.set_printoptions(precision=verdigitos)
algun_numero = [int,float,str,'Lambda L_k']
print('EDP Parabólica contínua a discreta')
for entrada in resultado:
    tipo = type(resultado[entrada])
    if tipo in algun_numero or entrada in algun_numero:
        print('',entrada,':',resultado[entrada])
    else:
        print('\n',entrada,':')
        sym.pprint(resultado[entrada])

EDP Parabólica [ contínua a discreta ] || [ sympy_explícito ] [sympy_ implícito ]

..


3. EDP Parabólica – Método explícito con Sympy

El desarrollo del método explícito para una ecuación diferencial parcial parabólica usando Sympy, requiere añadir instrucciones para procesar la expresión obtenida en el algoritmo anterior con los valores de i,j.

Se añaden las instrucciones para iterar la expresión discreta para cada i,j dentro de la matriz U creada para el ejercicio. Se aplican las condiciones iniciales a los bordes de la matriz, para luego proceder a iterar.

Como en las actividades del curso se requiere realizar al menos 3 iteraciones para las expresiones del algoritmo con papel y lápiz, solo se presentarán los resultados de las expresiones para la primera fila en j=0.

El algoritmo para la EDP del ejemplo muestra el siguiente resultado:

Método explícito EDP Parabólica
 discreta_itera : U(i, j + 1) = 0.5*U(i, j) + 0.25*U(i - 1, j) + 0.25*U(i + 1, j) 
 iterar i,j:
j:  0  ;  i:  1  ;  ecuacion: 0
U(1, 1) = 0.25*U(0, 0) + 0.5*U(1, 0) + 0.25*U(2, 0)
U(1, 1) = 33.75
j:  0  ;  i:  2  ;  ecuacion: 1
U(2, 1) = 0.25*U(1, 0) + 0.5*U(2, 0) + 0.25*U(3, 0)
U(2, 1) = 25.0
j:  0  ;  i:  3  ;  ecuacion: 2
U(3, 1) = 0.25*U(2, 0) + 0.5*U(3, 0) + 0.25*U(4, 0)
U(3, 1) = 25.0
j:  0  ;  i:  4  ;  ecuacion: 3
U(4, 1) = 0.25*U(3, 0) + 0.5*U(4, 0) + 0.25*U(5, 0)
U(4, 1) = 25.0
j:  0  ;  i:  5  ;  ecuacion: 4
U(5, 1) = 0.25*U(4, 0) + 0.5*U(5, 0) + 0.25*U(6, 0)
U(5, 1) = 25.0
j:  0  ;  i:  6  ;  ecuacion: 5
U(6, 1) = 0.25*U(5, 0) + 0.5*U(6, 0) + 0.25*U(7, 0)
U(6, 1) = 25.0
j:  0  ;  i:  7  ;  ecuacion: 6
U(7, 1) = 0.25*U(6, 0) + 0.5*U(7, 0) + 0.25*U(8, 0)
U(7, 1) = 25.0
j:  0  ;  i:  8  ;  ecuacion: 7
U(8, 1) = 0.25*U(7, 0) + 0.5*U(8, 0) + 0.25*U(9, 0)
U(8, 1) = 25.0
j:  0  ;  i:  9  ;  ecuacion: 8
U(9, 1) = 0.25*U(8, 0) + 0.5*U(9, 0) + 0.25*U(10, 0)
U(9, 1) = 28.75
... continua calculando

x: [0.  0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1. ]
t: [0.   0.01 0.02 0.03 0.04]
 j, U[i,j]
4 [60.   42.77 31.29 26.37 25.14 25.   25.06 25.59 27.7  32.62 40.  ]
3 [60.   40.86 29.38 25.55 25.   25.   25.   25.23 26.88 31.8  40.  ]
2 [60.   38.12 27.19 25.   25.   25.   25.   25.   25.94 30.62 40.  ]
1 [60.   33.75 25.   25.   25.   25.   25.   25.   25.   28.75 40.  ]
0 [60. 25. 25. 25. 25. 25. 25. 25. 25. 25. 40.]
>>> 

Instrucciones en Python

# Método explícito EDP Parabólica ------------
# ITERAR para cada i,j dentro de U 
# x[i] , y[j]  valor en posición en cada eje
xi = np.arange(a,b+dx,dx)
yj = np.arange(y0,n_iteraciones*dy+dy,dy)
n = len(xi)
m = len(yj)
# Matriz U
u_xy = np.zeros(shape=(n,m),dtype = float)
u_xy = u_xy*np.nan # valor inicial dentro de u
# llena u con valores en fronteras
u_xy[:,0]   = fxc(xi)  # inferior  Tc
u_xy[0,:]   = fya(yj)  # izquierda Ta
u_xy[n-1,:] = fyb(yj)  # derecha   Tb
u0 = np.copy(u_xy)     # matriz u inicial

# busca el término no conocido
u_explicito = sym.solve(discreta_L,buscar)
u_explicito = sym.Eq(buscar,u_explicito[0])
resultado['discreta_itera'] = u_explicito

print('\nMétodo explícito EDP Parabólica')
print('discreta_itera')
sym.pprint(u_explicito)
# iterando
print('\n iterar i,j:')
j_k = 0
for j_k in range(0,m-1,1):
    for i_k in range(1,n-1,1):
        eq_conteo = j_k*(n-2)+(i_k-1)
        discreta_ij = u_explicito.subs({i:i_k,j:j_k})
        valorij = 0 # calcula valor de nodo ij
        # separa cada término de suma
        term_suma = sym.Add.make_args(discreta_ij.rhs)
        for term_k in term_suma:
            term_ki = 1
            if term_k.is_Function:
                [i_f,j_f] = term_k.args
                term_ki = u_xy[i_f,j_f]
            elif term_k.is_Mul: # mas de un factor
                factor_Mul = sym.Mul.make_args(term_k)
                # separa cada factor de término 
                for factor_k in factor_Mul:
                    if factor_k.is_Function:
                        [i_f,j_f] = factor_k.args
                        term_ki = term_ki*u_xy[i_f,j_f]
                    else:
                        term_ki = term_ki*factor_k
            else:
                term_ki = term_k
            valorij = valorij + term_ki
        [i_f,j_f] = discreta_ij.lhs.args
        u_xy[i_f,j_f] = float(valorij) # actualiza matriz u
        # muestra las primeras m iteraciones
        if eq_conteo<(n-2):
            print('j: ',j_k,' ; ','i: ',i_k,
                  ' ; ','ecuacion:',eq_conteo)
            sym.pprint(discreta_ij)
            print(discreta_ij.lhs,'=',float(valorij))
print('... continua calculando')

# SALIDA
np.set_printoptions(precision=verdigitos)
print('\nx:',xi)
print('t:',yj)
print(' j, U[i,j]')
for j_k in range(m-1,-1,-1):
    print(j_k, (u_xy[:,j_k]))

EDP Parabólica [ contínua a discreta ] || [ sympy_explícito ] [sympy_ implícito ]
..


4. EDP Parabólica – Método implícito con Sympy

Para el ejercicio presentado en el numeral 1, se desarrolla el método implícito.

Las diferencias finitas centradas y hacia atras se definen en el diccionario dif_dividida.

Para desarrollar la expresión discreta, se incorpora la función edp_sustituye_L(), que busca el en cada término de la suma los componentes de λ y los sustituye por la variable L.

Para cada valor de j, se crean las ecuaciones desde la forma discreta de la EDP, moviendo a la derecha los valores conocidos para generar un sistema de ecuaciones A.X=B. Se resuelve el sistema de ecuaciones  y se actualiza la matriz u[i,j] de resultados.

Por la extensión de los pasos a mostrar, para las iteraciones en cada fila, solo se presentan para la primera fila. En caso de requerir observar más filas, puede cambiar el parámetro en el bloque de salida, en el condicional:

    elif (tipo==dict):
        if entrada<2

Los resultados para el algoritmo son:

EDP Parabólica contínua a discreta
 ordenDx : 2
 ordenDy : 1

 edp=0 :
  2                             
 d               d              
---(u(x, y)) - 4*--(u(x, y)) = 0
  2              dy             
dx                              
 K_ : 4.0

 Lambda_L :
0.25*Dy
-------
    2  
  Dx   
 Lambda L_k : 0.24999999999999997

 (discreta=0)*Dy**ordeny :
                             2*Dy*U(i, j)   Dy*U(i - 1, j)   Dy*U(i + 1, j)
-4*U(i, j) + 4*U(i, j - 1) - ------------ + -------------- + --------------
                                   2               2                2      
                                 Dx              Dx               Dx       

 discreta_L=0 :
4.0*L*U(i - 1, j) + 4.0*L*U(i + 1, j) + (-8.0*L - 4)*U(i, j) + 4*U(i, j - 1) =
 0

 discreta :
-6.0*U(i, j) + 4*U(i, j - 1) + 1.0*U(i - 1, j) + 1.0*U(i + 1, j) = 0

Método implícito EDP Parabólica
j: 1 ; i: 1
1.0*U(0, 1) + 4*U(1, 0) - 6.0*U(1, 1) + 1.0*U(2, 1) = 0
-6.0*U(1, 1) + 1.0*U(2, 1) = -160.0
j: 1 ; i: 2
1.0*U(1, 1) + 4*U(2, 0) - 6.0*U(2, 1) + 1.0*U(3, 1) = 0
1.0*U(1, 1) - 6.0*U(2, 1) + 1.0*U(3, 1) = -100.0
j: 1 ; i: 3
1.0*U(2, 1) + 4*U(3, 0) - 6.0*U(3, 1) + 1.0*U(4, 1) = 0
1.0*U(2, 1) - 6.0*U(3, 1) + 1.0*U(4, 1) = -100.0
j: 1 ; i: 4
1.0*U(3, 1) + 4*U(4, 0) - 6.0*U(4, 1) + 1.0*U(5, 1) = 0
1.0*U(3, 1) - 6.0*U(4, 1) + 1.0*U(5, 1) = -100.0
j: 1 ; i: 5
1.0*U(4, 1) + 4*U(5, 0) - 6.0*U(5, 1) + 1.0*U(6, 1) = 0
1.0*U(4, 1) - 6.0*U(5, 1) + 1.0*U(6, 1) = -100.0
j: 1 ; i: 6
1.0*U(5, 1) + 4*U(6, 0) - 6.0*U(6, 1) + 1.0*U(7, 1) = 0
1.0*U(5, 1) - 6.0*U(6, 1) + 1.0*U(7, 1) = -100.0
j: 1 ; i: 7
1.0*U(6, 1) + 4*U(7, 0) - 6.0*U(7, 1) + 1.0*U(8, 1) = 0
1.0*U(6, 1) - 6.0*U(7, 1) + 1.0*U(8, 1) = -100.0
j: 1 ; i: 8
1.0*U(7, 1) + 4*U(8, 0) - 6.0*U(8, 1) + 1.0*U(9, 1) = 0
1.0*U(7, 1) - 6.0*U(8, 1) + 1.0*U(9, 1) = -100.0
j: 1 ; i: 9
1.0*U(8, 1) + 4*U(9, 0) - 6.0*U(9, 1) + 1.0*U(10, 1) = 0
1.0*U(8, 1) - 6.0*U(9, 1) = -140.0
A :
[[-6.  1.  0.  0.  0.  0.  0.  0.  0.]
 [ 1. -6.  1.  0.  0.  0.  0.  0.  0.]
 [ 0.  1. -6.  1.  0.  0.  0.  0.  0.]
 [ 0.  0.  1. -6.  1.  0.  0.  0.  0.]
 [ 0.  0.  0.  1. -6.  1.  0.  0.  0.]
 [ 0.  0.  0.  0.  1. -6.  1.  0.  0.]
 [ 0.  0.  0.  0.  0.  1. -6.  1.  0.]
 [ 0.  0.  0.  0.  0.  0.  1. -6.  1.]
 [ 0.  0.  0.  0.  0.  0.  0.  1. -6.]]
B :
[-160. -100. -100. -100. -100. -100. -100. -100. -140.]
X :
[31.01 26.03 25.18 25.03 25.01 25.01 25.08 25.44 27.57]
j: 2 ; i: 1
1.0*U(0, 2) + 4*U(1, 1) - 6.0*U(1, 2) + 1.0*U(2, 2) = 0
-6.0*U(1, 2) + 1.0*U(2, 2) = -184.020210038146
...
Resultados para U(x,y)
xi: [0.  0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1. ]
yj: [0.   0.01 0.02 0.03 0.04]
 j, U[i,j]
4 [60.   40.67 30.61 26.75 25.51 25.18 25.24 25.76 27.41 31.71 40.  ]
3 [60.   38.34 29.06 26.09 25.28 25.09 25.13 25.47 26.74 30.72 40.  ]
2 [60.   35.25 27.49 25.55 25.12 25.03 25.05 25.24 26.07 29.39 40.  ]
1 [60.   31.01 26.03 25.18 25.03 25.01 25.01 25.08 25.44 27.57 40.  ]
0 [60. 25. 25. 25. 25. 25. 25. 25. 25. 25. 40.]
>>>

Instrucciones en Python

Para el método implícito se añade la sección donde se reemplazan los valores en los bordes y en la barra en la matriz U(x,y) en el algoritmo identificada como  u_xy . Para diferenciar si el valor existe se usa la matriz de estados u_mask con estados True/False.

También se incorpora la función edp_sustituyeValorU() que para cada término suma de la ecuación busca en la matriz u_xy si el valor U(i,j) existe y lo cambia. Luego pasa el valor al lado derecho de la ecuación par formar el sistema de ecuaciones.

Con los valores existentes antes de la iteración, se obtienen las ecuaciones por cada punto i,j de una fila para el sistema de ecuaciones en la forma A.x=B. El sistema se resuelve con las instrucciones de Numpy. El resultado se actualiza en la matriz u_xy y la matriz  u_mask para indicar que el valor ya está disponible.

Compare los resultados con el método numérico y mida los tiempos de ejecución de cada algoritmo para escribir sus observaciones.

# Ecuaciones Diferenciales Parciales Parabólicas
# EDP Parabólicas contínua a discreta con Sympy
import numpy as np
import sympy as sym

# funciones continuas y variables simbólicas usadas
t = sym.Symbol('t',real=True)
y = sym.Symbol('y',real=True)
x = sym.Symbol('x',real=True)
u = sym.Function('u')(x,y)
K = sym.Symbol('K',real=True)
f = sym.Function('f')(x,y) # funcion complemento
# funciones discretas y variables simbólicas usadas
i  = sym.Symbol('i',integer=True,positive=True)
j  = sym.Symbol('j',integer=True,positive=True)
Dt = sym.Symbol('Dt',real=True,positive=True)
Dx = sym.Symbol('Dx',real=True,positive=True)
Dy = sym.Symbol('Dy',real=True,positive=True)
L  = sym.Symbol('L',real=True)
U  = sym.Function('U')(i,j)

# INGRESO
fxy = 0*x+0*y  # f(x,y) = 0 , ecuacion complementaria
# ecuacion edp : LHS=RHS
LHS = sym.diff(u,x,2) 
RHS = 4*sym.diff(u,y,1) + fxy
edp = sym.Eq(LHS-RHS,0)

dif_dividida ={sym.diff(u,x,2): (U.subs(i,i-1)-2*U+U.subs(i,i+1))/(Dx**2),
               sym.diff(u,y,1): (U - U.subs(j,j-1))/Dy}

# Valores de frontera
fya = lambda y: 60 +0*y  # izquierda
fyb = lambda y: 40 +0*y  # derecha
fxc = lambda x: 25 +0*x  # inferior, función inicial

# dimensiones de la barra
a = 0  # longitud en x
b = 1
y0 = 0 # tiempo inicial, aumenta con dt en n iteraciones

# discretiza
x_tramos = 10
dx  = (b-a)/x_tramos
dy  = dx/10

verdigitos = 2 # para mostrar en tabla de resultados

n_iteraciones = 4 # iteraciones en tiempo

# PROCEDIMIENTO

def edp_coef_Dx(edp,x,ordenx):
    ''' Extrae el coeficiente de la derivada Dx de ordenx
    '''
    coeff_x = 1.0
    # separa cada término de suma
    term_suma = sym.Add.make_args(edp.lhs)
    for term_k in term_suma:
        if term_k.is_Mul: # mas de un factor
            factor_Mul = sym.Mul.make_args(term_k)
            coeff_temp = 1; coeff_usar=0
            # separa cada factor de término 
            for factor_k in factor_Mul:
                if not(factor_k.is_Derivative):
                    coeff_temp = coeff_temp*factor_k
                else: # factor con derivada de ordenx
                    partes = factor_k.args
                    if partes[1]==(x,ordenx):
                        coeff_usar = 1
            if coeff_usar==1:
                coeff_x = coeff_x*coeff_temp
    return(coeff_x)

def edp_sustituye_L(resultado):
    ''' sustituye lambda con Dy**ordeny/Dx**x/K_
    por L, al simplificar Lambda
    '''
    discreta = resultado['(discreta=0)*Dy**ordeny']
    ordenDy = resultado['ordenDy']
    ordenDx = resultado['ordenDx']
    discreta_L = 0
    # separa cada término de suma
    term_suma = sym.Add.make_args(discreta)
    for term_k in term_suma:
        # busca partes de L y cambia por valor L
        cambiar = 0 # por orden de derivada
        if term_k.has(Dx) and term_k.has(Dy):
            partes = term_k.args
            ordeny=1
            ordenx=1
            for unaparte in partes:    
                if unaparte.has(Dy):
                    if unaparte.is_Pow:
                        partey = unaparte.args
                        ordeny = partey[1]
                if unaparte.has(Dx):
                    if unaparte.is_Pow:
                        partey = unaparte.args
                        ordenx = partey[1]
            if (ordeny<=ordenDy and ordenx<=-ordenDx):
                cambiar=1
        if cambiar:
            term_k = term_k*L/resultado['Lambda_L']
        discreta_L = discreta_L + term_k
        # simplifica unos con decimal a entero 
        discreta_L = discreta_L.subs(1.0,1)
    return(discreta_L)

resultado={}
# expresión todo a la izquierda LHS (Left Hand side)
if not(edp.rhs==0):
    LHS = EDP.lhs
    RHS = EDP.rhs
    edp = sym.Eq(LHS,RHS)
# para orden de derivada por x, y
edp_x = edp.subs(x,0)
edp_y = edp.subs(y,0)
ordenDx = sym.ode_order(edp_x,u)
ordenDy = sym.ode_order(edp_y,u)
resultado['ordenDx'] = ordenDx
resultado['ordenDy'] = ordenDy
# convierte coeficiente d2u/dx2 a 1
coeff_x = edp_coef_Dx(edp,x,ordenDx)
if not(coeff_x==1):
    LHS = LHS/coeff_x
    RHS = RHS/coeff_x
edp = sym.Eq(LHS-RHS,0)
K_ = -edp_coef_Dx(edp,y,ordenDy)
resultado['edp=0']  = edp
resultado['K_']  = float(K_)

lamb = (Dy**ordenDy)/(Dx**ordenDx)
if ordenDy==1 and ordenDx==2:
        lamb = lamb/K_
resultado['Lambda_L'] = lamb
# valor de Lambda en ecuacion edp
L_k = lamb.subs([(Dx,dx),(Dy,dy),(1.0,1)])
resultado['Lambda L_k'] = float(L_k)

discreta = edp.lhs # EDP discreta
for derivada in dif_dividida:
    discreta = discreta.subs(derivada,dif_dividida[derivada])

# simplifica con Lambda L
discreta_L = sym.expand(discreta*(Dy**ordenDy),mul=True)
resultado['(discreta=0)*Dy**ordeny'] = discreta_L

discreta_L = edp_sustituye_L(resultado)
discreta_L = sym.collect(discreta_L,U)
discreta_L = sym.Eq(discreta_L,0)
resultado['discreta_L=0'] = discreta_L
        
# sustituye constantes en ecuación a iterar
discreta_L = discreta_L.subs([(Dx,dx),(Dy,dy),
                              (L,L_k),(1.0+0*x,1)])
resultado['discreta'] = discreta_L

# SALIDA
np.set_printoptions(precision=verdigitos)
algun_numero = [int,float,str,'Lambda L_k']
print('EDP Parabólica contínua a discreta')
for entrada in resultado:
    tipo = type(resultado[entrada])
    if tipo in algun_numero or entrada in algun_numero:
        print('',entrada,':',resultado[entrada])
    else:
        print('\n',entrada,':')
        sym.pprint(resultado[entrada])

# Método implícito EDP Parabólica ------------
# ITERAR para cada i,j dentro de U 
# x[i] , y[j]  valor en posición en cada eje
xi = np.arange(a,b+dx,dx)
yj = np.arange(y0,n_iteraciones*dy+dy,dy)
n = len(xi)
m = len(yj)
# Matriz U
u_xy = np.zeros(shape=(n,m),dtype = float)
u_xy = u_xy*np.nan # valor inicial dentro de u
# llena u con valores en fronteras
u_xy[:,0]   = fxc(xi)  # inferior  Tc
u_xy[0,:]   = fya(yj)  # izquierda Ta
u_xy[n-1,:] = fyb(yj)  # derecha   Tb
u0 = np.copy(u_xy)     # matriz u inicial

# u_mask[i,j] con valores iniciales o calculados:True
u_mask = np.zeros(shape=(n,m),dtype=bool)
u_mask[0,:]  = True # izquierda Ta
u_mask[-1,:] = True # derecha   Tb
u_mask[:,0]  = True # inferior  Tc

# Método implícito para EDP Parabólica
# a usar en iteraciones
discreta = resultado['discreta']

def edp_sustituyeValorU(discreta,xi,yj,u_xy,u_mask):
    '''Sustituye en edp discreta los valores conocidos de U[i,j]
    tomados desde u_xy, marcados con u_mask
    u_mask indica si el valor se ha calculado con edp.
    '''
    LHS = discreta.lhs # lado izquierdo de ecuacion
    RHS = discreta.rhs # lado derecho
    # sustituye U[i,j] con valores conocidos
    A_diagonal = [] # lista de i,j para matriz de coeficientes A
    # Separa términos suma
    term_suma = sym.Add.make_args(LHS)
    for term_k in term_suma:
        # busca U[i,j] y cambia por valor uxt[i,j]
        cambiar = 0 ; cambiar_valor = 0 ; cambiar_factor = 0
        # separa cada factor de término
        factor_Mul = sym.Mul.make_args(term_k)
        for factor_k in factor_Mul:
            # busca U[i,j] en matriz uxt[i,j]
            if factor_k.is_Function:
                [i_k,j_k] = factor_k.args
                if not(np.isnan(u_xy[i_k,j_k])):
                    cambiar = u_mask[i_k,j_k]
                    cambiar_factor = factor_k
                    cambiar_valor = u_xy[i_k,j_k]
                else:
                    A_diagonal.append([i_k,j_k,term_k/factor_k])
        # encontró valor U[i,j],term_k va a la derecha de ecuación
        if cambiar:
            LHS = LHS - term_k
            term_ki = term_k.subs(cambiar_factor,cambiar_valor)
            RHS = RHS - term_ki
    discreta = sym.Eq(LHS,RHS)
    B_diagonal = RHS
    resultado = [discreta,A_diagonal,B_diagonal]
    return (resultado)

# ITERAR para plantear las ecuaciones en [i,j]
resultado = {}
eq_itera = [] ; tamano = (n-2)
A = np.zeros(shape=(tamano,tamano),dtype=float)
B = np.zeros(tamano,dtype=float)
for j_k in range(1,m,1): # no usar valores en bordes
    resultado[j_k] ={}
    for i_k in range(1,n-1,1): 
        eq_conteo = i_k-1
        discreta_ij = discreta.subs({i:i_k,j:j_k,
                                     x:xi[i_k],y:yj[j_k]})
        resultado[j_k][eq_conteo]= {'j':j_k, 'i':i_k,
                                 'discreta_ij': discreta_ij}
        # usa valores de frontera segun u_mask con True
        discreta_k = edp_sustituyeValorU(discreta_ij,
                                        xi,yj,u_xy,u_mask)
        discreta_ij = discreta_k[0]
        A_diagonal  = discreta_k[1] # lista de (i,j,coeficiente) 
        B_diagonal  = discreta_k[2] # constante para B
        resultado[j_k][eq_conteo]['discreta_k'] = discreta_k[0]
        # añade ecuacion a resolver
        eq_itera.append(discreta_ij)
        # Aplica coeficientes de ecuacion en A y B:
        # A_diagonal tiene lista de (i,j,coeficiente) 
        for uncoeff in A_diagonal:
            columna = uncoeff[0]-1
            fila = eq_conteo
            A[fila,columna] = uncoeff[2]
        B[eq_conteo] = float(B_diagonal) # discreta_ij.rhs
    # resuelve el sistema A.X=B en Numpy
    X_k = np.linalg.solve(A,B)
    # actualiza uxt[i,j] , u_mask segun X_k
    u_xy[1:n-1,j_k] = X_k
    u_mask[1:n-1,j_k] = True
    # almacena resultados
    resultado[j_k]['A'] = np.copy(A)
    resultado[j_k]['B'] = np.copy(B)
    resultado[j_k]['X'] = np.copy(X_k)

# SALIDA
np.set_printoptions(precision=verdigitos)
algun_numero = [int,float,str,'Lambda L_k']
print('\nMétodo implícito EDP Parabólica')
for entrada in resultado:
    tipo = type(resultado[entrada])
    if tipo in algun_numero or entrada in algun_numero:
        print('',entrada,':',resultado[entrada])
    elif (tipo==dict):
        if entrada<2:
            eq_conteo = resultado[entrada].keys()
            for una_eq in eq_conteo:
                if una_eq=='A' or una_eq=='B' or una_eq=='X':
                    print(una_eq,':')
                    print(resultado[entrada][una_eq])
                else:
                    print('j:',resultado[entrada][una_eq]['j'],'; '
                          'i:',resultado[entrada][una_eq]['i'])
                    sym.pprint(resultado[entrada][una_eq]['discreta_ij'])
                    sym.pprint(resultado[entrada][una_eq]['discreta_k'])
    else:
        print('\n',entrada,': \n',resultado[entrada])
print('Resultados para U(x,y)')
print('xi:',xi)
print('yj:',yj)
print(' j, U[i,j]')
for j_k in range(m-1,-1,-1):
    print(j_k, (u_xy[:,j_k]))

EDP Parabólica [ contínua a discreta ] || [ sympy_explícito ] [sympy_ implícito ]

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

Referencia: Sympy ODE solver. Stewart James. Cálculo de varias variables. 17.2 p1148 pdf545. Lathi Ejemplo 2.1.a p155, Ejemplo 2.9 p175, Oppenheim 2.4.1 p117, ejemplo 1 de Modelo entrada-salida,

Continuando con el ejercicio de la sección anterior de una ecuación diferencial ordinaria, lineal y de coeficientes constantes, se obtuvo la solución homogénea o ecuación complementaria con dsolve().  Para el desarrollo analítico de la solución se requieren describir los pasos tales como la ecuación auxiliar, como se describe a continuación como un ejercicio de análisis de expresiones con Sympy.

Ejemplo 1. Ecuación diferencial de un circuito RLC, ecuación complementaria.circuito RLC

El circuito RLC con entrada x(t) de la figura tiene una corriente y(t) o salida del sistema que se representa  por medio de una ecuación diferencial.

\frac{d^{2}y(t)}{dt^{2}} + 3\frac{dy(t)}{dt} + 2y(t) = \frac{dx(t)}{dt}

Las condiciones iniciales del sistema a tiempo t=0 son y(0)=0 , y'(0)=-5

1. Ecuación diferencial y condiciones de frontera o iniciales

La ecuación diferencial del ejercicio con Sympy se escribe con el operador sym.diff(y,t,1). Indicando la variable independiente ‘t‘, que ‘y‘ es una función y el orden de la derivada con 1. esquema que se mantiene para la descripción de las condiciones iniciales.

# ecuacion: lado izquierdo = lado derecho
#           Left Hand Side = Right Hand Side
LHS = sym.diff(y(t),t,2) + 3*sym.diff(y(t),t,1) + 2*y(t)
RHS = sym.diff(x(t),t,1)
ecuacion = sym.Eq(LHS,RHS)

# condiciones de frontera o iniciales
y_cond = {y(0) : 0,
          sym.diff(y(t),t,1).subs(t,0) : -5}

2. Clasificar la ecuación diferencial ordinaria

Sympy permite revisar el tipo de EDO a resolver usando la instrucción:

>>> sym.classify_ode(ecuacion, y(t))
('factorable', 
 'nth_linear_constant_coeff_variation_of_parameters', 
 'nth_linear_constant_coeff_variation_of_parameters_Integral')
>>>

3. Ecuación homogenea

Para el análisis de la respuesta del circuito, se inicia haciendo la entrada x(t)=0, también conocida como ecuación para «respuesta a entrada cero» o ecuación homogenea.

\frac{d^{2}y(t)}{dt^{2}} + 3\frac{dy(t)}{dt} + 2y(t) = 0

En el algoritmo, La ecuación homogénea se obtiene al substituir x(t)=0 en cada lado de la ecuación, que también es un paso para encontrar la respuesta a entrada cero. La ecuacion homogenea se escribe de la forma f(t)=0.

# ecuación homogénea x(t)=0, entrada cero
RHSx0 = ecuacion.rhs.subs(x(t),0).doit()
LHSx0 = ecuacion.lhs.subs(x(t),0).doit()
homogenea = LHSx0 - RHSx0

4. Ecuación auxiliar o característica.

En la ecuación homogenea , se procede a sustituir en la ecuación el operador dy/dt de la derivada con una variable r elevada al orden de la derivada de cada término . El resultado es una ecuación algebraica que se analiza encontrando las raíces de r.

r ^2 + 3r +2 = 0

Los valores de ‘r’ que resuelven la ecuación permiten estimar la forma de la solución para y(t) conocida como la ecuación general

Se usan diferentes formas para mostrar la ecuación auxiliar, pues en algunos casos se requiere usar la forma de factores, en otros casos los valores de las raíces y las veces que se producen. Formas usadas para generar diagramas de polos y ceros, o expresiones de transformada de Laplace. Motivo por el que los resultados se los presenta en un diccionario, y asi usar la respuesta que sea de interés en cada caso.

homogenea :
                        2          
           d           d           
2*y(t) + 3*--(y(t)) + ---(y(t)) = 0
           dt           2          
                      dt           
auxiliar : r**2 + 3*r + 2
Q : r**2 + 3*r + 2
Q_factor : (r + 1)*(r + 2)
Q_raiz : {-1: 1, -2: 1}
>>> 

Instrucciones con Python

# Ecuación Diferencial Ordinaria EDO
# ecuación auxiliar o característica 
# Lathi 2.1.a pdf 155, (D^2+ 3D + 2)y = Dx
import sympy as sym

# INGRESO
t = sym.Symbol('t', real=True)
r = sym.Symbol('r')
y = sym.Function('y')
x = sym.Function('x')

# ecuacion: lado izquierdo = lado derecho
#           Left Hand Side = Right Hand Side
LHS = sym.diff(y(t),t,2) + 3*sym.diff(y(t),t,1) + 2*y(t)
RHS = sym.diff(x(t),t,1)
ecuacion = sym.Eq(LHS,RHS)

# condiciones de frontera o iniciales
y_cond = {y(0) : 0,
          sym.diff(y(t),t,1).subs(t,0) : -5}

# PROCEDIMIENTO
def edo_lineal_auxiliar(ecuacion,
                 t = sym.Symbol('t'),r = sym.Symbol('r'),
                 y = sym.Function('y'),x = sym.Function('x')):
    ''' ecuacion auxiliar o caracteristica de EDO
        t independiente
    '''
    # ecuación homogénea x(t)=0, entrada cero
    RHSx0 = ecuacion.rhs.subs(x(t),0).doit()
    LHSx0 = ecuacion.lhs.subs(x(t),0).doit()
    homogenea = LHSx0 - RHSx0
    homogenea = sym.expand(homogenea,t)

    # ecuación auxiliar o característica
    Q = 0*r
    term_suma = sym.Add.make_args(homogenea)
    for term_k in term_suma:
        orden_k = sym.ode_order(term_k,y)
        coef = 1 # coefientes del término suma
        factor_mul = sym.Mul.make_args(term_k)
        for factor_k in factor_mul:
            cond = factor_k.has(sym.Derivative)
            cond = cond or factor_k.has(y(t))
            if not(cond):
                coef = coef*factor_k
        Q = Q + coef*(r**orden_k)
               
    # Q factores y raices
    Q_factor = sym.factor(Q,r)
    Q_poly   = sym.poly(Q,r)
    Q_raiz   = sym.roots(Q_poly)
    
    auxiliar = {'homogenea' : sym.Eq(homogenea,0),
                'auxiliar'  : Q,
                'Q'         : Q,
                'Q_factor'  : Q_factor,
                'Q_raiz'    : Q_raiz }
    return(auxiliar)

# analiza la ecuación diferencial
edo_tipo = sym.classify_ode(ecuacion, y(t))
auxiliar = edo_lineal_auxiliar(ecuacion,t,r)

# SALIDA
print('clasifica EDO:')
for untipo in edo_tipo:
    print(' ',untipo)
print('ecuacion auxiliar:')
for entrada in auxiliar:
    print(' ',entrada,':',auxiliar[entrada])

Otra forma de mostrar el resultado es usando un procedimiento creado para mostrar elementos del diccionario según corresponda a cada tipo. el procedimiento se adjunta al final.

Las funciones se incorporan a los algoritmos del curso en matg1052.py


6. Ecuación diferencial, ecuación general y complementaria con Sympy-Python

La ecuación general se encuentra con la instrucción sym.dsolve(), sin condiciones iniciales, por que el resultado presenta constantes Ci por determinar.

    # solucion general de ecuación homogénea
    general = sym.dsolve(homogenea, y(t))
    general = general.expand()

el resultado para el ejercicio es

general :
           -t       -2*t
y(t) = C1*e   + C2*e

Para encontrar los valores de las constantes, se aplica cada una de las condiciones iniciales a la ecuación general, obteniendo un sistema de ecuaciones.

     # Aplica condiciones iniciales o de frontera
    eq_condicion = []
    for cond_k in y_cond: # cada condición
        valor_k = y_cond[cond_k]
        orden_k = sym.ode_order(cond_k,y)
        if orden_k==0: # condicion frontera
            t_k = cond_k.args[0] # f(t_k)
            expr_k = general.rhs.subs(t,t_k)
        else: # orden_k>0
            # f.diff(t,orden_k).subs(t,t_k)
            subs_param = cond_k.args[2] # en valores
            t_k = subs_param.args[0]  # primer valor
            dyk = general.rhs.diff(t,orden_k)
            expr_k = dyk.subs(t,t_k)
        eq_condicion.append(sym.Eq(valor_k,expr_k))

con el siguiente resultado:

eq_condicion :
0 = C1 + C2
-5 = -C1 - 2*C2

El sistema de ecuaciones se resuelve con la instrucción

constante = sym.solve(eq_condicion)

que entrega un diccionario con cada valor de la constante

constante : {C1: -5, C2: 5}

La ecuación complementaria se obtiene al sustituir en la ecuación general los valores de las constantes.

El resultado del algoritmo

homogenea :
                        2          
           d           d           
2*y(t) + 3*--(y(t)) + ---(y(t)) = 0
           dt           2          
                      dt           
general :
           -t       -2*t
y(t) = C1*e   + C2*e    
eq_condicion :
0 = C1 + C2
-5 = -C1 - 2*C2
constante : {C1: -5, C2: 5}
complementaria :
            -t      -2*t
y(t) = - 5*e   + 5*e    
>>>

Instrucciones con Python

Las instrucciones detalladas se presentan en el algoritmo integrado.

# Solución complementaria a una Ecuación Diferencial Ordinaria EDO
# Lathi 2.1.a pdf 155, (D^2+ 3D + 2)y = Dx
import sympy as sym
equivalentes = [{'DiracDelta': lambda x: 1*(x==0)},
                {'Heaviside': lambda x,y: np.heaviside(x, 1)},
                'numpy',]
import matg1052 as fcnm

# INGRESO
t = sym.Symbol('t', real=True)
r = sym.Symbol('r')
y = sym.Function('y')
x = sym.Function('x')

# ecuacion: lado izquierdo = lado derecho
#           Left Hand Side = Right Hand Side
LHS = sym.diff(y(t),t,2) + 3*sym.diff(y(t),t,1) + 2*y(t)
RHS = sym.diff(x(t),t,1)
ecuacion = sym.Eq(LHS,RHS)

# condiciones de frontera o iniciales
y_cond = {y(0) : 0,
          sym.diff(y(t),t,1).subs(t,0) : -5}

# PROCEDIMIENTO
def edo_lineal_complemento(ecuacion,y_cond):
    # ecuación homogénea x(t)=0, entrada cero
    RHSx0 = ecuacion.rhs.subs(x(t),0).doit()
    LHSx0 = ecuacion.lhs.subs(x(t),0).doit()
    homogenea = LHSx0 - RHSx0

    # solucion general de ecuación homogénea
    general = sym.dsolve(homogenea, y(t))
    general = general.expand()

    # Aplica condiciones iniciales o de frontera
    eq_condicion = []
    for cond_k in y_cond: # cada condición
        valor_k = y_cond[cond_k]
        orden_k = sym.ode_order(cond_k,y)
        if orden_k==0: # condicion frontera
            t_k    = cond_k.args[0] # f(t_k)
            expr_k = general.rhs.subs(t,t_k)
        else: # orden_k>0
            subs_param = cond_k.args[2] # en valores
            t_k = subs_param.args[0]  # primer valor
            dyk = sym.diff(general.rhs,t,orden_k)
            expr_k = dyk.subs(t,t_k)
        eq_condicion.append(sym.Eq(valor_k,expr_k))

    constante = sym.solve(eq_condicion)

    # ecuacion complementaria
    # reemplaza las constantes en general
    yc = general
    for Ci in constante:
        yc = yc.subs(Ci, constante[Ci])
    
    complemento = {'homogenea'      : sym.Eq(homogenea,0),
                   'general'        : general,
                   'eq_condicion'   : eq_condicion,
                   'constante'      : constante,
                   'complementaria' : yc}
    return(complemento)

# analiza ecuacion diferencial
edo_tipo = sym.classify_ode(ecuacion, y(t))
auxiliar = fcnm.edo_lineal_auxiliar(ecuacion,t,r)
complemento = edo_lineal_complemento(ecuacion,y_cond)
yc = complemento['complementaria'].rhs

# SALIDA
print('clasifica EDO:')
for elemento in edo_tipo:
    print(' ',elemento)

fcnm.print_resultado_dict(auxiliar)
fcnm.print_resultado_dict(complemento)

7. Gráfica de la ecuación complementaria yc

Para la gráfica se procede a convertir la ecuación yc en la forma numérica con sym.lambdify(). Se usa un intervalo de tiempo entre [0,5] y 51 muestras con el siguiente resultado:

solucion homogenea EDO

Instrucciones en Python

# GRAFICA ------------------
import numpy as np
import matplotlib.pyplot as plt

# INGRESO
t_a = 0 ; t_b = 5 # intervalo tiempo [t_a,t_b]
muestras = 51

# PROCEDIMIENTO
yt = sym.lambdify(t,yc, modules=equivalentes) 
ti = np.linspace(t_a,t_b,muestras)
yi = yt(ti)

# Grafica
plt.plot(ti,yi, color='orange', label='yc(t)')
untitulo = 'yc(t)=$'+ str(sym.latex(yc))+'$'
plt.title(untitulo)
plt.xlabel('t')
plt.ylabel('yc(t)')
plt.legend()
plt.grid()
plt.show()

8.  Mostrar el resultado del diccionario según el tipo de entrada

Para mostrar los resultados de una forma más fácil de leer, se usará sym.pprint() para las ecuaciones, y print() simple para los demás elementos del resultado

def print_resultado_dict(resultado):
    ''' print de diccionario resultado
        formato de pantalla
    '''
    for entrada in resultado:
        tipo = type(resultado[entrada])
        if tipo == sym.core.relational.Equality:
            print(entrada,':')
            sym.pprint(resultado[entrada])
        elif tipo==list or tipo==tuple:
            tipoelem = type(resultado[entrada][0])
            if tipoelem==sym.core.relational.Equality:
                print(entrada,':')
                for fila in resultado[entrada]:
                    sym.pprint(fila)
            else:
                print(entrada,':')
                for fila in resultado[entrada]:
                    print(' ',fila)
        else:
            print(entrada,':',resultado[entrada])
    return()

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

Referencia: Sympy ODE solver. Stewart James. Cálculo de varias variables. 17.2 p1148 pdf545. Lathi Ejemplo 2.1.a p155, Ejemplo 2.9 p175, Oppenheim 2.4.1 p117, ejemplo 1 de Modelo entrada-salida,

Los métodos analíticos para encontrar una solución y(t) a una ecuación diferencial ordinaria EDO requieren identificar el tipo de la ecuación para establecer el método de solución. Sympy incorpora librerías que pueden asistir en la solución con muchos de los métodos analíticos conocidos.

Ejemplo 1. Ecuación diferencial de un circuito RLC, ecuación complementaria.circuito RLC

El circuito RLC con entrada x(t) de la figura tiene una corriente y(t) o salida del sistema que se representa  por medio de una ecuación diferencial.

\frac{d^{2}y(t)}{dt^{2}} + 3\frac{dy(t)}{dt} + 2y(t) = \frac{dx(t)}{dt}

Las condiciones iniciales del sistema a tiempo t=0 son y(0)=0 , y'(0)=-5

Encuentre la solución considerando como entrada x(t)

x(t) = 10 e^{-3t} \mu(t)

1.1 Ecuación diferencial y condiciones de frontera o iniciales

La ecuación diferencial del ejercicio con Sympy se escribe con el operador sym.diff(y,t,k). Indicando la variable independiente, que ‘y‘ es una función y el orden de la derivada con k.

La ecuación puede escribirse como lado izquierdo (LHS) es igual al lado derecho (RHS). Una ecuación en Sympy se compone de las dos partes sym.Eq(LHS,RHS).

# INGRESO
t = sym.Symbol('t', real=True)
r = sym.Symbol('r')
y = sym.Function('y')
x = sym.Function('x')

# ecuacion: lado izquierdo = lado derecho
#           Left Hand Side = Right Hand Side
LHS = sym.diff(y(t),t,2) + 3*sym.diff(y(t),t,1) + 2*y(t)
RHS = sym.diff(x(t),t,1)
ecuacion = sym.Eq(LHS,RHS)

Las condiciones de frontera o iniciales y(0)=0, y'(0)=-5 se ingresan en un diccionario en el siguiente formato

# condiciones de frontera o iniciales
y_cond = {y(0) : 0,
          sym.diff(y(t),t,1).subs(t,0) : -5}

la entrada x(t) conocida se define como:

# ecuacion entrada x(t)
xp = 10*sym.exp(-3*t)*sym.Heaviside(t)

1.2 Clasificar la ecuación diferencial ordinaria

Para revisar el tipo de EDO a resolver se tiene la instrucción:

>>> sym.classify_ode(ecuacion, y(t))
('factorable', 
 'nth_linear_constant_coeff_variation_of_parameters', 
 'nth_linear_constant_coeff_variation_of_parameters_Integral')
>>>

Solución EDO como suma de solución complementaria y solución particular

La solución para una EDO puede presentarse también como la suma de una solución complementaria y una solución particular

y(t) = y_p(t) + y_c(t)

1.3 Solución complementaria yc(t)

La solución complementaria de la EDO se interpreta también como respuesta de entrada cero del circuito, donde la salida y(t) depende solo de las condiciones iniciales del circuito. Para el ejemplo, considera solo las energías internas almacenadas en el capacitor y el inductor. Para éste caso x(t) = 0

Es necesario crear la ecuación homogénea con x(t)=0 en la ecuación diferencial para encontrar la solución con las condiciones iniciales dadas en el enunciado del ejercicio.

# ecuación homogénea x(t)=0, entrada cero
RHSx0 = ecuacion.rhs.subs(x(t),0).doit()
LHSx0 = ecuacion.lhs.subs(x(t),0).doit()
homogenea = LHSx0 - RHSx0

# solucion general de ecuación homogénea
yc = sym.dsolve(homogenea, y(t),ics=y_cond)
yc = yc.expand()

el resultado de la ecuación complementaria es

solucion complementaria y_c(t): 
            -t      -2*t
y(t) = - 5*e   + 5*e    

1.4 Solución particular yp(t)

En el caso de la solución particular, se simplifica la ecuación diferencial al sustituir x(t) con la entrada particular xp(t). Las condiciones iniciales en este caso consideran que el circuito no tenía energía almacenada en sus componentes (estado cero), por lo que las condiciones iniciales no se consideran.

# ecuación particular x(t)=0, estado cero
RHSxp = ecuacion.rhs.subs(x(t),x_p).doit()
LHSxp = ecuacion.lhs.subs(x(t),x_p).doit()
particular = LHSxp - RHSxp

# solucion particular de ecuación homogénea
yp = sym.dsolve(particular, y(t))

Se aplica nuevamente la búsqueda de la solución con sym.dsolve() y se obtiene la solución particular que incluye parte del modelo de la ecuación homogénea con los coeficientes Ci

>>> sym.pprint(yp.expand())
           -t       -2*t      -t                    -2*t                    -3*t
y(t) = C1*e   + C2*e     - 5*e  *Heaviside(t) + 20*e    *Heaviside(t) - 15*e    *Heaviside(t)

>>> yp.free_symbols
{C2, C1, t}
>>>

En éste punto se incrementan los términos de las constantes C1 y C2 como símbolos, que para tratar exclusivamente la solución particular, se reemplazan con ceros. Para obtener las variables de la ecuación se usa la instrucción yp.free_symbols que entrega un conjunto. El conjunto se descarta el símbolo t y se sustituye todos los demás por cero en la ecuación.

    # particular sin terminos Ci
    y_Ci = yp.free_symbols
    y_Ci.remove(t) # solo Ci
    for Ci in y_Ci: 
        yp = yp.subs(Ci,0)

quedando la solución particular como:

            -t                    -2*t                    -3*t             
y(t) = - 5*e  *Heaviside(t) + 20*e    *Heaviside(t) - 15*e    *Heaviside(t)

1.5 Solución EDO como suma de complementaria + particular

La solución de la ecuación diferencial ordinaria, ante la entrada x(t) y condiciones iniciales es la suma de los componentes complementaria y particular:

# solucion total = complementaria + particular
ytotal = yc.rhs + yp.rhs

El resultado de todo el algoritmo  permite interpretar la respuesta con los conceptos de las respuestas del circuito a entrada cero y estado cero, que facilitan el análisis de las soluciones en ejercicios de aplicación.

ecuación diferencial a resolver: 
                        2                 
           d           d          d       
2*y(t) + 3*--(y(t)) + ---(y(t)) = --(x(t))
           dt           2         dt      
                      dt                  
condiciones iniciales
  y(0)  =  0
  Subs(Derivative(y(t), t), t, 0)  =  -5
clasifica EDO:
  factorable
  nth_linear_constant_coeff_variation_of_parameters
  nth_linear_constant_coeff_variation_of_parameters_Integral
x(t): 
    -3*t             
10*e    *Heaviside(t)
solucion complementaria yc(t): 
            -t      -2*t
y(t) = - 5*e   + 5*e    
solucion particular yp(t): 
       /     -t       -2*t       -3*t\             
y(t) = \- 5*e   + 20*e     - 15*e    /*Heaviside(t)
solucion y(t) = yc(t) +yp(t): 
/     -t       -2*t       -3*t\                   -t      -2*t
\- 5*e   + 20*e     - 15*e    /*Heaviside(t) - 5*e   + 5*e    
>>> 

Instrucciones en Python

# Solución de unaEcuación Diferencial Ordinaria EDO
# Lathi 2.1.a pdf 155, (D^2+ 3D + 2)y = Dx
import sympy as sym

# INGRESO
t = sym.Symbol('t', real=True)
y = sym.Function('y')
x = sym.Function('x')

# ecuacion: lado izquierdo = lado derecho
#           Left Hand Side = Right Hand Side
LHS = sym.diff(y(t),t,2) + 3*sym.diff(y(t),t,1) + 2*y(t)
RHS = sym.diff(x(t),t,1)
ecuacion = sym.Eq(LHS,RHS)

# condiciones de frontera o iniciales
y_cond = {y(0) : 0,
          sym.diff(y(t),t,1).subs(t,0) : -5}

# ecuacion entrada x(t)
xp = 10*sym.exp(-3*t)*sym.Heaviside(t)

# PROCEDIMIENTO
edo_tipo = sym.classify_ode(ecuacion, y(t))

# ecuación homogénea x(t)=0, entrada cero
RHSx0 = ecuacion.rhs.subs(x(t),0).doit()
LHSx0 = ecuacion.lhs.subs(x(t),0).doit()
homogenea = LHSx0 - RHSx0

# solucion general de ecuación homogénea
yc = sym.dsolve(homogenea, y(t),ics=y_cond)
yc = yc.expand()

def edo_lineal_particular(ecuacion,xp):
    ''' edo solucion particular con entrada x(t)
    '''
    # ecuación particular x(t)=0, estado cero
    RHSxp = ecuacion.rhs.subs(x(t),xp).doit()
    LHSxp = ecuacion.lhs.subs(x(t),xp).doit()
    particular = LHSxp - RHSxp

    # solucion particular de ecuación homogénea
    yp = sym.dsolve(particular, y(t))

    # particular sin terminos Ci
    y_Ci = yp.free_symbols
    y_Ci.remove(t) # solo Ci
    for Ci in y_Ci: 
        yp = yp.subs(Ci,0)
        
    # simplifica y(t) y agrupa por escalon unitario
    yp = sym.expand(yp.rhs,t)
    lista_escalon = yp.atoms(sym.Heaviside)
    yp = sym.collect(yp,lista_escalon)
    yp = sym.Eq(y(t),yp)
    
    return(yp)

yp = edo_lineal_particular(ecuacion,xp)
# solucion total
ytotal = yp.rhs + yc.rhs


# SALIDA
print('ecuación diferencial a resolver: ')
sym.pprint(ecuacion)

# condiciones iniciales
print('condiciones iniciales')
for caso in y_cond:
    print(' ',caso,' = ', y_cond[caso])

print('clasifica EDO:')
for elemento in edo_tipo:
    print(' ',elemento)

print('x(t): ')
sym.pprint(xp)

print('solucion complementaria yc(t): ')
sym.pprint(yc)

print('solucion particular yp(t): ')
sym.pprint(yp)

print('solucion y(t) = yc(t) +yp(t): ')
sym.pprint(ytotal)

1.6  Grafica de resultados

Se adjunta la gráfica de los resultados de las ecuaciones complementaria, particular y total

EDO lineal Complementaria Particular 01

Instrucciones adicionales en Python

# GRAFICA ------------------
import numpy as np
import matplotlib.pyplot as plt
equivalentes = [{'DiracDelta': lambda x: 1*(x==0)},
                {'Heaviside': lambda x,y: np.heaviside(x, 1)},
                'numpy',]

# INGRESO
t_a = 0 ; t_b = 5 # intervalo tiempo [t_a,t_b]
muestras = 51

# PROCEDIMIENTO
y_c = sym.lambdify(t,yc.rhs, modules=equivalentes)
y_p = sym.lambdify(t,yp.rhs, modules=equivalentes)
y_cp = sym.lambdify(t,ytotal, modules=equivalentes)
ti = np.linspace(t_a,t_b,muestras)
yci = y_c(ti)
ypi = y_p(ti)
ycpi = y_cp(ti)

# Grafica
plt.plot(ti,yci, color='orange', label='yc(t)')
plt.plot(ti,ypi, color='dodgerblue', label='yp(t)')
plt.plot(ti,ycpi, color='green', label='y(t)')
untitulo = 'y(t)=$'+ str(sym.latex(ytotal))+'$'
plt.title(untitulo)
plt.xlabel('t')
plt.legend()
plt.grid()
plt.show()

Un desarrollo semejante del ejercicio aplicando conceptos del curso de señales y sistemas de proporciona en:

LTI CT Respuesta del Sistema Y(s)=ZIR+ZSR con Sympy-Python

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

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 con Python

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