s1Eva_IT2016_T3_MN Tasa interés anual

Ejercicio: 1Eva_IT2016_T3_MN Tasa interés anual

Propuesta de Solución  empieza con el planteamiento del problema, luego se desarrolla con el método de Bisección y método del Punto Fijo solo con el objetivo de comparar resultados.


Planteamiento del problema

La fórmula del enunciado para el problema es:

A = P \frac{i(1+i)^{n}}{(1+i)^{n} -1}

que con los datos dados se convierte a:

5800 = 35000 \frac{i(1+i)^8}{(1+i)^8 -1} 35000 \frac{i(1+i)^8}{(1+i)^8 -1}-5800 =0

que es la forma de f(x) = 0

f(i)=35000 \frac{i(1+i)^8}{(1+i)^8 -1}-5800

Intervalo de búsqueda

Como el problema plantea la búsqueda de una tasa de interés, consideramos que:

  • Las tasas de interés no son negativas. ⌉(i<0)
  • Las tasas de interés no son cero en las instituciones bancarias (i≠0)
  • Las tasas muy grandes 1 = 100/100 = 100% tampoco tienen mucho sentido

permite acotar la búsqueda a un intervalo (0,1].
Sin embargo tasas demasiado altas tampoco se consideran en el problema pues el asunto es regulado (superintendencia de bancos), por lo que se podría intentar entre un 1% = 0.01l y la mitad del intervalo 50%= 0.5 quedando

[0.01,0.5]

Tolerancia

si la tolerancia es de tan solo menor a un orden de magnitud que el valor buscado, se tiene que las tasas de interés se representan por dos dígitos después del punto decimal, por lo que la tolerancia debe ser menor a eso.

Por ejemplo: tolerancia < 0.001 o aumentando la precisión

tolera = 0.0001


Método de la Bisección

itera = 1

a = 0.01, b = 0.5

c = \frac{a+b}{2} = \frac{0.01+0.5}{2} = 0.255 f(0.01) = 35000 \frac{0.01(1+0.01)^8}{(1+0.01)^8-1}-5800 = -1225.83 f(0.5) = 35000 \frac{0.5(1+0.5i)^8}{(1+0.5)^8-1}-5800 = 12410.54

con lo que se verifica que existe cambio de signo al evaluar f(x) en el intervalo y puede existir una raíz.

f(0.255) = 35000 \frac{0.255(1+0.255)^8}{(1+0.255)^8-1}-5800 = 4856.70

lo que muestra que f(x) tiene signos en a,c,b de (-) (+) (+), seleccionamos el intervalo izquierdo para continuar la búsqueda [0.01, 0.255]

El tramo permite estimar el error, reduce el intervalo a:

tramo = b-a = 0.255-0.01 = 0.245

valor que todavía es más grande que la tolerancia de 10-4, por lo que hay que continuar las iteraciones.

itera = 2

a = 0.01, b = 0.255

c = \frac{0.01+0.255}{2} = 0.1325 f(0.01) = -1225.83 f(0.255) = 4856.70 f(0.1325 ) = 35000 \frac{0.1325(1+0.1325)^8}{(1+0.1325)^8-1}-5800 = 1556.06

lo que muestra que f(x) tiene signos en a,c,b de (-) (+) (+), seleccionamos el intervalo izquierdo para continuar la búsqueda [0.01, 0.1325]

El tramo permite estimar el error, reduce el intervalo a:

tramo = b-a = 0.1325-0.01 = 0.1225

valor que todavía es más grande que la tolerancia de 10-4, por lo que hay que continuar las iteraciones.

itera = 3

a = 0.01, b = 0.1325

c = \frac{0.01+0.1225}{2} = 0.071 f(0.01) = -1225.83 f(0.1325) = 1556.06 f(0.071 ) = 35000 \frac{0.071(1+0.071)^8}{(1+0.071)^8-1}-5800 = 89.79

lo que muestra que f(x) tiene signos en a,c,b de (-) (+) (+), seleccionamos el intervalo izquierdo para continuar la búsqueda [0.01, 0.071]

El tramo permite estimar el error, reduce el intervalo a:

tramo = b-a = 0.071-0.01 = 0.061

valor que todavía es más grande que la tolerancia de 10-4, por lo que hay que continuar las iteraciones.

Para una evaluación del tema en forma escrita es suficiente para mostrar el objetivo de aprendizaje, el valor final se lo encuentra usando el algoritmo.

       raiz en:  0.06724243164062502
error en tramo:  5.981445312500111e-05
iteraciones:  13
>>>

Algoritmo en Python

# 1Eva_IT2016_T3_MN Tasa interés anual
import numpy as np
import matplotlib.pyplot as plt

# INGRESO
fx = lambda x: 35000*(x*(1+x)**8)/((1+x)**8 -1) -5800
a = 0.01
b = 0.5
tolera = 0.0001

# PROCEDIMIENTO
cuenta = 0
np.set_printoptions(precision=3)
tramo = b-a
while not(tramo<tolera):
    c = (a+b)/2
    fa = fx(a)
    fb = fx(b)
    fc = fx(c)
    cambia = np.sign(fa)*np.sign(fc)
    if cambia < 0: a = a b = c if cambia > 0:
        a = c
        b = b
    tramo = b-a
    cuenta = cuenta+1

# SALIDA
print('       raiz en: ', c)
print('error en tramo: ', tramo)
print('iteraciones: ',cuenta)


Método del Punto Fijo

El planteamiendo del punto fijo se realiza con x= g(x), por lo que se reordena la ecuación a:

35000 \frac{i(1+i)^8}{(1+i)^8 -1}-5800 =0 35000 \frac{i(1+i)^8}{(1+i)^8 -1}=5800 \frac{i(1+i)^8}{(1+i)^8 -1}=\frac{5800}{35000} i=\frac{58}{350} \frac{(1+i)^8 -1} {i(1+i)^8}

con lo que g(x) es:

g(i)=\frac{58}{350} \frac{(1+i)^8 -1} {(1+i)^8}

valor inicial de búsqueda

Para el punto inicial i0 se puede usar uno de los extremos del intervalo propuesto en la sección de planteamiento. Para reducir aun más la búsqueda se pude seleccionar el punto intermedio

 i0 = 0.255

Itera = 1

 i0 = 0.255

g(0.255i)=\frac{58}{350} \frac{(1+0.255)^8 -1} {(1+0.255)^8} = 0.1387

el error se estrima como el tramo recorrido entre el valor nuevo y el valor inicial

tramo = | nuevo-antes| = |0.1387 - 0.255| = 0.1163

como el tramo o error es aún mayor que el valor de tolera, se continúa con la siguiente iteración.

Itera = 2

i1 = g(i0 ) = 0.1387

g(0.1387)=\frac{58}{350} \frac{(1+0.1387)^8 -1} {(1+0.1387)^8} = 0.1071

el error se estrima como el tramo recorrido entre el valor nuevo y el valor inicial

tramo = | nuevo-antes| = |0.1071 - 0.1387| = 0,0316

como el tramo o error es aún mayor que el valor de tolera, se continúa con la siguiente iteración.

Itera = 3

i2 = g(i1 ) = 0.1071

g(0.1071)=\frac{58}{350} \frac{(1+0.1071)^8 -1} {(1+0.1071)^8} = 0.0922

el error se estrima como el tramo recorrido entre el valor nuevo y el valor inicial

tramo = | nuevo-antes| = |0.0922 - 0.1071| = 0,0149

como el tramo o error es aún mayor que el valor de tolera, se continúa con la siguiente iteración.

Observación: Como el error disminuye entre cada iteración, se considera que el método converge, si se realizan suficientes iteraciones se cumplierá con el valor de tolerancia y se habrá llegado a la precisión requerida.

Usando el algoritmo se tiene:

iteraciones: 17
    raiz en: 0.06754389199556779
>>>

Algoritmo en Python

# Algoritmo de Punto Fijo
# x0 valor inicial de búsqueda
# error = tolera

import numpy as np
def puntofijo(gx,antes,tolera, iteramax=50):
    itera = 1 # iteración
    nuevo = gx(antes)
    tramo = abs(nuevo-antes)
    while(tramo>=tolera and itera<=iteramax ):
        antes = nuevo
        nuevo = gx(antes)
        tramo = abs(nuevo-antes)
        itera = itera+1
    respuesta = nuevo
    # Validar respuesta
    if (itera>=iteramax ):
        respuesta = np.nan
    print('iteraciones:',itera)
    return(respuesta)

# PROGRAMA ---------

# INGRESO
gx = lambda i: (58/350)*((1+i)**8-1)/((1+i)**8)

a = 0       # intervalo
b = 0.5

x0 = 0.255
tolera = 0.0001
iteramax  = 50      # itera máximo

# PROCEDIMIENTO
respuesta = puntofijo(gx,0.255,tolera,iteramax)

# SALIDA
print('    raiz en:',respuesta)

s3Eva_2020PAOI_T2 Modelo epidemiológico no letal

El ejercicio representa un sistema de ecuaciones diferenciales ordinarias, que serán resueltas usando Runge-Kutta de 2do Orden.

De compararse con la curva de contagios de Covid-19 se tienen diferencias en la población recuperada, pues el modelo se considera no letal por lo que no se contabiliza el número de fallecidos.

El módelo es el más básico y permite cambiar por ejemplo la tasa de infección, y se ve los cambios en la curva de infectados. Se puede observar lo que se indicaba como objetivo de «aplanar la curva» al disminuir la población expuesta mediante cambiar la tasa de infección al exponer más o menos población al contagio por iteacción entre «suceptibles» e «infectados.


Desarrollo analítico

Las fórmulas para el algoritmo se identifican como:

binfecta = 1.4
grecupera = 1/4
# Ecuaciones
fS = lambda t,S,I,R : -binfecta*S*I
gI = lambda t,S,I,R : binfecta*S*I - grecupera*I
qR = lambda t,S,I,R : grecupera*I

que luego se usan en cada iteración que se registra en la tabla empezando con las condiciones iniciales

itera Si Ii Ri
0 1 0.001 0
1 0.997797107 0.002809143 0.00039375
2 0.9916392093856501 0.007862023500353846 0.0014987671139961277
3 Tarea

itera= 1

K1S = h * fS(ti,Si,Ii,Ri) 
    = 1(-1.4*1*0.001) = -0.0014
K1I = h * gI(ti,Si,Ii,Ri) 
    = 1(1.4*1*0.001 - (1/4)*0.001) = 0.00115
K1R = h * qR(ti,Si,Ii,Ri) 
    = 1((1/4)*0.001)  = 0.00025

K2x = h * fS(ti+h, Si + K1S, Ii+K1I, Ri +K1R) 
    = 1(-1.4*(1-0.0014)(0.001+0.00115)) = -0.003005786
K2y = h * gI(ti+h, Si + K1S, Ii+K1I, Ri +K1R) 
    = 1(1.4*(1-0.0014)*(0.001+0.00115)
        -(1/4)*(0.001+0.00115)) 
    = 0.002468286
K2z = h * qR(ti+h, Si + K1S, Ii+K1I, Ri +K1R) 
    = 1( (1/4)*(0.001+0.00115)) = 0.0005375

Si = Si + (1/2)*(K1S+K2S) 
   = 1 + (-0.0014 -0.003005786)/2 = 0.997797107
Ii = Ii + (1/2)*(K1I+K2I) 
   = 0.001 + (0.00115+0.002468286)/2 = 0.002809143
Ri = Ri + (1/2)*(K1R+K2R) 
   = 0 + (0.00025 + 0.0005375)/2 = 0.00039375
ti = ti + h = 1 + 1 = 2

itera = 2

K1S = 1(-1.4*0.997797107* 0.002809143)
    = -0.003924136661969021
K1I = 1(1.4*0.997797107* 0.002809143 
        - (1/4)*0.002809143)
    = 0.003221850911969021
K1R = 1((1/4)*0.002809143)
    = 0.00070228575

K2S = 1(-1.4*(0.997797107-0.003924136661969021)*
        *(0.002809143+0.003221850911969021))
    = -0.008391658566730926
K2I = 1(1.4*(0.997797107-0.003924136661969021)*
        *(0.002809143+0.003221850911969021) 
- (1/4)*(0.002809143+0.003221850911969021))
    = 0.006883910088738671
K2R = 1( (1/4)*(0.002809143+0.003221850911969021))
    = 0.0015077484779922553

Si = 0.997797107 + 
     + (-0.003924136661969021 -0.008391658566730926)/2
   = 0.9916392093856501
Ii =  0.002809143 + 
      + (0.003221850911969021+0.006883910088738671)/2
   = 0.007862023500353846
Ri = 0.00039375 + 
     + (0.00070228575 + 0.0015077484779922553)/2
   = 0.0014987671139961277
ti = ti + h = 2 + 1 = 3

itera 3 – TAREA


Instrucciones en Python

# 3Eva_2020PAOI_T2 Modelo epidemiológico no letal
# Sistemas EDO con Runge Kutta de 2do Orden
import numpy as np

def rungekutta2_fgq(f,g,q,t0,x0,y0,z0,h,muestras):
    tamano = muestras +1
    tabla = np.zeros(shape=(tamano,4),dtype=float)
    tabla[0] = [t0,x0,y0,z0]
    ti = t0
    xi = x0
    yi = y0
    zi = z0
    for i in range(1,tamano,1):
        K1x = h * f(ti,xi,yi,zi)
        K1y = h * g(ti,xi,yi,zi)
        K1z = h * q(ti,xi,yi,zi)
        
        K2x = h * f(ti+h, xi + K1x, yi+K1y, zi +K1z)
        K2y = h * g(ti+h, xi + K1x, yi+K1y, zi +K1z)
        K2z = h * q(ti+h, xi + K1x, yi+K1y, zi +K1z)

        xi = xi + (1/2)*(K1x+K2x)
        yi = yi + (1/2)*(K1y+K2y)
        zi = zi + (1/2)*(K1z+K2z)
        ti = ti + h
        
        tabla[i] = [ti,xi,yi,zi]
    tabla = np.array(tabla)
    return(tabla)

# Programa
# Parámetros de las ecuaciones

binfecta = 1.4
grecupera = 1/4
# Ecuaciones
fS = lambda t,S,I,R : -binfecta*S*I
gI = lambda t,S,I,R : binfecta*S*I - grecupera*I
qR = lambda t,S,I,R : grecupera*I
# Condiciones iniciales
t0 = 0
S0 = 1.0
I0 = 0.001
R0 = 0.0

# parámetros del algoritmo
h = 1.0
muestras = 51

# PROCEDIMIENTO
tabla = rungekutta2_fgq(fS,gI,qR,t0,S0,I0,R0,h,muestras)
ti = tabla[:,0]
Si = tabla[:,1]
Ii = tabla[:,2]
Ri = tabla[:,3]
# SALIDA
np.set_printoptions(precision=6)
print(' [ ti, Si, Ii, Ri]')
print(tabla)

# Grafica tiempos vs población
import matplotlib.pyplot as plt
plt.plot(ti,Si, label='S')
plt.plot(ti,Ii, label='I')
plt.plot(ti,Ri, label='R')
plt.title('Modelo SIR')
plt.xlabel('t tiempo')
plt.ylabel('población')
plt.legend()
plt.grid()
plt.show()

s1Eva_IIT2009_T1 Movimiento de partícula en plano

a. Planteamiento del problema

Las ecuaciones expresan las trayectorias de dos partículas,

x(t) = 3 \sin ^{3}(t)-1 y(t) = 4 \sin (t)\cos (t)

que para que se encuentren o choquen, sus coordenadas deberían ser iguales.

x(t) = y(t) 3 \sin ^{3}(t)-1 = 4 \sin (t)\cos (t)

Se reordena la igualdad, de tal manera que uno de los lados sea cero, de la forma f(t) = 0 que es usada en el algoritmo de búsqueda de raíces.

3 \sin ^{3}(t)-1 - 4 \sin (t)\cos (t) = 0 f(t) = 3 \sin ^{3}(t)-1 - 4 \sin (t)\cos (t)

b. Intervalo de búsqueda de raíz

Como la variable independiente es tiempo, el evento a buscar se suponte sucede en tiempos positivos t>=0, por lo que el valor inicial a la izquierda del intervalo será a=0

Para el caso de b, a la derecha, se usa lo indicado en el enunciado para la pregunta dell literal b), donde se indica t ∈ [0, π/2], por lo que b = π/2

[0, π/2]

verificando que exista cambio de signo entre f(a) y f(b)

f(0) = 3 \sin ^{3}(0)-1 - 4 \sin (0)\cos (0) = -1 f(\pi /2) = 3 \sin ^{3}(\pi /2)-1 - 4 \sin (\pi /2)\cos (\pi /2) = 3 (1)^3-1 - 4 (1)(0) = 2

con lo que se comprueba que al existir cambio de signo, debe existir una raíz en el intervalo.


c. Método de Falsa Posición


Desarrollo analítico con lápiz y papel

se trata de mostrar los pasos en al menos tres iteraciones del método, usando las siguientes expresiones:

f(t) = 3 \sin ^{3}(t)-1 - 4 \sin (t)\cos (t) c = b - f(b) \frac{a-b}{f(a)-f(b)}

[0, π/2]

iteración 1

a = 0 , b = \pi/2

tomando los datos al validar los puntos extremos

f(0) = -1 f(\pi /2) = 2 c = \pi/2 - 2 \frac{0-\pi/2}{-1-2} = \pi/6 f(\pi/6) = 3 \sin ^{3}(\pi/6)-1 - 4 \sin (\pi/6)\cos (\pi/6) = -2.3570

el signo de f(c) es el mismo que f(a), se ajusta el lado izquierdo

tramo = |c-a| = |\pi /6 - 0| = \pi/6 a = c = \pi/6 , b = \pi/2

iteración 2

a = \pi/6 , b = \pi/2 f(\pi/6) = -2.3570 f(\pi /2) = 2 c = \pi/2 - 2 \frac{\pi/6-\pi/2}{-1-2} = 1.0901 f(1.0901) = 3 \sin ^{3}(1.0901)-1 - 4 \sin (1.0901)\cos (1.0901) = -0.5486

el signo de f(c) es el mismo que f(a), se ajusta el lado izquierdo

tramo = |c-a| = | 1.0901-\pi/6| = 1.0722 a = c = 1.0901 , b = \pi/2

iteración 3

a = 1.0901 , b = \pi/2 f(1.0901) = -0.5486 f(\pi /2) = 2 c = \pi/2 - 2 \frac{-0.5486-\pi/2}{-0.5486-2} = 1.19358 f(1.19358) = 3 \sin ^{3}(1.19358)-1 - 4 \sin (1.19358)\cos (1.19358) = 0.0409

el signo de f(c) es el mismo que f(b), se ajusta el lado derecho

tramo = |b-c| = | \pi/2- 1.19358| = 0.3772 a = 1.0901 , b = 1.19358


Algoritmo con Python

Los parámetros aplicados en el algoritmo son los desarrollados en el planteamiento del problema e intervalo de búsqueda, con lo que se obtiene los siguientes resultados:

 raiz: 1.1864949811467547
error: 9.919955449055884e-05

las instrucciones en python son:

# 1Eva_IIT2009_T1 Movimiento de partícula en plano
import numpy as np
import matplotlib.pyplot as plt

#INGRESO
xt = lambda t: 3*(np.sin(t)**3)-1
yt = lambda t: 4*np.sin(t)*np.cos(t)
fx = lambda t: 3*(np.sin(t)**3)-1 - 4*np.sin(t)*np.cos(t) 

a = 0
b = np.pi/2
tolera = 0.001

# intervalo para gráfica
La = a
Lb = b
muestras = 21

# PROCEDIMIENTO
# Posicion Falsa
tramo = abs(b-a)
fa = fx(a)
fb = fx(b)
while not(tramo<=tolera):
    c = b - fb*(a-b)/(fa-fb)
    fc = fx(c)
    cambio = np.sign(fa)*np.sign(fc)
    if cambio>0:
        # actualiza en izquierda
        tramo = abs(c-a)
        a = c
        b = b
        fa = fc
    else:
        # actualiza en derecha
        tramo = abs(b-c)
        a = a
        b = c
        fb = fc

# para grafica
ti = np.linspace(La,Lb,muestras)
xi = xt(ti)
yi = yt(ti)
fi = fx(ti)

# SALIDA
print(' raiz:', c)
print('error:', tramo)
plt.plot(ti,xi, label='x(t)')
plt.plot(ti,yi, label='y(t)')
plt.plot(ti,fi, label='f(t)')
plt.plot(c,fx(c),'ro')
plt.axhline(0, color='green')
plt.axvline(c, color='magenta')
plt.legend()
plt.xlabel('t')
plt.title('Método de Falsa Posición')
plt.show()

s1Eva_IIT2010_T1 Aproximar con polinomio

Desarrollo Analítico

Ejemplo para Lápiz  y Papel


Tarea 01 Semana 01 Fecha: año/mes/dia
Apellidos Nombres
Referencia: 1Eva_IIT2010_T1 Aproximar con polinomio

Para el ejemplo, supondremos que x0=0

El polinomio de Taylor requerido es de grado 2

P_{n}(x) = f(x_0)+\frac{f'(x_0)}{1!} (x-x_0) + + \frac{f''(x_0)}{2!}(x-x_0)^2 +

función f(x) y sus derivadas:

f(x) = e^x \cos (x) +1
f'(x) = e^x \cos (x) - e^x \sin(x) f'(x) = e^x (\cos (x) - \sin(x))
f''(x) = e^x( \cos (x) - \sin(x))+ + e^x (-\sin(x) - \cos(x)) f''(x) = -2 e^x \sin(x))

Punto x0 = 0 (ejemplo), dentro del intervalo.

Observación: escriba las expresiones, reemplazando los valores, asi si en la lección o examen no tuvo tiempo para usar la calculadora, se puede evaluar si realizaba las operaciones con el punto de referencia y expresiones correctas.

f(0) = e^0 \cos (0) +1 = 2 f'(0) = e^0(\cos (0) - \sin(0)) = 1 f''(0) = -2 e^0 \sin(0)) = 0

Sustitución en fórmula de polinomio:

P_{2}(x) = 2+\frac{1}{1} (x-0) + + \frac{0}{2}(x-0)^2 + P_{2}(x) = 2+ x

Desarrollo con Python

Algoritmo desarrollado en clase, usado como taller, modificado para el problema planteado.

Observación: Se reordena el algoritmo para mantener ordenados y separados los bloques de ingreso, procedimiento y salida. Así los bloques pueden ser convertidos fácilmente a funciones algoritmicas def-return.

Observe que la variable n se interprete correctamente como «términos» o «grados» del polinomio de Taylor.

# Aproximación Polinomio de Taylor alrededor de x0
# f(x) en forma simbólica con sympy
import numpy as np
import sympy as sym
import matplotlib.pyplot as plt

# INGRESO --------------------
x = sym.Symbol('x')
fx = sym.exp(x)*sym.cos(x) + 1

x0 = 0
n  = 3 # grado de polinomio

# Intervalo para Gráfica
a = 0
b = np.pi
muestras = 21

# PROCEDIMIENTO  -------------
# construye polinomio Taylor
k = 0 # contador de términos
polinomio = 0
while (k <= n):
    derivada   = fx.diff(x,k)
    derivadax0 = derivada.subs(x,x0)
    divisor   = np.math.factorial(k)
    terminok  = (derivadax0/divisor)*(x-x0)**k
    polinomio = polinomio + terminok
    k = k + 1

# forma lambda para evaluación numérica
fxn = sym.lambdify(x,fx,'numpy')
pxn = sym.lambdify(x,polinomio,'numpy')

# evaluar en intervalo para gráfica
xi = np.linspace(a,b,muestras)
fxi = fxn(xi)
pxi = pxn(xi)

# SALIDA  --------------------
print('polinomio p(x)=')
print(polinomio)
print()
sym.pprint(polinomio)

# Gráfica
plt.plot(xi,fxi,label='f(x)')
plt.plot(xi,pxi,label='p(x)')
# franja de error
plt.fill_between(xi,pxi,fxi,color='yellow')
plt.xlabel('xi')
plt.axvline(x0,color='green', label='x0')
plt.axhline(0,color='grey')
plt.title('Polinomio Taylor: f(x) vs p(x)')
plt.legend()
plt.show()

Resultado del algoritmo

Revisar si el polinomio es concordante con lo realizado a lápiz y papel, de no ser así revisar el algoritmo o los pasos realizados en papel, deben ser iguales.
Comprobando que el algoritmo esté correcto y pueda ser usado en otros ejercicios.

 RESTART: D:\MATG1052Ejemplos\Taylot01Tarea01.py 
polinomio p(x)=
-x**3/3 + x + 2

   3        
  x         
- -- + x + 2
  3         

Resultados gráficos para x0=0

Continuar con el ejercicio con x0 = π y luego con el siguiente punto x0 = π/2.

Comparar resultados y presentar: Observaciones  y recomendaciones semejantes a las indicadas durante el desarrollo de la clase.

s3Eva_IIT2019_T4 completar polinomio de interpolación

Para visualizar la solución, se plantea graficar los polinomios que están completos S0(x) y S2(x).

S_0(x) = 1 + 1.1186x + 0.6938 x^3

 0.0 ≤ x ≤ 0.4

S_1(x) = 1.4918 + 1.4516(x-0.4) + c(x-0.4)^2 +d(x-0.4)^3

0.4 ≤ x ≤ 0.6

S_2(x) = 1.8221 + 1.8848(x-0.6) + +1.3336(x-0.6)^2 - 1.1113(x-0.6)^3

0.6 ≤ x ≤ 1.0

Como en trazadores cúbicos, lo que se usa es un polinomio por cada tramo muestreado para una curva contínua, etc. se tiene que los polinomios deben tener valores iguales en los puntos el eje x = 0.4 y 0.6

Por lo que se evalua con  los polinomios completos:

S_0(0.4) = 1 + 1.1186(0.4) + 0.6938 (0.4)^3 = 1.4918432 S_2(0.6) = 1.8221 + 1.8848(0.6-0.6) + +1.3336(0.6-0.6)^2-1.1113(0.6-0.6)^3=1.8221

Opción 1

Valores que se usan en los extremos del polinomio S1(x) para crear un sistema de dos ecuaciones y determinar los valores de c y d, completando el polinomio.

S_1(0.4) = 1.4918 + 1.4516(0.4-0.4) + c(0.4-0.4)^2 +d(0.4-0.4)^3

como los términos de c y d se hacen cero, hace falta una ecuación.

S_1(0.6) = 1.4918 + 1.4516(0.6-0.4) + + c(0.6-0.4)^2 +d(0.6-0.4)^3 = 1.8221 1.8221 = 1.4918 + 1.4516(0.2) + c(0.2)^2 +d(0.2)^3 0,03998 = 0.04c +0,008d

la otra ecuación se podría obtener usando la propiedad que las primeras derivadas de los polinomios deben ser iguales en los puntos x=0,4 y x= 0.6

S’0(0.2) =S’1(0.2)

Tarea: Desarrollar la siguiente ecuación y resolver


Opción 2

Si no recuerda la propiedad anterior, puede optar por usar otros conceptos para aproximar el resultado.

Si para el tramo en que se busca el polinomio se puede retroceder un tamaño de paso x = 0.2 y evualuar usando S0(0.2), se obtiene otropunto de referencia para crear un polinomio que pase por los mismos puntos.

S_0(0.2) = 1 + 1.1186*(0.2) + 0.6938 (0.2)^3

Se aplica lo mismo para un tamaño de paso más adelante de x = 0.6 es x = 0.8m se evalua S2(0.8) y se tienen suficientes puntos para usar cualquier método de interpolación y determinar el polinomio para el tramo faltante.

S_2(0.8) = 1.8221 + 1.8848(0.8-0.6) + +1.3336(0.8-0.6)^2 - 1.1113(0.8-0.6)^3
xi     = [0.        0.2        0.4      ]
S0(xi) = [1.        1.2292704  1.4918432]
xi     = [0.6       0.8        1.       ]
S2(xi) = [1.8221    2.2435136  2.7182728]

que permite hacer una tabla de puntos, y usando por ejemplo el método de interpolación de Lagrange  con x entre [0.2, 0.8] se obtiene otra forma del polinomio buscado:

p(x)=1.2292704 \frac{(x-0.4)(x-0.6)(x-0.8)}{(0.2-0.4)(0.2-0.6)(0.2-0.8)} + +1.4918432\frac{(x-0.2)(x-0.6)(x-0.8)}{0.4-0.2)(0.4-0.6)(0.4-0.8)}+ +1.8221\frac{(x-0.2)(x-0.4)(x-0.8)}{(0.6-0.2)(0.6-0.4)(0.6-0.8)} + + 2.2435136\frac{(x-0.2)(x-0.4)(x-0.6)}{(0.8-0.2)(0.8-0.4)(0.8-0.6)}

Tarea: continuar con el desarrollo


El literal b, requiere usar un metodo de búsqueda de raíces, para el cual se puede usar incluso bisección.

Tarea: continuar con el desarrollo

s3Eva_IIT2019_T3 Preparación de terreno en refineria

Se requiere usar el nivel inicial en la matriz, para restar del nivel requerido que es constante 220

Nivel inicio (m) 0 50 100 150 200
0 241 239 238 236 234
25 241 239 237 235 233
50 241 239 236 234 231
75 242 239 236 232 229
100 243 239 235 231 227

lo que genera la matriz de diferencias. El valor es positivo indica remoción, el valor negativo indica por rellenar.

Diferencia (m) 0 50 100 150 200
0 21 19 18 16 14
25 21 19 17 15 13
50 21 19 16 14 11
75 22 19 16 12 9
100 23 19 15 11 7

El volumen se puede calcular por un método en cada fila, y luego los resultados por columnas por otro método o el mismo.
Por ejemplo Simpson de 1/3

I= \frac{hx}{3}(f(x_0) +4f(x_1)+f(x_2))

con lo que se obtiene:

I_{fila}(0) = \frac{50}{3}(21 +4(19)+18) +\frac{50}{3}(18 +4(16)+14) =24750 I_{fila}(25) = = \frac{50}{3}(21 +4(19)+17) + \frac{50}{3}(17 +4(15)+13) = 23383,33 I_{fila}(50) = \frac{50}{3}(21 +4(19)+16) + \frac{50}{3}(16 +4(14)+11) = 22000 I_{fila}(75) = \frac{50}{3}(22 +4(19)+16) + \frac{50}{3}(16 +4(12)+5) =21850 I_{fila}(100) = \frac{50}{3}(23 +4(19)+15) + \frac{50}{3}(15 +4(11)+7) = 20483,33

y usando el otro eje, se completa el volumen usando dos veces simpson:

Volumen = \frac{h_y}{3}(f(x_0) +4f(x_1)+f(x_2)) Remover = \frac{25}{3}(24750 +4(23383,33)+22000) + + \frac{25}{3}(22000 +4(21850)+20483,33)=2251388,89

El signo lo trae desde la diferencia, y muestra el sentido del desnivel.

Se adjunta la gráfica de superficie en azul como referencia del signo,  respecto al nivel requerido en color verde.

Error de truncamiento

la cota del error de truncamiento se estima como O(h5)

error_{trunca} = -\frac{h^5}{90} f^{(4)}(z)

para un valor de z entre [a,b]

para cuantificar el valor, se puede usar la diferencia finita Δ4f, pues con la derivada sería muy laborioso.

s3Eva_IIT2019_T2 Diferenciación, valor en frontera

La ecuación de problema de valor de frontera, con h = 1/4 = 0.25:

y'' = -(x+1)y' + 2y + (1-x^2) e^{-x}

0 ≤ x ≤ 1
y(0) = -1
y(1) = 0

Se interpreta en la tabla como los valores que faltan por encontrar:

i 0 1 2 3 4
xi 0 1/4 1/2 3/4 1
yi -1 0

Por ejemplo, se usan las derivadas en diferencias finitas divididas centradas para segunda derivada y hacia adelante para primera derivada.

Semejante al procedimiento aplicado para métodos con EDO.

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

Se formula entonces la ecuación en forma discreta, usando solo los índices para los puntos yi:

\frac{y[i+1]-2y[i]+y[i-1]}{h^2} = -(x_i +1) \frac{y[i+1]-y[i]}{h} + 2y[i] + (1-x_i ^2) e^{-x_i}

se multiplica todo por h2

y[i+1]-2y[i]+y[i-1] = -h(x_i +1)(y[i+1]-y[i]) + + 2 h^2 y[i] + h^2 (1-x_i ^2) e^{-x_i}
y[i+1]-2y[i]+y[i-1] = -h(x_i +1)y[i+1] + + h(x_i +1)y[i] + 2 h^2 y[i] + h^2 (1-x_i ^2) e^{-x_i}
y[i+1](1 +h(x_i +1)) + y[i](-2 - h(x_i +1)-2 h^2)+y[i-1] = = h^2 (1-x_i ^2) e^{-x_i}

Ecuación que se aplica en cada uno de los puntos desconocidos con i =1,2,3


i = 1

y[2](1 +h(x_1 +1)) + y[1](-2 - h(x_1 +1)-2 h^2)+y[0] = = h^2 (1-x_1 ^2) e^{-x_1} y[2]\Big(1 +\frac{1}{4}\Big(\frac{1}{4} +1\Big)\Big) + y[1]\Big(-2 - \frac{1}{4}\Big(\frac{1}{4} +1\Big)-2 \Big(\frac{1}{4}\Big)^2\Big) -1 = = \Big(\frac{1}{4}\Big)^2 \Big(1-\Big(\frac{1}{4}\Big) ^2\Big) e^{-1/4} y[2]\Big(1 +\frac{5}{16} \Big) + y[1]\Big(-2 - \frac{5}{16}- \frac{2}{16}\Big) = = 1+ \frac{1}{16} \Big(1-\frac{1}{16}\Big) e^{-1/4} \frac{21}{16} y[2]- \frac{39}{16}y[1] = 1+ \frac{15}{16^2} e^{-1/4} - \frac{39}{16}y[1] + \frac{21}{16} y[2] = 1+ \frac{15}{16^2} e^{-1/4}

multiplicando ambos lados de la ecuacion por 16 y reordenando

- 39 y[1] + 21 y[2] = 16+ \frac{15}{16} e^{-1/4}

i = 2

y[3](1 +h(x_2 +1)) + y[2](-2 - h(x_2 +1)-2 h^2)+y[1] = = h^2 (1-x_2 ^2) e^{-x_2} y[3]\Big(1 +\frac{1}{4}\Big(\frac{1}{2} +1\Big)\Big) + y[2]\Big(-2 - \frac{1}{4}\Big(\frac{1}{2} +1\Big)-2 \Big(\frac{1}{4}\Big)^2\Big)+y[1] = = \Big(\frac{1}{4}\Big)^2 \Big(1-\Big(\frac{1}{2}\Big) ^2\Big) e^{-\frac{1}{2}} y[3]\Big(1 +\frac{3}{8}\Big) + y[2]\Big(-2 - \frac{1}{4}\Big(\frac{3}{2}\Big)- \frac{1}{8}\Big)+y[1] = = \Big(\frac{1}{16}\Big) \Big(1-\frac{1}{4}\Big) e^{-\frac{1}{2}} \frac{11}{8} y[3] - \frac{21}{8} y[2]+y[1] = \frac{1}{16} \Big(\frac{3}{4}\Big) e^{-\frac{1}{2}}

multiplicando ambos lados por 8 y reordenando,

8y[1] - 21 y[2] + 11 y[3] = \frac{3}{8} e^{-\frac{1}{2}}

i = 3

y[4](1 +h(x_3 +1)) + y[3](-2 - h(x_3 +1)-2 h^2)+y[2] = = h^2 (1-x_3 ^2) e^{-x_3} (0) \Big(1 +h\Big(x_3 +1\Big)\Big) + y[3]\Big(-2 - \frac{1}{4}\Big(\frac{3}{4} +1\Big)-2 \frac{1}{16}\Big)+y[2] = = \frac{1}{16}\Big (1-\Big(\frac{3}{4}\Big) ^2 \Big) e^{-3/4} y[3]\Big(-2 - \frac{7}{16} - \frac{2}{16})+y[2] = \frac{1}{16}\Big (1-\frac{9}{16}\Big) e^{-3/4}

multiplicando todo por 16 y reordenando:

16 y[2] - 41 y[3] = \frac{7}{16} e^{-3/4}

Con lo que se puede crear la forma matricial del sistema de ecuaciones Ax=B

\begin{bmatrix} -39 && 21 && 0\\8 && 21 && 11 \\ 0 && 16 && 41 \end{bmatrix}\begin{bmatrix} y[1]\\ y[2] \\y[3] \end{bmatrix} =\begin{bmatrix} 16+ \frac{15}{16} e^{-1/4} \\ \frac{3}{8} e^{-\frac{1}{2}} \\ \frac{7}{16} e^{-3/4} \end{bmatrix}

con lo que resolviendo la matriz se obtienen los valroes para y[1], y[2], y[3]

solución de matriz: 
[-0.59029143 -0.29958287 -0.12195088]

con lo que se completan los puntos de la tabla,

Solución de ecuación
x[i] =  [ 0.   0.25        0.5         0.75        1.  ]
y[i] =  [-1.  -0.59029143 -0.29958287 -0.12195088  0.  ]

con la siguiente gráfica:


Algoritmo para solución de matriz con Python

tarea: modificar para cambiar el valor del tamaño de paso.

# Problema de Frontera
import numpy as np
import matplotlib.pyplot as plt

# INGRESO

h = 1/4
y0 = -1
y1 = 0

xi = np.arange(0,1+h,h)
n = len(xi)
yi = np.zeros(n,dtype=float)
yi[0] = y0
yi[n-1] = y1
    
A = np.array([[-39.0, 21,  0],
              [  8.0, -21, 11],
              [  0.0, 16, -41]])
B = np.array([16+(15/16)*np.exp(-1/4),
              (3/8)*np.exp(-1/2),
              (7/16)*np.exp(-3/4)])

# PROCEDIMIENTO
x = np.linalg.solve(A,B)

for i in range(1,n-1,1):
    yi[i] = x[i-1]

# SALIDA
print('solución de matriz: ')
print(x)
print('Solución de ecuación')
print('x[i] = ',xi)
print('y[i] = ',yi)

# Grafica
plt.plot(xi,yi,'ro')
plt.plot(xi,yi)
plt.show()

s3Eva_IIT2019_T1 Lanzamiento de Cohete

A partir de la tabla del enunciado  se realiza la tabla de diferencias finitas.

i ti fi Δfi Δ2fi Δ3fi Δ4fi Δ5fi
1 0 0 32 -6 0 0 0
2 25 32 26 -6 0 0
3 50 58 20 -6 0
4 75 78 14 -6
5 100 92 8
6 125 100

Observando que a partir de la tercera diferencia finita  los valores son cero, por lo que se usa la fórmula general de diferencias finitas divididas hasta el polinomio de grado 2.

p_2 (x) = f_0 + \frac{\Delta f_0}{h} (x - x_0) + + \frac{\Delta^2 f_0}{2!h^2} (x - x_0)(x - x_1)

al sustituir los valores conocidos, se convierte en,

p_2 (t) =0 + \frac{32}{25} (t -0) + + \frac{-6}{2(25)^2} (t -0)(t - 25) =\frac{32}{25}t + \frac{-3}{(25)^2} (t^2 - 25t) =\frac{32}{25}t + \frac{-3}{(25)^2} t^2 - \frac{-3}{(25)^2}25t =\frac{7}{5}t - \frac{3}{625} t^2 y(t) =p_2 (t) =1.4 t - 0.0048 t^2

Con lo que se puede obtener la velocidad:

y'(t) = 1.4 - 0.0096 t

y luego la aceleración:

y''(t) = - 0.0096

Si el error es el próximo término del polinomio Δ3fi  entonces se estima en cero.

Tarea:  Evaluar la velocidad y aceleración para cada punto de la tabla

La gráfica del polinomio encontrado es:

Algoritmo en Python

El algoritmo realizado en Python entrega los siguientes resultados:

[[  i,  ti,  fi, df1, df2, df3, df4, df5,  df6]]
[[  1.   0.   0.  32.  -6.   0.   0.   0.   0.]
 [  2.  25.  32.  26.  -6.   0.   0.   0.   0.]
 [  3.  50.  58.  20.  -6.   0.   0.   0.   0.]
 [  4.  75.  78.  14.  -6.   0.   0.   0.   0.]
 [  5. 100.  92.   8.   0.   0.   0.   0.   0.]
 [  6. 125. 100.   0.   0.   0.   0.   0.   0.]]
polinomio:
-0.0048*t**2 + 1.4*t

las instrucciones en Python son:

# Diferencias finitas avanzadas para polinomio interpolación
# http://blog.espol.edu.ec/analisisnumerico/5-1-1-diferencias-finitas-avanzadas-polinomio/
# Referencia Rodriguez 6.6.4 Pdf.221
# Tarea: Verificar tamaño de vectores
#        considerar puntos no equidistantes en eje t
import numpy as np
import matplotlib.pyplot as plt
import sympy as sym

# INGRESO , Datos de prueba
ti = np.array([0.0, 25, 50, 75, 100, 125])
fi = np.array([0.0, 32, 58, 78, 92, 100])

# PROCEDIMIENTO
# Tabla de diferencias finitas
titulo = ['i','ti','fi']
n = len(ti)
# cambia a forma de columnas
i = np.arange(1,n+1,1)
i = np.transpose([i])
ti = np.transpose([ti])
fi = np.transpose([fi])
# Añade matriz de diferencias
dfinita = np.zeros(shape=(n,n),dtype=float)
tabla = np.concatenate((i,ti,fi,dfinita), axis=1)
# Sobre matriz de diferencias, por columnas
[n,m] = np.shape(tabla)
c = 3
diagonal = n-1
while (c<m):
    # Aumenta el título para cada columna
    titulo.append('df'+str(c-2))
    # calcula cada diferencia por fila
    f = 0
    while (f < diagonal):
        tabla[f,c] = tabla[f+1,c-1]-tabla[f,c-1]
        f = f+1
    
    diagonal = diagonal - 1
    c = c+1

# POLINOMIO con diferencias finitas
# caso: puntos en eje t equidistantes
dfinita = tabla[:,3:]
n = len(dfinita)
t = sym.Symbol('t')
h = ti[1,0]-ti[0,0]
polinomio = fi[0,0]
for c in range(1,n,1):
    denominador = np.math.factorial(c)*(h**c)
    factor = dfinita[0,c-1]/denominador
    termino=1
    for f  in range(0,c,1):
        termino = termino*(t-ti[f])
    polinomio = polinomio + termino*factor
# simplifica polinomio, multiplica los (t-ti)
polinomio = polinomio.expand()
# para evaluacion numérica
pt = sym.lambdify(t,polinomio)

# Puntos para la gráfica
a = np.min(ti)
b = np.max(ti)
muestras = 101
ti_p = np.linspace(a,b,muestras)
fi_p = pt(ti_p)

# SALIDA
print([titulo])
print(tabla)
print('polinomio:')
print(polinomio)

# Gráfica
plt.title('Interpolación polinómica')
plt.plot(ti,fi,'o', label = 'Puntos')
plt.plot(ti_p,fi_p, label = 'Polinomio')
plt.legend()
plt.show()

s2Eva_IIT2019_T4 Integrar con Cuadratura de Gauss

f(x) = x ln(x)

1 ≤x≤4

se requiere:

I = \int_1^4 x ln(x) dx

literal a. Usando el método de Cuadratura de Gauss con 2 términos

x_a = \frac{b+a}{2} + \frac{b-a}{2}x_0 = \frac{4+1}{2} + \frac{4-1}{2}\Big(\frac{-1}{\sqrt{3}} \Big)

xa =1.6339745962155612

x_b = \frac{b+a}{2} + \frac{b-a}{2}x_1 = \frac{4+1}{2} + \frac{4-1}{2}\Big(\frac{1}{\sqrt{3}} \Big)

xb =3.366025403784439

I \cong \frac{b-a}{2}(f(x_a) + f(x_b)) I \cong \frac{4-1}{2}(x_a ln(x_a) + x_b ln(x_b))

I = 7.33164251999249

literal b.  De la fórmula , despejar el valor del error<0.0001

\Big|\frac{(b-a)}{180}h^4 f^{(4)} (\xi)\Big| <0.0001; \xi \in[a,b] h^4 <0.0001\frac{180}{(4-1)}\frac{1}{f^{(4)} (\xi)} h^4 < 0.006\frac{1}{f^{(4)} (\xi)} h <\Big(0.006\frac{1}{f^{(4)} (\xi)}\Big)^{1/4}

obteniendo la 4ta derivada de la función:

f(x) = x ln(x) f'(x) = ln(x) + x\Big(\frac{1}{x} \Big) = ln(x) +1 f''(x) = \frac{1}{x} f'''(x) = -\frac{1}{x^2} f^{(4)}(x) = 2\frac{1}{x^3}

se tiene que:

h <\Big(0.006\frac{1}{f^{(4)} (\xi)}\Big)^{1/4} h <\Big(0.006\frac{1}{2\frac{1}{\xi^3}}\Big)^{1/4} h <\Big(0.003\xi^3\Big)^{1/4} h <(0.003)^{1/4}\xi^{3/4}

en el peor de los csis, se toma el valor menor de ξ =1

h <(0.003)^{1/4} h<0.2340347319320716

 

s2Eva_IIT2019_T3 EDP elíptica, placa en (1,1)

dada la ecuación del problema:

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

los tamaños de paso en ambos ejes son de igual valor, se simplifica la ecuación

\lambda= \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\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)

Iteraciones

que permite plantear las ecuaciones para cada punto en posición [i,j]


i=1, j=1

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

i=2, j=1

u[1,1]-4u[2,1]+u[3,1]+ +0.5 ln(0.5)+u[2,2]=0.15625 u[1,1]-4u[2,1]+u[3,1]+u[2,2]=0.15625-0.5 ln(0.5) u[1,1]-4u[2,1]+u[3,1]+u[2,2]=0.8493

i=3, j=1

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

Tarea: continuar con el ejercicio hasta plantear todo el sistema de ecuaciones.

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(
[0.125 - 2(0.25)ln(0.25),
 0.15625 - 0.5ln(0.5),
 0.20833 - 0.75 ln(0.75),
 ...])

B = [0.8181,
     0.8493,
     0.4240,
     ...]

Algoritmo con Python

Con valores para la matriz solución:

iteraciones:  15
error entre iteraciones:  6.772297286980838e-05
solución para u: 
[[0.        0.2789294 0.6081976 0.9793276 1.3862943]
 [0.2789294 0.6978116 1.1792239 1.7127402 2.2907268]
 [0.6081976 1.1792239 1.8252746 2.5338403 3.2958368]
 [0.9793276 1.7127402 2.5338403 3.4280053 4.3846703]
 [1.3862943 2.2907268 3.2958368 4.3846703 5.5451774]]
>>>

y algoritmo detallado:

# 2Eva_IIT2019_T2 EDO, problema de valor inicial
# 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 la mitad en intervalo x,
# mitad en intervalo y, para menos iteraciones
mx = int(n/2)
my = int(m/2)
promedio = (u[mx,0]+u[mx,m-1]+u[0,my]+u[n-1,my])/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
            fij = (dy**2)*fxy(xi[i],yj[j])
            u[i,j]=(u[i-1,j]+u[i+1,j]+u[i,j-1]+u[i,j+1]-fij)/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()