s2Eva_IIT2018_T1 Masa entra o sale de un reactor

Ejercicio: 2Eva_IIT2018_T1 Masa entra o sale de un reactor

a) Se pueden combinar los métodos para realizar la integral. Se usa el método de Simpson 1/3 para los primeros dos tramos y Simpson 3/8 para los 3 tramos siguientes.  Siendo f(x) equivalente a Q(t)C(t). El tamaño de paso h es constante para todo el ejercicio con valor 5.

a.1 Simpson 1/3, tramos 2, puntos 3:

I_1 \cong \frac{h}{3}[f(x_0)+4f(x_1) + f(x_2)] I_1 \cong \frac{5}{3}[(10)(4)+4(18)(6) + (27)(7)] I_1 \cong 1101,66

a.2 Simpson de 3/8, tramos 3, puntos 4:

I_2 \cong \frac{3h}{8}[f(x_0)+3f(x_1) +3 f(x_2)+f(x_3)] I_2 \cong \frac{3(5)}{8}[(27)(7)+3(35)(6) +3(40)(5)+(30)(5)] I_2 \cong 2941,88 I_1 + I_2 \cong 4043,54

b) El error se calcula por tramo y se acumula.

b.1 se puede estimar como la diferencia entre la parábola del primer tramo y simpson 1/3
b.2 siguiendo el ejemplo anterior, como la diferencia entre la interpolación de los tramos restantes y simpson 3/8.

s2Eva_IT2010_T2 Movimiento angular

Ejercicio: 2Eva_IT2010_T2 Movimiento angular

Para resolver, se usa Runge-Kutta_fg de segundo orden como ejemplo

y'' + 10 \sin (y) =0

se hace

y' = z = f(t,y,z)

y se estandariza:

y'' =z'= -10 \sin (y) = g(t,y,z)

teniendo como punto de partida t0=0, y0=0 y z0=0.1

y(0)=0, y'(0)=0.1

Se desarrolla el algotitmo para obtener los valores:

 [ t, 		 y, 	 dyi/dti=z]
[[ 0.          0.          0.1       ]
 [ 0.2         0.02        0.08000133]
 [ 0.4         0.03200053  0.02401018]
 [ 0.6         0.03040355 -0.04477916]
 [ 0.8         0.01536795 -0.09662411]
 [ 1.         -0.00703034 -0.10803459]]

que permiten generar la gráfica de respuesta:


Algoritmo en Python

# 2Eva_IT2010_T2 Movimiento angular
import numpy as np
import matplotlib.pyplot as plt

def rungekutta2_fg(f,g,x0,y0,z0,h,muestras):
    tamano = muestras + 1
    estimado = np.zeros(shape=(tamano,3),dtype=float)
    # incluye el punto [x0,y0,z0]
    estimado[0] = [x0,y0,z0]
    xi = x0
    yi = y0
    zi = z0
    for i in range(1,tamano,1):
        K1y = h * f(xi,yi,zi)
        K1z = h * g(xi,yi,zi)
        
        K2y = h * f(xi+h, yi + K1y, zi + K1z)
        K2z = h * g(xi+h, yi + K1y, zi + K1z)

        yi = yi + (K1y+K2y)/2
        zi = zi + (K1z+K2z)/2
        xi = xi + h
        
        estimado[i] = [xi,yi,zi]
    return(estimado)

# INGRESO theta = y
ft = lambda t,y,z: z
gt = lambda t,y,z: -10*np.sin(y)

t0 = 0
y0 = 0
z0 = 0.1
h=0.2
muestras = 5

# PROCEDIMIENTO
tabla = rungekutta2_fg(ft,gt,t0,y0,z0,h,muestras)

# SALIDA
print(' [ t, \t\t y, \t dyi/dti=z]')
print(tabla)

# Grafica
ti = np.copy(tabla[:,0])
yi = np.copy(tabla[:,1])
zi = np.copy(tabla[:,2])
plt.subplot(121)
plt.plot(ti,yi)
plt.xlabel('ti')
plt.title('yi')
plt.subplot(122)
plt.plot(ti,zi, color='green')
plt.xlabel('ti')
plt.title('dyi/dti')
plt.show()

s1Eva_IIT2018_T4 Tasa de interés en hipoteca

Ejercicio: 1Eva_IIT2018_T4 Tasa de interés en hipoteca

literal a

Siguiendo el desarrollo analítico tradicional, para adecuar la ecuación para los algoritmo de búsquda de raíces de ecuaciones,  se reemplazan los valores en la fórmula.

P = A\Big(\frac{1-(1+i)^{-n}}{i} \Big) 70000 = 1200\Big(\frac{1-(1+i)^{-300}}{i} \Big)

Como ambos lados de la ecuación deben ser iguales, si se restan ambos se obtiene una ecuación que tiene como resultado cero, que es la forma ideal para usar en el algoritmo que representa f(x) o en este caso f(i)

70000 - 1200\Big(\frac{1-(1+i)^{-300}}{i} \Big) = 0

Para evitar inconvenientes con la división para cero en caso que i tome el valor de cero, dado se multiplica toda la ecuación por i:

i \Big[70000 - 1200\Big(\frac{1-(1+i)^{-300}}{i} \Big) \Big]= i (0) 70000 i - 1200 (1-(1+i)^{-300}) = 0

La ecuación es la utilizada en el algoritmo de búsqueda de raíces pueden ser:

fx(i) = 70000 - 1200\Big(\frac{1-(1+i)^{-300}}{i} \Big) fx(i) = 70000i - 1200(1-(1+i)^{-300})

literal b

El intervalo de existencia correspondería a la tasa de interés mínimo y el interés máximo.

[izquierda, derecha] = [a,b]

Para el intervalo se deben tomar en cuenta algunas consideraciones descritas a continuación:

izquierda:

En el extremo izquierdo, las tasas no son negativas, lo que se interpreta en que un banco paga por que le presten dinero.

Tampoco tiene mucho sentido el valor cero, que son prestamos sin intereses. A menos que sean sus padres quienes le dan el dinero.

Un valor inicial para el interés puede ser por ejemplo 1% ó 0.01, siempre que se cumpla que existe cambio de signo en la función a usar.

derecha:

En el extremo derecho, si se propone por ejemplo i con 100%, o 1.00, no tendría mucho sentido un préstamo con intereses al 100% anual, que resulta en el doble del valor inicial en tan solo un periodo o año.

La tasa de interés de consumo que son de las más alto valor, se encuentran reguladas. En Ecuador es un valor alrededor del 16% anuales o 0.16.

Considerando las observaciones iniciales del problema, se propone empezar el análisis para la búsqueda de la raíz en el intervalo en un rango más amplio:

[ 0.01, 0.50]

Ser realiza la comprobación que existe cambio de signo en los extremos del intervalo.

fx(0.001) =- 43935.86

fx(0.50) = 67600.0

Para el ejercicio se hace notar que la es tasa nominal anual, pero los pagos son mensuales. Por lo que se debe unificar las tasas de interes a mensuales. Una aproximación es usar las tasas anuales divididas para los 12 meses del año.

Tolerancia al error

La tolerancia se considera en éste ejercicio como el valor de diferencias  (tramo) entre iteraciones con precisión satisfactoria.

Por ejemplo si no negociaremos más con el banco por variaciones de tasas del 0.1% , entonces la tolerancia será de 0.001.

Las publicaciones de tasas en el mercado incluyen dos decimales, por lo que para el ejercicio aumentamos la precisión a : 0.0001

tolera = 1×10-4


Literal c


Se presentan dos formas se solución para el litera c:

– c.1 la requerida en el enunciado con Newton-Raphson

– c.2 una alterna con el método de la Bisección.


c.1. Desarrollo del ejercicio con el método del enunciado Newton-Raphson


Para el método de Newton-Raphson se tiene que:

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

Se requiere la derivada de la función planteada en el literal a:

fx(i) = 70000i - 1200(1-(1+i)^{-300}) f'x(i) = 70000 + 1200(300)(1+i)^{-301})

tomando como valor inicial xi = 0.16/12 ≈ 0.013

Se realizan las iteraciones suponiendo que tolera = 1×10-4

iteración 1

fx(0.013) = 70000(0.013) - 1200(1-(1+0.013)^{-300})

 = -265.0914

f'x(0.013) = 70000 + 1200(300)(1+0.013)^{-301})

= 62623.3454

x_{2} = 0.013 - \frac{-265.0914}{63318.8456} = 0.01723

error = |0.013 – 0.01723| = 0.004229

iteración 2

fx(0.01723) = 70000i - 1200(1-(1+0.0.01723)^{-300})

= 13.2356

f'x(0.01723) = 70000 + 1200(300)(1+0.01723)^{-301}

= 67895.5656

x_{3} = 0.01723 - \frac{13.2356}{67895.5656} = 0.01703

error = |0.01723 – 0.01703| = 0.0001999

cuyo valr de error está casi dentro del valor de tolerancia,

que permite tomar el último valor como respuesta de tasa mensual

raiz = tasa mensual = 0.01703

Convirtiendo a la tasa tasa anual que es la publicada por las instituciones financieras se tiene que:

tasa anual nominal =  0.01703*12 = 0.2043

Teniendo como resultado una tasa anual de 20.43%


Algoritmo en Python

# 1ra Evaluación II Término 2018
# Tema 4. Tasa de interes para hipoteca
import numpy as np
import matplotlib.pyplot as plt

def newton_raphson(fx, dfx, xi, tolera):
    '''
    funciónx y fxderiva son de forma numérica
    xi es el punto inicial de búsqueda
    '''
    tramo = abs(2*tolera)
    while (tramo>=tolera):
        xnuevo = xi - fx(xi)/dfx(xi)
        tramo = abs(xnuevo-xi)
        print(xi,xnuevo, tramo)
        xi = xnuevo
    return(xi)

# INGRESO
P = 70000.00
A = 1200.00
n = 25*12
fx  = lambda i: P*i - A*(1-(1+i)**(-n))
dfx = lambda i: P + A*(-n)*(i+1)**(-n-1)
a = 0.01/12
b = 0.25/12
muestras = 51
tolera = 0.0001

# PROCEDIMIENTO
tasa = np.linspace(a,b,muestras)
fi   = fx(tasa)

raiz   = newton_raphson(fx, dfx, b, tolera)
tanual = 12*raiz

# SALIDA
print('raiz encontrada en: ', raiz)
print('tasa anual: ',tanual)

# Gráfica
plt.plot(tasa*12,fi)
plt.title('tasa anual de interes para Hipoteca')
plt.xlabel('tasa')
plt.ylabel('fx(tasa)')
plt.axhline(0, color='green')
plt.show()

c.2. Desarrollo con el método de la Bisección


Desarrollo Analítico con Bisección

Como parte del desarrollo del ejercicio se presenta las iteraciones para el algoritmo, tradicionalmente realizadas con una calculadora.

fx(i) = 70000 - 1200\Big(\frac{1-(1+i)^{-300}}{i} \Big)

iteración 1

a = 0.01, b = 0.5 c = \frac{a+b}{2} = \frac{0.01+0.5}{2} = 0.255 fx(0.01) = 70000 - 1200\Big(\frac{1-(1+(0.01))^{-300}}{0.01} \Big) = -43935.86 fx(0.255) = 70000 - 1200\Big(\frac{1-(1+(0.255))^{-300}}{0.255} \Big) = 65294.11 fx(0.5) = 70000 - 1200\Big(\frac{1-(1+(0.5))^{-300}}{0.5} \Big) = 67600.0 tramo = 0.5-0.01 =0.49

cambio de signo a la izquierda

a = 0.01, b=0.255

iteración 2

a = 0.01, b = 0.225 c = \frac{a+b}{2} = \frac{0.01+0.225}{2} = 0.1325 fx(0.01) = -43935.86 fx(0.1325) = 70000 - 1200\Big(\frac{1-(1+(0.1325))^{-300}}{0.1325} \Big) = 60943.39 fx(0.225) = 65294.11 tramo = 0.225-0.01 =0.215

cambio de signo a la izquierda

a = 0.01, b=0.1325

iteración 3

a = 0.01, b = 0.1325 c = \frac{a+b}{2} = \frac{0.01+0.1325}{2} = 0.07125 fx(0.01) = -43935.86 fx(0.07125) = 70000 - 1200\Big(\frac{1-(1+(0.07125))^{-300}}{0.07125} \Big) = 53157.89 fx(0.1325) = 60943.39 tramo = 0.1325-0.01 =0.1225

cambio de signo a la izquierda

a = 0.01, b=0.07125

y se continuaría con las iteraciones, hasta cumplir que tramo<=tolera

Tabla de datos obtenidos

tabla para Bisección
i a c b f(a) f(c) f(b) tramo
1 0.01 0.255 0.5 -43935.86 65294.11 67600.0 0.49
2 0.01 0.1325 0.255 -43935.86 60943.39 65294.11 0.215
3 0.01 0.07125 0.1325 -43935.86 53157.89 60943.39 0.1225

hasta lo calculado la raiz se encontraría en el intervalo [0.01,0.07125] con error esitmado de 0.06125, aún por mejorar con más iteraciones.

Algoritmo en Python para Bisección

  • El algoritmo bisección usa las variables a y b, por lo que los limites en el intervalo usados son [La,Lb]
  • para el problema la variable ‘i’ se usa en el eje x.
  • La selección de cambio de rango [a,b] se hace usando solo el signo del valor.
  • El algoritmo presentado es tal como se explica en la parte conceptual

Se deja como tarea convertir el algoritmo a funcion def-return de Python.

# 1Eva_IIT2018_T4 Tasa de interés en hipoteca
import numpy as np
import matplotlib.pyplot as plt

# INGRESO
P = 70000.00
A = 1200.00
n = 25*12
fi = lambda i: P - A*(1-((1+i)**-n))/i

# Intervalo de observación
# e inicio de Bisección
La = 0.01
Lb = 0.50

tolera = 0.0001 #grafica

muestras = 21

# PROCEDIMIENTO

# Método de Bisección
a = La
b = Lb
c = (a+b)/2
tramo = np.abs(b-a)
while (tramo>tolera):
    fa = fi(a)
    fb = fi(b)
    fc = fi(c)
    cambio = np.sign(fc)*np.sign(fa)
    if (cambio>0):
        a = c
        b = b
    else:   
        b = c
        a = a
    c = (a+b)/2
    tramo = np.abs(b-a)

# Para la gráfica
tasa = np.linspace(La,Lb,muestras)
fr = fi(tasa)

# SALIDA
print('a, f(a):', a,fa)
print('c, f(c):', c,fc)
print('b, f(b):', b,fb)
print('la raiz esta entre: \n',a,b)
print('con un error de: ', tramo)
print('raiz es tasa buscada: ', c)
print('tasas anual buscada: ',c*12)

# Gráfica
plt.plot(tasa,fr)
plt.axhline(0, color='green')
plt.title('tasa de interes mensual')
plt.show()

la ejecución del algoritmo da como resultado

>>> 
 RESTART: D:/MATG1052Ejemplos/HipotecaInteres.py 
a, f(a): 0.016998291015625 -385.52828922150366
c, f(c): 0.0170281982421875 -145.85350695741363
b, f(b): 0.01705810546875 92.28034212642524
la raiz esta entre: 
 0.016998291015625 0.01705810546875
con un error de:  5.981445312500111e-05
raiz es tasa buscada:  0.0170281982421875
tasas anual buscada:  0.20433837890625

y la gráfica obtenida es:

s1Eva_IIT2018_T1 Interpolar velocidad del paracaidista

Ejercicio: 1Eva_IIT2018_T1 Interpolar velocidad del paracaidista

El ejercicio tiene dos partes: la interpolación y el integral.

Literal a

No se especifica el método a seguir, por lo que se puede seleccionar el de mayor preferencia.

Usando Lagrange, con los puntos primero, medio y último:

p_2(t) = 0\frac{(t-4)(t-8)}{(0-4)(0-8)} + + 27.77\frac{(t-0)(t-8)}{(4-0)(4-8)} + + 41.10\frac{(t-0)(t-4)}{(8-0)(8-4)} p_2(t) = 0 + 27.77\frac{t(t-8)}{-16}) + + 41.10\frac{t(t-4)}{32} p_2(t) = -1.73(t^2-8t) + 1.28(t^2-4t) p_2(t) = -0.45 t^2 + 8.72t


Literal b

El tema de integración para primera evaluación se realiza de forma analítica.

Como introducción a la Unidad 7 de la 2da Evaluación, se usa el método del trapecio para el polinomio encontrado en el literal anterior:

v(t) =  -0.45125*t**2 + 8.7475*t
distancia recorrida:  202.8912640000001

Realice la integración analítica del polinomio encontrado y determine el error entre los métodos.


Desarrollo con Python

las instrucciones en Python para el ejercicio son:

# 1ra Evaluación II Término 2018
# Tema 1. Interpolar velocidad del paracaidista

import numpy as np
import matplotlib.pyplot as plt
import sympy as sym

# Literal a)
def interpola_lagrange(xi,yi):
    '''
    Interpolación con método de Lagrange
    resultado: polinomio en forma simbólica
    '''
    # PROCEDIMIENTO
    n = len(xi)
    x = sym.Symbol('x')
    # Polinomio
    polinomio = 0
    for i in range(0,n,1):
        # Termino de Lagrange
        termino = 1
        for j  in range(0,n,1):
            if (j!=i):
                termino = termino*(x-xi[j])/(xi[i]-xi[j])
        polinomio = polinomio + termino*yi[i]
    # Expande el polinomio
    polinomio = polinomio.expand()
    return(polinomio)

# INGRESO
t = [0.0, 2, 4, 6, 8]
v = [0.0, 16.40, 27.77, 35.64, 41.10]

xi = [t[0],t[2],t[4]]
yi = [v[0],v[2],v[4]]

muestras = 51

# PROCEDIMIENTO
polinomio = interpola_lagrange(xi,yi)
velocidad = polinomio.subs('x','t')

# Para graficar
vt = sym.lambdify('t',velocidad)
a = t[0]
b = t[-1]
ti = np.linspace(a, b, muestras)
vi = vt(ti)

# SALIDA
print('v(t) = ', velocidad)
# Grafica
plt.plot(t,v,'ro')
plt.plot(ti,vi)
plt.title('Interpolar velocidad de paracidista')
plt.xlabel('t')
plt.ylabel('v')
plt.show()


# Literal b
def integratrapecio(funcionx,a,b,tramos):
    h = (b-a)/tramos
    x = a
    suma = funcionx(x)
    for i in range(0,tramos-1,1):
        x = x+h
        suma = suma + 2*funcionx(x)
    suma = suma + funcionx(b)
    area = h*(suma/2)
    return(area)

# INGRESO
# El ingreso es el polinomio en forma lambda
# se mantienen las muestras
tramos = muestras-1
# PROCEDIMIENTO
distancia = integratrapecio(vt,a,b,tramos)

# SALIDA
print('distancia recorrida: ', distancia)

s1Eva_IIT2018_T3 Interpolar con sistema de ecuaciones

Ejercicio: 1Eva_IIT2018_T3 Interpolar con sistema de ecuaciones

El tema es semejante al tema 1, cambiando el método de interpolación.


Literal a

Se usan los puntos de las posiciones 0, 3 y 5.
en la fórmula:

p_2(x) = b_0 + b_1x + b_2 x^2

en la fórmula:

punto x[0] = 1, y[0]= 1.84

1.84 = b_0 + b_1(1) + b_2 (1)^2 1.84 = b_0 + b_1 + b_2

punto x[3] = 1.5, y[3]= 2.28

2.28 = b_0 + b_1(1.5) + b_2 (1.5)^2 2.28 = b_0 + 1.5 b_1 + 2.25 b_2

punto x[5] = 2.1, y[5]= 3.28

3.28= b_0 + b_1(2.1) + b_2 (2.1)^2 3.28= b_0 + 2.1 b_1 + 4.41 b_2

se obtiene el sistema de ecuaciones:

b_0 + b_1 + b_2 = 1.84 b_0 + 1.5 b_1 + 2.25 b_2 = 2.28 b_0 + 2.1 b_1 + 4.41 b_2 = 3.28

Con lo que se plantea la forma Ax=B:

A = \begin{bmatrix} 1 & 1 & 1\\ 1 & 1.5 & 2.25 \\1 & 2.1 & 4.41 \end{bmatrix} B = \begin{bmatrix} 1.84\\ 2.28 \\ 3.28 \end{bmatrix}

y se obtiene el resultado de la interpolación.


Literal b

Se requiere calcular una norma de suma de filas. es suficiente para demostrar el conocimiento del concepto el usar A.

Se adjunta el cálculo del número de condición y la solución al sistema de ecuaciones:

suma de columnas:  [3.   4.75 7.51]
norma A:  7.51
numero de condicion:  97.03737354737129
solucion: 
[ 2.03272727 -0.90787879  0.71515152]

El comentario importante corresponde al número de condición, que es un número muy alto para usar un método iterativo, por lo que la solución debe ser un método directo.
Se puede estimar será un número mucho mayor que 1, pues la matriz no es diagonal dominante.


Instrucciones en Python

# 1ra Evaluación II Término 2018
# Tema 3. Interpolar con sistema de ecuaciones

import numpy as np
import matplotlib.pyplot as plt

# --------------------------
# forma matricial para interpolar
A = np.array([[1, 1. , 1.  ],
              [1, 1.5, 2.25],
              [1, 2.1, 4.41]])

B = np.array([1.84, 2.28, 3.28])

# literal b
sumacolumnas = np.sum(A, axis =1)
norma = np.max(sumacolumnas)
print('suma de columnas: ', sumacolumnas)
print('norma A: ', norma)

numerocondicion = np.linalg.cond(A)
print('numero de condicion: ', numerocondicion)

solucion = np.linalg.solve(A,B)
print('solucion: ')
print(solucion)

s1Eva_IIT2018_T2 Distancia mínima a un punto

Ejercicio: 1Eva_IIT2018_T2 Distancia mínima a un punto

Literal a

Se requiere analizar la distancias entre una trayectoria y el punto = [1,1]

Al analizar las distancias de ex y el punto [1,1] se trazan lineas paralelas a los ejes desde el punto [1,1], por lo que se determina que el rango de x = [a,b] para distancias se encuentra en:

a > 0, a = 0.1
b < 1, b = 0.7

El ejercicio usa la fórmula de distancia entre dos puntos:

d = \sqrt{(x_2-x_1)^2+(y_2- y_1)^2}

en los cuales:

[x1,y1] = [1,1]
[x2,y2] = [x, ex]

que al sustituir en la fórmula se convierte en:

d = \sqrt{(x-1)^2+(e^x- 1)^2}

que es lo requerido en el literal a


Literal b

Para encontrar el punto más cercano, se debe encontrar el mínimo de la distancia, se podría derivar la función y encontrar la raiz en cero.

Considere simplificar la función a un polinomio, donde tiene dos opciones:

b.1 Polinomio de Taylor, que también requiere derivadas (descartado)

b.2 evaluar la función en varios puntos, interpolar y obtener un polinomio.

Al reutilizar el algoritmo del tema 1 se obtiene lo planteado en b.2, usando un polinomio de grado 3, con muestras de 4 puntos equidistantes en el eje x.

polinomio
0.867192074184622*x**3 + 1.22015957396232*x**2 - 1.21861610672236*x + 1.01491694350023
derivada polinomio:
2.60157622255387*x**2 + 2.44031914792464*x - 1.21861610672236

Al aplicar un método para encontrar raíces se tiene que:

en distancia mínima en x= 0.3644244280922699
con y = 1.4396851273165785
distancia =  0.7728384816953889


Instrucciones en Python

# 1ra Evaluación II Término 2018
# Tema 2. Distancia mínima a un punto
import numpy as np
import matplotlib.pyplot as plt
import sympy as sym

# PRESENTA PROBLEMA
# INGRESO
punto = [1,1]
trayectoria = lambda x: np.exp(x)

a = 0
b = 1
muestras = 51

# PROCEDIMIENTO
xif = np.linspace(a,b,muestras)
trayecto = trayectoria(xif)

# SALIDA
plt.plot(xif,trayecto, label = 'trayectoria')
plt.plot(punto[0],punto[1],'ro')
plt.axhline(punto[0], color='grey')
plt.axvline(punto[1], color='grey')
plt.title('distancia a un punto')
plt.legend()
plt.xlabel('x')
plt.ylabel('y')
plt.show()

# ------------------------
# PARA ANALIZAR DISTANCIAS

def interpola_lagrange(xi,yi):
    '''
    Interpolación con método de Lagrange
    resultado: polinomio en forma simbólica
    '''
    # PROCEDIMIENTO
    n = len(xi)
    x = sym.Symbol('x')
    # Polinomio
    polinomio = 0
    for i in range(0,n,1):
        # Termino de Lagrange
        termino = 1
        for j  in range(0,n,1):
            if (j!=i):
                termino = termino*(x-xi[j])/(xi[i]-xi[j])
        polinomio = polinomio + termino*yi[i]
    # Expande el polinomio
    polinomio = polinomio.expand()
    return(polinomio)

def posicionfalsa(fx,a,b,tolera):
    fa = fx(a)
    fb = fx(b)
    c = b - fb*(a-b)/(fa-fb)
    
    tramo = abs(c-a)
    while (tramo > tolera):
        fc = fx(c)
        cambia = np.sign(fa)*np.sign(fc)
        if (cambia > 0):
            tramo = abs(c-a)
            a=c
            fa=fc
        else:
            tramo = abs(c-b)
            b=c
            fb=fc
        c = b - fb*(a-b)/(fa-fb)
        print(tramo)
    respuesta = c
    
    # Valida respuesta
    fa = fx(a)
    fb = fx(b)
    cambio = np.sign(fa)*np.sign(fb)
    if (cambio>0):
        respuesta = np.nan
        
    return(respuesta)

# INGRESO
# Trayectoria y punto tomados de sección PRESENTA PROBLEMA
y = lambda x: np.sqrt((x-punto[0])**2+(trayectoria(x)-punto[1])**2)

a = 0.1
b = 0.7
muestras = 51

tolera = 0.0001
muestrap = 4 # Para el polinomio

# PROCEDIMIENTO
yif = y(xif)
xip = np.linspace(a,b,muestrap)
yip = y(xip)

x = sym.Symbol('x')

polinomio = interpola_lagrange(xip,yip)

px = sym.lambdify('x',polinomio)
pxi = px(xif)

dpx = polinomio.diff('x',1)
dpxn = sym.lambdify('x',dpx)
fxi = dpxn(xif)

raiz = posicionfalsa(dpxn,a,b,tolera)

# SALIDA
print('polinomio')
print(polinomio)
print('derivada polinomio:')
print(dpx)
print('en distancia mínima en x=',raiz)
print('con y =', trayectoria(raiz))
print('distancia = ', y(raiz))

# GRAFICA
plt.plot(xif,yif, label = 'distancia')
plt.plot(xif,pxi, label = 'p(x)')
plt.plot(xif,fxi, label = 'dpx(x)')
plt.axhline(0, color = 'grey')
plt.axvline(raiz,  color = 'grey')
plt.legend()
plt.title('Distancia mínima de punto a f(x)')
plt.xlabel('x')
plt.ylabel('distancia')
plt.show()

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]

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.

De la gráfica se usa el ‘zoom’ y se puede aproximar los valores para la intersección de las curvas estimando raices 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 ecuacion 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 raices 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


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

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 ídentificar 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 raiz 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 busqueda [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 busqueda de desplaza dx/10 cada xc hacia el lado que amplia el subintervalo de búsqueda.

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,  adicionalmete se topan en un solo punto, el algoritmo anterior indica que se usan todos los sectores de los círculos, dando como resultado cuatro raices 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 intrucciones 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]]
>>>

s2Eva_IT2018_T3 EDP Eliptica

Ejercicio: 2Eva_IT2018_T3 EDP Eliptica

Generar las ecuaciones a resolver usando diferencias finitas divididas centradas:

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

por facilidad se sustituye también en la forma discreta de la ecuación:

f[i,j] = f(x_i,y_j) = 2(x_{i} ^2+y_{j} ^2)
\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} = f[i,j]
\frac{\Delta y^2}{\Delta x^2}\Big(u[i-1,j]-2u[i,j]+u[i+1,j]\Big) + + u[i,j-1]-2u[i,j]+u[i,j+1] = \Delta y^2 f[i,j]

dado que hx = hy = 1/3

\frac{\Delta y^2}{\Delta x^2} = 1
(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 f[i,j]
u[i-1,j]-4u[i,j]+u[i+1,j] + + u[i,j-1]+u[i,j+1] = \Delta y^2 f[i,j]
4u[i,j] = u[i-1,j]+u[i+1,j] + + u[i,j-1]+u[i,j+1]-\Delta y^2 f[i,j]
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 f[i,j] \Big)

La última ecuación puede ser usada de forma iterativa, para lo cual hay que definir los valores iniciales de la matriz u.

Al conocer el rango de operación para los ejes x, y, hx, hy se realizan los cálculos para:

1. Evaluar los valores para cada eje x[i], y[j]

x[i]:
[ 0.    0.33  0.67  1.  ]
y[j]:
[ 0.    0.33  0.67  1.  ]

2. Evaluar en cada punto generando una matriz f(i,j):

f[i,j]:
[[ 0.    0.22  0.89  2.  ]
 [ 0.22  0.44  1.11  2.22]
 [ 0.89  1.11  1.78  2.89]
 [ 2.    2.22  2.89  4.  ]]

3. Se evaluan las funciones indicadas para la frontera y se tiene la matriz inicial para u:

matriz inicial u[i,j]:
[[ 1.    1.33  1.67  2.  ]
 [ 1.33  0.    0.    2.44]
 [ 1.67  0.    0.    3.11]
 [ 2.    2.44  3.11  4.  ]]

con lo que se puede trabajar cada punto i,j de forma iterativa, teniendo como resultado para la matriz u:

resultado para u, iterando: 
converge =  1
[[ 1.    1.33  1.67  2.  ]
 [ 1.33  1.68  2.05  2.44]
 [ 1.67  2.05  2.53  3.11]
 [ 2.    2.44  3.11  4.  ]]

La gráfica usando una mayor resolución para tener una idea de la solución:


Los resultados se obtienen usando las siguientes instrucciones:

# 2da Evaluación I Término 2018
# Tema 3. EDP Eliptica
import numpy as np

# INGRESO
# ejes x,y
x0 = 0 ; xn = 1 ; hx = (1/3)# (1/3)/10
y0 = 0 ; yn = 1 ; hy = (1/3) # (1/3)/10
# Fronteras
fux0 = lambda x: x+1
fu0y = lambda y: y+1
fux1 = lambda x: x**2 + x + 2
fu1y = lambda y: y**2 + y + 2

fxy = lambda x,y: 2*(x**2+y**2)

# PROCEDIMIENTO
xi = np.arange(x0,xn+hx,hx)
yj = np.arange(y0,yn+hy,hy)
n = len(xi)
m = len(yj)
# funcion f[xi,yi]
fij = np.zeros(shape=(n,m), dtype = float)
for i in range(0,n,1):
    for j in range(0,m,1):
        fij[i,j]=fxy(xi[i],yj[j])
# matriz inicial u[i,j]
u = np.zeros(shape=(n,m), dtype = float)
u[:,0] = fux0(xi)
u[0,:] = fu0y(yj)
u[:,m-1] = fux1(xi)
u[n-1,:] = fu1y(yj)

uinicial = u.copy()

# Calcular de forma iterativa
maxitera = 100
tolera = 0.0001
# valor inicial de iteración
promedio = (np.max(u)+np.min(u))/2
u[1:n-1,1:m-1] = promedio
# iterar puntos interiores
itera = 0
converge = 0
erroru = 2*tolera # para calcular al menos una matriz
while not(erroru=maxitera):
    itera = itera +1
    nueva = np.copy(u)
    for i in range(1,n-1):
        for j in range(1,m-1):
            u[i,j] = (u[i-1,j]+u[i+1,j]+u[i,j-1]+u[i,j+1]-(hy**2)*fij[i,j])/4
    diferencia = nueva-u
    erroru = np.linalg.norm(np.abs(diferencia))
if (erroru<tolera):
    converge=1

# SALIDA
np.set_printoptions(precision=2)
print('x[i]:')
print(xi)
print('y[j]:')
print(yj)
print('f[i,j]:')
print(fij)
print('matriz inicial u[i,j]:')
print(uinicial)
print('resultado para u, iterando: ')
print('converge = ', converge)
print('iteraciones = ', itera)
print(u)

para obtener la gráfica se debe añadir:

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

X, Y = np.meshgrid(xi, yj)
figura = plt.figure()
ax = Axes3D(figura)
U = np.transpose(u) # ajuste de índices fila es x
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()