Categoría: Sol_3Eva 2011-2020

  • s3Eva_IT2018_T2 Drenaje de estanque

    Ejercicio: 3Eva_IT2018_T2 Drenaje de estanque

    literal a

    Se usa interpolación para encontrar los polinomios que pasan por los puntos seleccionados.

    El error de A(5) se obtiene como la diferencia entre el valor de la tabla y el polinomio del tramo [4,6] evaluado en el punto.

    ordenado:  [6 5 4 3 2 1 0]
    hi:  [0 1 2 3 4 5 6]
    Ai:  [ 0.02  0.18  0.32  0.45  0.67  0.97  1.17]
    
    puntos seleccionados:
    h1:  [0, 2, 4, 6]
    A1:  [ 0.02  0.32  0.67  1.17]
    
    Polinomios por tramos: 
     x = [0,2]
    0.000416666666666669*x**3 + 0.148333333333333*x + 0.02
     x = [2,4]
    0.00416666666666666*x**3 - 0.0224999999999999*x**2 + 0.193333333333333*x - 0.00999999999999984
     x = [4,6]
    -0.00458333333333333*x**3 + 0.0824999999999999*x**2 - 0.226666666666666*x + 0.549999999999999
    
    error en px(5):  0.0637499999999998
    

    se observa que la evaluación se realiza para el polinomio entre [4,6]

    drenaje estanque 02 Ah

    Desarrollo en Python

    # 3ra Evaluación I Término 2018
    # Tema 2. Drenaje de Estanque
    
    import numpy as np
    import matplotlib.pyplot as plt
    import sympy as sym
    
    def traza3natural(xi,yi):
        # Trazador cúbico natural, splines
        # resultado: polinomio en forma simbólica
        n = len(xi)
        # Valores h
        h = np.zeros(n-1, dtype = float)
        for j in range(0,n-1,1):
            h[j] = xi[j+1] - xi[j]
        
        # Sistema de ecuaciones
        A = np.zeros(shape=(n-2,n-2), dtype = float)
        B = np.zeros(n-2, dtype = float)
        S = np.zeros(n, dtype = float)
        A[0,0] = 2*(h[0]+h[1])
        A[0,1] = h[1]
        B[0] = 6*((yi[2]-yi[1])/h[1] - (yi[1]-yi[0])/h[0])
        for i in range(1,n-3,1):
            A[i,i-1] = h[i]
            A[i,i] = 2*(h[i]+h[i+1])
            A[i,i+1] = h[i+1]
            B[i] = 6*((yi[i+2]-yi[i+1])/h[i+1] - (yi[i+1]-yi[i])/h[i])
        A[n-3,n-4] = h[n-3]
        A[n-3,n-3] = 2*(h[n-3]+h[n-2])
        B[n-3] = 6*((yi[n-1]-yi[n-2])/h[n-2] - (yi[n-2]-yi[n-3])/h[n-3])
        
        # Resolver sistema de ecuaciones
        r = np.linalg.solve(A,B)
        # S
        for j in range(1,n-1,1):
            S[j] = r[j-1]
        S[0] = 0
        S[n-1] = 0
        
        # Coeficientes
        a = np.zeros(n-1, dtype = float)
        b = np.zeros(n-1, dtype = float)
        c = np.zeros(n-1, dtype = float)
        d = np.zeros(n-1, dtype = float)
        for j in range(0,n-1,1):
            a[j] = (S[j+1]-S[j])/(6*h[j])
            b[j] = S[j]/2
            c[j] = (yi[j+1]-yi[j])/h[j] - (2*h[j]*S[j]+h[j]*S[j+1])/6
            d[j] = yi[j]
        
        # Polinomio trazador
        x = sym.Symbol('x')
        polinomio = []
        for j in range(0,n-1,1):
            ptramo = a[j]*(x-xi[j])**3 + b[j]*(x-xi[j])**2 + c[j]*(x-xi[j])+ d[j]
            ptramo = ptramo.expand()
            polinomio.append(ptramo)
        
        return(polinomio)
    
    # PROGRAMA -------------------------
    
    hi = np.array([6, 5, 4, 3, 2, 1, 0])
    Ai = np.array([1.17, 0.97, 0.67, 0.45, 0.32, 0.18, 0.02])
    xk = 5
    
    # PROCEDIMIENTO LITERAL a
    # reordena en forma ascendente
    ordenado = np.argsort(hi)
    hi = hi[ordenado]
    Ai = Ai[ordenado]
    
    # Selecciona puntos
    xi = [0,2,4,6]
    fi = Ai[xi]
    n = len(xi)
    
    polinomio = traza3natural(xi,fi)
    
    # literal a, estima error
    px = polinomio[2]
    pxk = px.subs('x',xk)
    errado = np.abs(Ai[xk] - pxk)
    
    # SALIDA
    print('ordenado: ', ordenado)
    print('hi: ', hi)
    print('Ai: ', Ai)
    print('puntos seleccionados:')
    print('h1: ', xi)
    print('A1: ', fi)
    
    print('Polinomios por tramos: ')
    for tramo in range(1,n,1):
        print(' x = ['+str(xi[tramo-1])+','+str(xi[tramo])+']')
        print(str(polinomio[tramo-1]))
    
    print('error en px(5): ', errado)
    
    # GRAFICA
    # Puntos para grafica en cada tramo
    resolucion = 10 # entre cada par de puntos
    xtrazado = np.array([])
    ytrazado = np.array([])
    tramo = 1
    while not(tramo>=n):
        a = xi[tramo-1]
        b = xi[tramo]
        xtramo = np.linspace(a,b,resolucion)
        
        ptramo = polinomio[tramo-1]
        pxtramo = sym.lambdify('x',ptramo)
        ytramo = pxtramo(xtramo)
        
        xtrazado = np.concatenate((xtrazado,xtramo))
        ytrazado = np.concatenate((ytrazado,ytramo))
        tramo = tramo + 1
    
    # GRAFICA
    # puntos originales
    plt.plot(hi,Ai,'o',label = 'Ai')
    # Trazador cúbico
    plt.plot(xtrazado,ytrazado, label = 'p(h)')
    plt.plot(xi,fi,'o', label = 'Apx')
    plt.title('Trazador cúbico natural (splines)')
    plt.xlabel('x')
    plt.ylabel('y')
    plt.legend()
    plt.grid()
    plt.show()
    

    Literal b

    TAREA

  • s3Eva_IT2018_T1 Intersección de dos círculos

    Ejercicio: 3Eva_IT2018_T1 Intersección de dos círculos

    Para la solución se presentan dos secciones:

    1. Solución particular de intersección de círculos

    2. Solución General de intersección de círculos

    _


    1. Solución Particular de intersección de círculos

    La solución particular se enfoca en el enunciado del ejercicio presentado

    Literal a

    Se grafica las funciones usando Python, para encontrar el rango de búsqueda de raíces.

    intersecta circulos 01

    De la gráfica se usa el 'zoom' y se puede aproximar los valores para la intersección de las curvas estimando raíces en x=1.80 y x=3.56

    Desarrollo numérico

    Se usan las ecuaciones para encontrar la diferencia entre las funciones.

    (x-4)^2 + (y-4)^2 = 5 x^2 + y^2 = 16

    Se despeja la variable y para la primera ecuación:

    (y-4)^2 = 5 - (x-4)^2 y-4 = \sqrt{5 - (x-4)^2} f(x) = y = \sqrt{5 - (x-4)^2} + 4

    la segunda ecuación se transforma en

    x^2 + y^2 = 16 y^2 = 16 - x^2 g(x) = y = \sqrt{16 - x^2}

    La intersección se obtiene restando las ecuaciones, para f(x) se usa la parte inferior del circulo y para g(x) la parte superior de circulo.

    Para buscar las raíces se analiza en el rango de existencia entre las dos funciones:

    [-4,4]\text{ y } [4 -\sqrt{5} ,4 + \sqrt{5}] [-4,4] \text{ y } [1.7639 , 6.2360]

    por lo que la diferencia existe en el rango:

    [1.7639 ,4] \text{diferencia}(x) = f(x)-g(x)

    que es el que se usa para el literal b

    intersecta circulos 02 diferencia


    Literal b

    Las ecuaciones para la diferencia entre las funciones son :

    f_{2} (x) = -\sqrt{5-(x-4)^2}+4 g_{1} (x) = \sqrt{16-x^2}

    Para el método de Newton-Raphson se requieren las derivadas:

    \frac{d f_2}{dx} = \frac{x-4}{ \sqrt{5-(x-4)^2} } \frac{d g_{1}}{dx} = \frac{-x}{ \sqrt{16-x^2} }

    por lo que:

    \frac{d \text{diferencia}}{dx} = \frac{d f_{2}}{dx} - \frac{d g_{1}}{dx}

    Usando el algoritmo con Python se obtienen las raices:

     usando Newton-Raphson
    raices en:  1.80582463574 3.56917099898

    Desarrollo en Python:

    El desarrollo se realiza por partes, en el mismo orden del planteamiento de  los literales

    # 3ra Evaluación I Término 2018
    # Tema 1. Intersección de círculos
    import numpy as np
    import matplotlib.pyplot as plt
    
    # literal a
    
    fx1 = lambda x: np.sqrt(5-(x-4)**2)+4
    fx2 = lambda x: -np.sqrt(5-(x-4)**2)+4
    gx1 = lambda x: np.sqrt(16-x**2)
    gx2 = lambda x: -np.sqrt(16-x**2)
    
    # Rango inicial de análisis (visual)
    a = -5; b = 7
    muestras = 501
    
    # PROCEDIMIENTO
    # Evalua los puntos en el rango
    xi = np.linspace(a,b,muestras)
    fx1i = fx1(xi)
    fx2i = fx2(xi)
    gx1i = gx1(xi)
    gx2i = gx2(xi)
    
    # SALIDA - Gráfica
    plt.plot(xi,fx1i)
    plt.plot(xi,fx2i)
    plt.plot(xi,gx1i)
    plt.plot(xi,gx2i)
    plt.xlabel('x')
    plt.ylabel('y')
    plt.title('Intersección de círculos')
    plt.grid()
    plt.show()
    
    # GRAFICAR las diferencias
    a = 4 - np.sqrt(5)
    b = 4 + np.sqrt(5)
    # PROCEDIMIENTO
    xi = np.linspace(a,b,muestras)
    diferencia = fx2(xi) - gx1(xi)
    # GRAFICA
    plt.plot(xi,diferencia)
    plt.axhline(0)
    plt.xlabel('x')
    plt.ylabel('y')
    plt.title('diferencia entre círculos')
    plt.grid()
    plt.show()
    
    # literal b -----------------------
    def newton_raphson(funcionx, fxderiva, xi, tolera):
        # funciónx y fxderiva en forma numérica
        # xi es el punto inicial de búsqueda
        tramo = abs(2*tolera)
        while (tramo>=tolera):
            xnuevo = xi - funcionx(xi)/fxderiva(xi)
            tramo = abs(xnuevo-xi)
            xi = xnuevo
        return(xi)
    
    funcionx = lambda x: fx2(x) - gx1(x)
    fxderiva = lambda x: (x-4)/np.sqrt(5-(x-4)**2)+x/np.sqrt(16-x**2)
    
    tolera = 0.001
    xi1 = a + tolera
    xi2 = 3.5
    
    raiz1 = newton_raphson(funcionx, fxderiva, xi1, tolera)
    raiz2 = newton_raphson(funcionx, fxderiva, xi2, tolera)
    
    # SALIDA
    print('\n usando Newton-Raphson')
    print('raices en: ', raiz1,raiz2)
    

    _


    2. Solución General de intersección de círculos

    Una solución más general de la intersección de círculos, considerada como para una actividad de mayor duración, revisa previamente si existe un cruce de áreas entre los dos círculos y estima el intervalo donde se encuentran las raíces [xa,xb].

    De existir esa posibilidad, con el intervalo anterior  [xa,xb] busca por un método de búsqueda de raíces las coordenadas de la intersección de las circunferencias.

    2.1 Buscar cruce de áreas entre dos círculos

    El cruce de áreas entre dos círculos se determina comparando si la distancia entre la suma de los radios es mayor o igual a la distancia entre los centros de los círculos.

    De cumplirse la condición anterior, es posible encontrar las intersecciones de los círculos. El valor xa se obtiene como el mayor entre los límites x hacia la izquierda de cada círculo, mientras que xb se obtiene como el límite x hacia la derecha entre los círculos.

    intersecta Circulos 00

    Lo siguiente que hay que reconocer es cuál de las partes (superior e inferior) de cada círculo es necesario usar para encontrar las intersecciones. Esta sección es necesaria puesto que la fórmula que describe el círculo contiene una raiz cuadrada que puede se positiva o negativa, generando dos segmentos en cada círculo.

    Por ejemplo, partiendo de la fórmula general de un círculo con centro en [x1,y1] y radio r1:

    (x-x_1)^2 + (y-y_1)^2 = r_1^2 (y-y_1)^2 = r_1^2 - (x-x_1)^2 \sqrt{(y-y_1)^2} = \sqrt{r_1^2 - (x-x_1)^2} y = \sqrt{r_1^2 - (x-x_1)^2} + y_1

    Con lo que se muestra la necesidad de identificar para cada círculo el sector arriba y abajo que interviene para encontrar las intersecciones. El orden del sector se establece con las posibles combinaciones de:

    tabla de signos en raíz cuadrada para círculo
    círculo 2 abajo círculo2 arriba
    círculo 1 abajo [-1,-1] [-1,1]
    círculo 1 arriba [ 1,-1] [ 1,1]

    El uso de cada combinación se estrablece en el vector de 1 y 0 con el siguiente orden:

    sector = [ abajo1*abajo2,  abajo1*arriba2,
              arriba1*abajo2, arriba1*arriba2]

    las instrucciones en Python para lo descrito se muestran como una función:

    import numpy as np
    import scipy.optimize as sp
    def cruce2circulos(x1,y1,r1,x2,y2,r2):
        ''' Revisa intervalo de area de cruce
            entre dos círculos de centro y radio
            x1,y1,r1 // x2,y2,r2
        '''
        intersecta = []
        dx = x2 - x1
        dy = y2 - y1
        d_centros = np.sqrt(dx**2 + dy**2)
        d_cruce   = r2 + r1
        
        # los circulos se cruzan o tocan
        if d_cruce >= d_centros:
    
            # intervalos de cruce
            xa = np.max([x1-r1,x2-r2])
            xb = np.min([x1+r1,x2+r2])
            ya = np.max([y1-r1,y2-r2])
            yb = np.min([y1+r1,y2+r1])
            
            # cada circulo arriba, abajo
            abajo1 = 0 ; arriba1 = 0
            abajo2 = 0 ; arriba2 = 0
            if ya<=y1:
                abajo1  = 1
            if yb>=y1:
                arriba1 = 1
            if ya<=y2:
                abajo2  = 1
            if yb>=y2:
                arriba2 = 1
            sector  = [ abajo1*abajo2, abajo1*arriba2,
                       arriba1*abajo2, arriba1*arriba2]
            uncruce = [xa,xb,ya,yb,sector]
        return(uncruce)
    

    El resultado para los círculos del ejercicio son:

    >>> x1=4; y1=4; r1=np.sqrt(5)
    >>> x2=0; y2=0; r2=np.sqrt(16)
    >>> uncruce = cruce2circulos(x1,y1,r1,x2,y2,r2)
    >>> uncruce
    [1.7639320225002102, 4.0, 
     1.7639320225002102, 2.23606797749979, 
    [0, 1, 0, 0]]
    >>> 
    

    2.2 Raíces como coordenadas de intersección entre dos círculos

    Las coordenadas de intersección entre dos círculos se obtienen aplicando un método de búsqueda de raíces. Por ejemplo bisección, que para esta parte se usa el algoritmo de SciPy con la instrucción sp.bisect(fx,xa,xb,xtol=2e-12).

    Para el caso más general, donde existen dos raíces que buscar, se divide el intervalo de búsqueda [xa,xb] en dos medios segmentos [xa,xc] y [xc,xb]. Se aplica un método de búsqueda de raíces para cada subintervalo. Para minimizar errores de truncamiento, en cada búsqueda de desplaza dx/10 cada xc hacia el lado que amplia el subintervalo de búsqueda.

    intersecta Circulos 01

    Para el caso donde los círculos solo tienen un punto de contacto, se realiza una revisión considerando que el intervalo de búsqueda podría ser menor al valor de tolerancia del radio.

    Por ejemplo, cuando la linea que une los centros de los círculos resulta paralelos al eje de las x,  adicionalmente se topan en un solo punto, el algoritmo anterior indica que se usan todos los sectores de los círculos, dando como resultado cuatro raíces iguales. El caso se corrige realizando la parte de sectores solo cuando la distancia entre [xa,xb] es mayor a cero.

    El resultado se presenta como los vectores raizx y raizy.

    Las instrucciones en Python para esta sección se describen a continuación:

    def raices2circulos(x1,y1,r1,x2,y2,r2,tolera=2e-12):
        ''' busca las intersección entre 2 circulos
            de centro y radio: x1,y1,r1 || x2,y2,r2
            revisa con cruce2circulos()
        '''
        uncruce = cruce2circulos(x1,y1,r1,x2,y2,r2)
        raizx = []; raizy = []
    
        # si hay cruce de circulos
        if len(uncruce)>0:
            sectores = [[-1,-1],[-1,1], 
                        [ 1,-1],[ 1,1]]
            [xa,xb,ya,yb,sector] = uncruce
            xc = (xa+xb)/2
            dx = np.abs(xb-xa)
            dy = np.abs(yb-ya)
            k = 1    # se tocan en un punto
            if dx>0: # se tocan en mas de un punto
                k = len(sector)
            for j in range(0,k,1):
                if sector[j]==1:
                    s1 = sectores[j][0]
                    s2 = sectores[j][1]
                    fx1 = lambda x: s1*np.sqrt(r1**2-(x-x1)**2)+y1
                    fx2 = lambda x: s2*np.sqrt(r2**2-(x-x2)**2)+y2
                    fx  = lambda x: fx1(x)-fx2(x)
                    fa = fx(xa)
                    fb = fx(xb)
                    raiz1 = np.nan
                    raiz2 = np.nan
                    
                    # intervalo/2 izquierda
                    xc = xc + dx/10
                    fc = fx(xc)
                    cambio = np.sign(fa)*np.sign(fc)
                    if cambio<0:
                        raiz1 = sp.bisect(fx,xa,xc,xtol=tolera)
                        
                    # intervalo/2 derecha
                    xc = xc - 2*dx/10
                    fc = fx(xc)
                    cambio = np.sign(fc)*np.sign(fb)
                    if cambio<0:
                        raiz2 = sp.bisect(fx,xc,xb,xtol=tolera)
                        
                    # si hay contacto en un borde
                    if dx<tolera*r1 and dy>0:
                        raiz1 = xa
                    if dy<tolera*r1 and dx>0:
                        raiz1 = x1
                        
                    # Añade si existe raiz
                    if not(np.isnan(raiz1)):
                        raizx.append(raiz1)
                        raizy.append(fx1(raiz1))
                    if not(np.isnan(raiz2)):
                        raizx.append(raiz2)
                        raizy.append(fx1(raiz2))
            raices = [raizx,raizy]
        return(raices)
    

    El resultado del algoritmo para el ejercicio es:

    >>> raices = raices2circulos(x1,y1,r1,x2,y2,r2,tolera=2e-12)
    >>> raices
    [[1.805829001269906, 3.569170998730207],
     [3.569170998734088, 1.8058290012706681]]
    >>>
  • s3Eva_IT2017_T3 Sustancia en lago

    Ejercicio: 3Eva_IT2017_T3 Sustancia en lago

    El ejercicio se divide en dos partes: sección transversal con la derivada y concentración promedio con integrales.

    Sección transversal

    Se calcula la derivada con  una aproximación básica con error O(h)

    f'(x_i) = \frac{f(x_{i+1})-f(x_i)}{h} + O(h)

    repidiendo la fórmula entre cada par de puntos consecutivos

    dv/dz: [-1.1775  -0.7875  -0.39175 -0.09825  0.     ]
    

    Concentración promedio

    Para los integrales usamos la regla del trapecio:

    I = (b-a) \frac{f(a)+f(b)}{2}
    numerador:  224.38960000000003
    denominador:  29.852
    concentracion promedio:  7.516735897092323
    

    Aplicando los algoritmos en Python para todos los puntos:

    # 3Eva_IT2017_T3 Sustancia en lago
    import numpy as np
    import matplotlib.pyplot as plt
    
    # INGRESO
    zi = np.array([0.  , 4   , 8   , 12    , 16])
    vi = np.array([9.82, 5.11, 1.96,  0.393,  0.])
    ci = np.array([10.2, 8.5 , 7.4 ,  5.2  ,  4.1])
    
    # PROCEDIMIENTO
    n = len(zi)
    # primera derivada hacia adelante con error O(h)
    dv = np.zeros(n,dtype=float)
    for i in range(0,n-1,1):
        h = zi[i+1]-zi[i]
        dv[i]=(vi[i+1]-vi[i])/h
    
    As = -dv*zi
    
    # integrales por rectángulo
    numerador = 0
    for i in range(0,n-1,1):
        altura = (ci[i]*As[i]+ci[i+1]*As[i+1])/2
        numerador = numerador +altura*(zi[i+1]-zi[i])
    
    denominador = 0
    for i in range(0,n-1,1):
        altura = (As[i]+As[i+1])/2
        denominador = denominador +altura*(zi[i+1]-zi[i])
    
    cpromedio = numerador/denominador
    
    # SALIDA
    print('dv/dz: ')
    print(dv)
    print('numerador: ',numerador)
    print('denominador: ',denominador)
    print('concentracion promedio: ',cpromedio)
    
    # Grafica
    plt.subplot(121)
    plt.plot(zi,vi)
    plt.plot(zi,vi,'bo')
    plt.xlabel('profundidad z')
    plt.ylabel('Volumen')
    plt.grid()
    plt.subplot(122)
    plt.plot(zi,ci, color = 'orange')
    plt.plot(zi,ci,'ro')
    plt.xlabel('profundidad z')
    plt.ylabel('concentración')
    plt.grid()
    plt.show()
    

  • s3Eva_IT2017_T4 EDP elíptica, placa desplazada

    Ejercicio: 3Eva_IT2017_T4 EDP elíptica, placa desplazada

    La ecuación del problema en forma contínua:

    \frac{\delta ^2 U}{\delta x^2} + \frac{\delta ^2 U}{\delta y^2} = \frac{x}{y} + \frac{y}{x}

    1 <  x < 2
    1 <  y < 2

    Se convierte a la versión discreta usando diferencias divididas centradas


    \frac{u[i-1,j]-2u[i,j]+u[i+1,j]}{\Delta x^2} + + \frac{u[i,j-1]-2u[i,j]+u[i,j+1]}{\Delta y^2} = \frac{x_i}{y_j} + \frac{y_j}{x_i}

    Se agrupan los términos Δx, Δy semejante a formar un λ al multiplicar todo por Δy2

    \frac{\Delta y^2}{\Delta x^2}\Big(u[i-1,j]-2u[i,j]+u[i+1,j] \Big) + + \frac{\Delta y^2}{\Delta y^2}\Big(u[i,j-1]-2u[i,j]+u[i,j+1]\Big) = =\Delta y^2\Big( \frac{x_i}{y_j} + \frac{y_j}{x_i}\Big)
    \lambda= \frac{\Delta y^2}{\Delta x^2} = 1

    por ser los tamaños de paso iguales en ambos ejes, se simplifica la ecuación a usar:


    u[i-1,j]-2u[i,j]+u[i+1,j] + + u[i,j-1]-2u[i,j]+u[i,j+1] = =\Delta y^2\Big( \frac{x_i}{y_j} + \frac{y_j}{x_i}\Big)
    u[i-1,j]-4u[i,j]+u[i+1,j] + + u[i,j-1]+u[i,j+1] =\Delta y^2\Big( \frac{x_i}{y_j} + \frac{y_j}{x_i}\Big)

    Por simplicidad se usará el método iterativo en el ejercicio, por lo que se despeja la ecuación del centro del rombo formado por los puntos,


    4u[i,j] = u[i-1,j]+u[i+1,j] + + u[i,j-1]+u[i,j+1] -\Delta y^2\Big( \frac{x_i}{y_j} + \frac{y_j}{x_i}\Big)
    u[i,j] = \frac{1}{4}\Big( u[i-1,j]+u[i+1,j] + + u[i,j-1]+u[i,j+1] -\Delta y^2\Big( \frac{x_i}{y_j} + \frac{y_j}{x_i}\Big)\Big)

    Iteraciones:

    Se utiliza una matriz de ceros para la iteración inicial. En el ejercicio se muestran cálculos para 3 nodos, el resto se realiza con el algoritmo en Python.

    Para varias iteraciones se usa Δx =Δy = 1/4 = 0.25

    y las ecuaciones para los valores en las fronteras o bordes de la placa

    U(x,1)= x \ln (x), U(x,2) = x \ln (4x^{2}),1 \lt x \lt 2 U(1,y)= y \ln(y), U(2,y) = 2y \ln (2y), 1 \lt x \lt 2

    i=1, j=1


    u[1,1] = \frac{1}{4}\Big( u[0,1]+u[2,1] + + u[1,0]+u[1,2] -(0.25)^2\Big( \frac{1.25}{1.25} + \frac{1.25}{1.25}\Big)\Big)
    u[1,1] = \frac{1}{4}\Big(1.25 \ln (1.25)+0 + + 1.25 \ln(1.25) + 0 -(0.25)^2\Big( \frac{1.25}{1.25} + \frac{1.25}{1.25}\Big)\Big)

    i = 2, j =1


    u[2,1] = \frac{1}{4}\Big( u[1,1]+u[3,1] + + u[2,0]+u[2,2] -(0.25)^2\Big( \frac{1.5}{1.5} + \frac{1.5}{1.5}\Big)\Big)

    Método iterativo

    usando el método iterativo se obtiene los siguientes resultados:

    iteraciones:  15
    error entre iteraciones:  6.772297286980838e-05
    solución para u: 
    [[0.         0.27892944 0.60819766 0.97932763 1.38629436]
     [0.27892944 0.69781162 1.1792239  1.7127402  2.29072683]
     [0.60819766 1.1792239  1.8252746  2.53384036 3.29583687]
     [0.97932763 1.7127402  2.53384036 3.42800537 4.38467039]
     [1.38629436 2.29072683 3.29583687 4.38467039 5.54517744]]
    >>> 
    

    Algoritmo en Python

    # 3Eva_IT2017_T4 EDP elíptica, placa desplazada
    # método iterativo
    import numpy as np
    
    # INGRESO
    # longitud en x
    a = 1
    b = 2
    # longitud en y
    c = 1
    d = 2
    # tamaño de paso
    dx = 0.25
    dy = 0.25
    # funciones en los bordes de la placa
    abajo     = lambda x,y: x*np.log(x)
    arriba    = lambda x,y: x*np.log(4*(x**2))
    izquierda = lambda x,y: y*np.log(y)
    derecha   = lambda x,y: 2*y*np.log(2*y)
    # función de la ecuación
    fxy = lambda x,y: x/y + y/x
    
    # control de iteraciones
    maxitera = 100
    tolera = 0.0001
    
    # PROCEDIMIENTO
    # tamaño de la matriz
    n = int((b-a)/dx)+1
    m = int((d-c)/dy)+1
    # vectores con valore de ejes
    xi = np.linspace(a,b,n)
    yj = np.linspace(c,d,m)
    # matriz de puntos muestra
    u = np.zeros(shape=(n,m),dtype=float)
    
    # valores en los bordes
    u[:,0]   = abajo(xi,yj[0])
    u[:,m-1] = arriba(xi,yj[m-1])
    u[0,:]   = izquierda(xi[0],yj)
    u[n-1,:] = derecha(xi[n-1],yj)
    
    # valores interiores
    # para menos iteraciones
    mitadx = int(n/2)
    mitady = int(m/2)
    promedio = (u[mitadx,0]+u[mitadx,m-1]+u[0,mitady]+u[n-1,mitady])/4
    u[1:n-1,1:m-1] = promedio
    
    # método iterativo
    itera = 0
    converge = 0
    while not(itera>=maxitera or converge==1):
        itera = itera +1
        # copia u para calcular errores entre iteraciones
        nueva = np.copy(u)
        for i in range(1,n-1):
            for j in range(1,m-1):
                # usar fórmula desarrollada para algoritmo
                u[i,j] = (u[i-1,j]+u[i+1,j]+u[i,j-1]+u[i,j+1]-(dy**2)*fxy(xi[i],yj[j]))/4 
        diferencia = nueva-u
        erroru = np.linalg.norm(np.abs(diferencia))
        if (erroru<tolera):
            converge=1
    
    # SALIDA
    print('iteraciones: ',itera)
    print('error entre iteraciones: ',erroru)
    print('solución para u: ')
    print(u)
    
    # Gráfica
    import matplotlib.pyplot as plt
    from matplotlib import cm
    from mpl_toolkits.mplot3d import Axes3D
    
    # matrices de ejes para la gráfica 3D
    X, Y = np.meshgrid(xi, yj)
    U = np.transpose(u) # ajuste de índices fila es x
    figura = plt.figure()
    
    grafica = Axes3D(figura)
    grafica.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()