7.1.4 EDP Parabólicas método implícito – analítico con Sympy-Python

Referencia:  Chapra 30.3 p893 pdf917, Burden 9Ed 12.2 p729, Rodriguez 10.2.4 p415

Para el ejercicio desarrollado en forma analítica como EDP Parabólicas método implícito, se desarrollan los pasos esta vez  usando Sympy

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

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

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.

 EDP=0  :
                    2             
    d              d              
- K*--(u(x, t)) + ---(u(x, t)) = 0
    dt              2             
                  dx              

 Lambda L  :
  Dt 
-----
  2  
Dx *K

 discreta = 0  :
-2*U(i, j) + U(i - 1, j) + U(i + 1, j)   K*(U(i, j) - U(i, j - 1))
-------------------------------------- - -------------------------
                   2                                 Dt           
                 Dx                                               

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

Ecuaciones de iteraciones j:
 j = 1
-1.5*U(1, 1) + 0.25*U(2, 1) = -40.0
0.25*U(1, 1) - 1.5*U(2, 1) + 0.25*U(3, 1) = -25.0
0.25*U(2, 1) - 1.5*U(3, 1) + 0.25*U(4, 1) = -25.0
0.25*U(3, 1) - 1.5*U(4, 1) + 0.25*U(5, 1) = -25.0
0.25*U(4, 1) - 1.5*U(5, 1) + 0.25*U(6, 1) = -25.0
0.25*U(5, 1) - 1.5*U(6, 1) + 0.25*U(7, 1) = -25.0
0.25*U(6, 1) - 1.5*U(7, 1) + 0.25*U(8, 1) = -25.0
0.25*U(7, 1) - 1.5*U(8, 1) + 0.25*U(9, 1) = -25.0
0.25*U(8, 1) - 1.5*U(9, 1) = -35.0
 j = 2
-1.5*U(1, 2) + 0.25*U(2, 2) = -46.0050525095364
0.25*U(1, 2) - 1.5*U(2, 2) + 0.25*U(3, 2) = -26.0303150572187
0.25*U(2, 2) - 1.5*U(3, 2) + 0.25*U(4, 2) = -25.1768378337756
0.25*U(3, 2) - 1.5*U(4, 2) + 0.25*U(5, 2) = -25.030711945435
0.25*U(4, 2) - 1.5*U(5, 2) + 0.25*U(6, 2) = -25.0074338388344
0.25*U(5, 2) - 1.5*U(6, 2) + 0.25*U(7, 2) = -25.0138910875712
0.25*U(6, 2) - 1.5*U(7, 2) + 0.25*U(8, 2) = -25.0759126865931
0.25*U(7, 2) - 1.5*U(8, 2) + 0.25*U(9, 2) = -25.4415850319874
0.25*U(8, 2) - 1.5*U(9, 2) = -37.5735975053312
 j = 3
-1.5*U(1, 3) + 0.25*U(2, 3) = -50.2512763902537
0.25*U(1, 3) - 1.5*U(2, 3) + 0.25*U(3, 3) = -27.4874483033763
0.25*U(2, 3) - 1.5*U(3, 3) + 0.25*U(4, 3) = -25.5521532011296
0.25*U(3, 3) - 1.5*U(4, 3) + 0.25*U(5, 3) = -25.1181195682987
0.25*U(4, 3) - 1.5*U(5, 3) + 0.25*U(6, 3) = -25.0337164269226
0.25*U(5, 3) - 1.5*U(6, 3) + 0.25*U(7, 3) = -25.0544436378994
0.25*U(6, 3) - 1.5*U(7, 3) + 0.25*U(8, 3) = -25.2373810501889
0.25*U(7, 3) - 1.5*U(8, 3) + 0.25*U(9, 3) = -26.0661919168616
0.25*U(8, 3) - 1.5*U(9, 3) = -39.3934303230311
 j = 4
-1.5*U(1, 4) + 0.25*U(2, 4) = -53.3449104170748
0.25*U(1, 4) - 1.5*U(2, 4) + 0.25*U(3, 4) = -29.0643569414341
0.25*U(2, 4) - 1.5*U(3, 4) + 0.25*U(4, 4) = -26.0914380180247
0.25*U(3, 4) - 1.5*U(4, 4) + 0.25*U(5, 4) = -25.2756583621956
0.25*U(4, 4) - 1.5*U(5, 4) + 0.25*U(6, 4) = -25.0900338819544
0.25*U(5, 4) - 1.5*U(6, 4) + 0.25*U(7, 4) = -25.1296792218401
0.25*U(6, 4) - 1.5*U(7, 4) + 0.25*U(8, 4) = -25.4702668974885
0.25*U(7, 4) - 1.5*U(8, 4) + 0.25*U(9, 4) = -26.7423979623353
0.25*U(8, 4) - 1.5*U(9, 4) = -40.7193532090766

Solución u:
[[60.         60.         60.         60.         60.        ]
 [25.         31.00505251 35.25127639 38.34491042 40.6652135 ]
 [25.         26.03031506 27.4874483  29.06435694 30.61163935]
 [25.         25.17683783 25.5521532  26.09143802 26.74719483]
 [25.         25.03071195 25.11811957 25.27565836 25.50577757]
 [25.         25.00743384 25.03371643 25.09003388 25.18483712]
 [25.         25.01389109 25.05444364 25.12967922 25.24310963]
 [25.         25.07591269 25.23738105 25.4702669  25.75510376]
 [25.         25.44158503 26.06619192 26.74239796 27.40644533]
 [25.         27.57359751 29.39343032 30.71935321 31.71397636]
 [40.         40.         40.         40.         40.        ]]
>>> 

Instrucciones en Python

El desarrollo de la ecuación en forma discreta se realizó en EDP Parabólica método explícito – analítico con Sympy-Python, por lo que se incorpora a éste algoritmo como una función, desde donde se inicia los parsos para el método implicito.

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
# método implícito con Sympy
import numpy as np
import sympy as sym

# variables continuas
t = sym.Symbol('t',real=True)
x = sym.Symbol('x',real=True)
u = sym.Function('u')(x,t)
K = sym.Symbol('K',real=True)
# variables discretas
i  = sym.Symbol('i',integer=True,positive=True)
j  = sym.Symbol('j',integer=True,positive=True)
Dx = sym.Symbol('Dx',real=True,positive=True)
Dt = sym.Symbol('Dt',real=True,positive=True)
L  = sym.Symbol('L',real=True)
U  = sym.Function('U')(i,j)

# INGRESO
# ecuacion: LHS=RHS
LHS = sym.diff(u,x,2)
RHS = K*sym.diff(u,t,1)
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,t,1): (U - U.subs(j,j-1))/Dt}

# parametros
Ta = 60 ; Tb = 40 ; T0 = 25
K_ = 4
x_a = 0 ; x_b = 1
x_tramos = 10
t_a = 0 # t_b depende de t_tramos y dt
t_tramos = 4

dx  = (x_b-x_a)/x_tramos
dt  = dx/10
L_k = dt/(dx**2)/K_

# PROCEDIMIENTO
def edp_discreta(EDP,dif_dividida):
    ''' EDP contínua y diferencias divididas a usar
    '''
    desarrollo={}
    # expresión a la izquierda LHS (Left Hand side)
    if EDP.rhs!=0:
        LHS = EDP.lhs
        RHS = EDP.rhs
        EDP = sym.Eq(LHS-RHS,0)
    desarrollo['EDP=0'] = EDP

    lamb = Dt/K/(Dx**2)
    desarrollo['Lambda L'] = lamb

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

    # simplifica con lambda L
    discreta_L = sym.expand(discreta*Dt/K)
    discreta_L = discreta_L.subs(Dt/K/(Dx**2),L)
    discreta_L = sym.collect(discreta_L,U)
    discreta_L = sym.Eq(discreta_L,0)

    desarrollo['discreta_L = 0'] = discreta_L
    return(desarrollo)

edp_discreta = edp_discreta(EDP,dif_dividida)

# u[x,t] resultados, valores iniciales
xi = np.arange(x_a,x_b+dx,dx)
ti = np.arange(t_a,t_a+t_tramos*dt+dt,dt)

uxt = np.zeros(shape=(x_tramos+1,t_tramos+1),
               dtype=float)
uxt = uxt*np.nan
uxt[0,:] = Ta
uxt[1:x_tramos+1,0] = T0
uxt[x_tramos,:] = Tb


# EDP parabólicas método implícito 
# conoce valor de uxt[i,j] 
conoce = np.zeros(shape=(x_tramos+1,t_tramos+1),
                  dtype=bool)
conoce[:,0]  = 1 # inferior
conoce[0,:]  = 1 # izquierda
conoce[-1,:] = 1 # derecha

# ecuaciones por iteracion j
eq_itera = []
for jk in range(1,t_tramos+1,1):
    eq_jk = []
    for ik in range(1,x_tramos,1):
        LHS = edp_discreta['discreta_L = 0'].lhs
        RHS = 0*t
        LHS = LHS.subs([(i,ik),(j,jk),(L,L_k)])
        # revisa terminos de suma
        term_suma = sym.Add.make_args(LHS)
        for term_k in term_suma:
            cambiar = False
            cambiar_valor = 0 ; cambiar_factor = 0
            # revisa factores en termino de suma
            factor_Mul = sym.Mul.make_args(term_k)
            for factor_k in factor_Mul:
                if factor_k.is_Function: # es U(i,j)
                    [i_k,j_k] = factor_k.args
                    cambiar = conoce[i_k,j_k]
                    cambiar_factor = factor_k
                    cambiar_valor = uxt[i_k,j_k]
            # U(i,j)conocidos pasan a la derecha
            if cambiar: 
                RHS = RHS-term_k
                RHS = RHS.subs(cambiar_factor,
                               cambiar_valor)
                LHS = LHS-term_k
        # agrupa ecuaciones de iteración i
        eq_jk.append(sym.Eq(LHS,RHS))

    # resuelve sistema ecuaciones A.X=B
    X_jk = sym.solve(eq_jk)[0]
    
    # actualiza uxt(i,j) , conoce[i,j]
    for nodoij in X_jk:
        [i_k,j_k] = nodoij.args
        uxt[i_k,j_k] = X_jk[nodoij]
        conoce[i_k,j_k] = True

    # agrupa ecuaciones por cada j
    eq_itera.append(eq_jk)
        
# SALIDA
for k in edp_discreta:
    print('\n',k,' :')
    sym.pprint(edp_discreta[k])
print('\nEcuaciones de iteraciones j:')
for jk in range(0,t_tramos,1):
    print(' j =',jk+1)
    for una_eq in eq_itera[jk]:
        sym.pprint(una_eq)
print('\nSolución u:')
print(uxt)

 

7.1.2 EDP Parabólica método explícito – analítico con Sympy-Python

Desarrollo analítico del método explícito de una ecuación diferencial parcial con Sympy

Ejemplos:

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

El algoritmo desarrolla la parte analítica de una EDP modelo con el siguiente resultado:

EDP=0  :
                    2             
    d              d              
- K*--(u(x, t)) + ---(u(x, t)) = 0
    dt              2             
                  dx              

 Lambda L  :
  dt 
-----
    2
K*dx 

 discreta = 0  :
  K*(-U(i, j) + U(i, j + 1))   -2*U(i, j) + U(i - 1, j) + U(i + 1, j)
- -------------------------- + --------------------------------------
              dt                                  2                  
                                                dx                   

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

 U(i, j + 1)  :
L*U(i - 1, j) + L*U(i + 1, j) + (1 - 2*L)*U(i, j)

Instrucciones en Phyton

# Ecuaciones Diferenciales Parciales Parabólicas con Sympy
import sympy as sym

# variables continuas
t = sym.Symbol('t',real=True)
x = sym.Symbol('x',real=True)
u = sym.Function('u')(x,t)
K = sym.Symbol('K',real=True)
# variables discretas
i = sym.Symbol('i',integer=True,positive=True)
j = sym.Symbol('j',integer=True,positive=True)
dx = sym.Symbol('dx',real=True,positive=True)
dt = sym.Symbol('dt',real=True,positive=True)
L = sym.Symbol('L',real=True)
U = sym.Function('U')(i,j)

# INGRESO
# ecuacion: LHS=RHS
LHS = sym.diff(u,x,2)
RHS = K*sym.diff(u,t,1)
ecuacion = 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,t,1): (U.subs(j,j+1)-U)/dt}

buscar = U.subs(j,j+1)

# PROCEDIMIENTO
desarrollo={}
desarrollo['EDP=0'] = ecuacion

lamb = dt/K/(dx**2)
desarrollo['Lambda L'] = lamb

discreta = ecuacion.lhs # forma discreta EDP
for factor_k in dif_dividida:
    discreta = discreta.subs(factor_k,dif_dividida[factor_k])
desarrollo['discreta = 0'] = discreta

# simplifica con Lambda L
discreta_L = sym.expand(discreta*dt/K)
discreta_L = discreta_L.subs(dt/K/(dx**2),L)
discreta_L = sym.collect(discreta_L,U)
discreta_L = sym.Eq(discreta_L,0)

desarrollo['discreta_L = 0'] = discreta_L

# busca el término no conocido
u_explicito = sym.solve(discreta_L,buscar)
u_explicito = sym.collect(u_explicito[0],U)

desarrollo[buscar] = u_explicito
        
# SALIDA

for k in desarrollo:
    print('\n',k,' :')
    sym.pprint(desarrollo[k])

6.4.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.4 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 ordinara EDO requieren identificar el tipo de la ecuación para establecer el método de solución. Sympy incorpora librerias 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 se x(t) = 0

Es necesario crear la ecuación homogenea 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 homogenea 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

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

2021OctubreEstMetorologica.csv


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

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

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

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

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

Desarrollo

literal a

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

literal b

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

literal c

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

Resultados obtenidos:

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

Instrucciones en Python

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Tarea

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

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

8.3 Regresión polinomial de grado m

Referencia: Chapra 17.2 p482, Burden 8.1 p501

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

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

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

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

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

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

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

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

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

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

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

El error estándar se obtiene como:

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

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

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

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

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

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

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

Ejercicio Chapra 17.5 p484

Los datos de ejemplo para la referencia son:

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

resultado de algoritmo:

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

Instrucciones en Python

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

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

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

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

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

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

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

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

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

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

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

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

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

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

8.2 Regresión por Mínimos cuadrados con Python

Referencia: Chapra 17.1.2 p 469. Burden 8.1 p499

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

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

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

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

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

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

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

simplificando,

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

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

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

cuya solución es:

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

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

Coeficiente de correlación

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

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

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

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

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

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


Algoritmo en Python

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

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

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

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

con las instrucciones:

# mínimos cuadrados, regresión con polinomio grado 1
import numpy as np
import sympy as sym
import matplotlib.pyplot as plt

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

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

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

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

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

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

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

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

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

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

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

Coeficiente de correlación con Numpy

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

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

8.1 Regresión vs interpolación

Referencia: Chapra 17.1 p 466. Burden 8.1 p498

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

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

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

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

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


Descripción con la media yi

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

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

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

Instrucciones Python

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

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

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

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

# SALIDA
print('ymedia = ', ym)

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

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

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

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

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

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

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

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

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

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

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

plt.show()

7.3 EDP hiperbólicas

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Convergencia:

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

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

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

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

j=0

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

Aplicando solo cuando j = 0

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

Algoritmo en Python:

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

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

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

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

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

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

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

# GRAFICA
import matplotlib.pyplot as plt

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

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

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

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

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

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

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

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

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

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

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

7.2.2 EDP Elípticas método implícito

con el resultado desarrollado en EDP elípticas para:

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

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

se puede plantear que:

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

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

j=1, i =1

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

j=1, i =2

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

j=1, i=3

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

j=2, i=1

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

j = 2, i = 2

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

j = 2, i = 3

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

j=3, i = 1

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

j=3, i = 2

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

j=3, i=3

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

con las ecuaciones se arma una matriz:

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

que al resolver el sistema de ecuaciones se obtiene:

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

ingresando los resultados a la matriz u:

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

Algoritmo usado para resolver el problema:

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

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

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


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

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

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

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

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

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

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

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

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

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