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

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

= 62623.3

x_{2} = 0.013 - \frac{-265.09}{62623.34} = 0.017233

error = |0.013 – 0.01723| = 0.004331

iteración 2

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

= 13.446

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

= 67897.5

x_{3} = 0.017233 - \frac{13.446}{67897.5} = 0.017031

error = |0.017233 – 0.017031| = 0.000198

cuyo valor 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%

E2_IIT2018_T4 Tasa Interes Hipoteca 01


Algoritmo en Python

El resultado con el algoritmo es:

método de Newton-Raphson
i ['xi', 'fi', 'dfi', 'xnuevo', 'tramo']
0 [ 1.30000e-02 -2.65091e+02  6.26233e+04  1.72331e-02  4.23311e-03]
1 [1.72331e-02 1.34468e+01 6.78975e+04 1.70351e-02 1.98045e-04]
2 [1.70351e-02 1.24433e-02 6.77706e+04 1.70349e-02 1.83609e-07]
raiz encontrada en:  0.017034880749732726
tasa anual:  0.20441856899679273

Instrucciones en Python

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

def newton_raphson(fx,dfx,xi, tolera, iteramax=100,
                   vertabla=False, precision=4):
    '''fx y dfx en forma numérica lambda
    xi es el punto inicial de búsqueda
    '''
    itera=0
    tramo = abs(2*tolera)
    if vertabla==True:
        print('método de Newton-Raphson')
        print('i', ['xi','fi','dfi', 'xnuevo', 'tramo'])
        np.set_printoptions(precision)
    while (tramo>=tolera):
        fi = fx(xi)
        dfi = dfx(xi)
        xnuevo = xi - fi/dfi
        tramo = abs(xnuevo-xi)
        if vertabla==True:
            print(itera,np.array([xi,fi,dfi,xnuevo,tramo]))
        xi = xnuevo
        itera = itera + 1

    if itera>=iteramax:
        xi = np.nan
        print('itera: ',itera,
              'No converge,se alcanzó el máximo de iteraciones')

    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)

x0 = 0.013 # 0.16/12
tolera = 0.0001

# PROCEDIMIENTO
raiz   = newton_raphson(fx, dfx, x0, tolera, vertabla=True, precision=5)
tanual = 12*raiz

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

# GRAFICA
import matplotlib.pyplot as plt
a = 0.01/12
b = 0.25/12
muestras = 21

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

plt.plot(tasa*12,fi, label="tasa anual")
plt.axhline(0, color='green')
plt.title('tasa anual de interes para Hipoteca')
plt.xlabel('tasa')
plt.ylabel('fx(tasa)')
plt.grid()
plt.legend()
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

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

cambio de signo a la izquierda

a = 0.01, b=0.1325

iteración 3

c = \frac{a+b}{2} = \frac{0.01+0.1325}{2} = 0.07125 fx(0.07125) = 70000 - 1200\Big(\frac{1-(1+(0.07125))^{-300}}{0.07125} \Big) = 53157.89 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 estImado 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

t [s] 0 2 4 6 8
v(t) [m/s] 0.0 16.40 27.77 35.64 41.10

https://www.dreamstime.com/stock-photo-skydiving-formation-group-people-image62015024No se especifica el método a seguir, por lo que se puede seleccionar el de mayor preferencia.

Por ejemplo. usando el método de Lagrange, con los puntos primero, medio y último, para cubrir todo el intervalo:

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.74t

2_IIT2018_T1 Interpola Paracaidista 01


Literal b

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

Una de las formas, que es independiente si se resolvió el literal a, es usar los datos proporcionados en la tabla el ejercicio:

t [s] 0 2 4 6 8
v(t) [m/s] 0.0 16.40 27.77 35.64 41.10

Se podría usar el método de Simpson de 1/3, puesto que los tamaños de paso en t son equidistantes se puede aplicar: h=2-0=2

\int_0^8 v(t)dt = \frac{2}{3}\Big( 0+ 4(16.40)+27.77\Big) + \frac{2}{3}\Big( 27.77+ 4(35.64)+41.10\Big) =203.2

con error del orden de O(h5) que al considerar h=2 no permite hacer una buena estimación del error. Sin embargo la respuesta es bastante cercana si se usa el método el trapecio con el algoritmo:

    valores de fi:  [ 0.   27.77 41.1 ]
divisores en L(i):  [ 32. -16.  32.]

Polinomio de Lagrange, expresiones
-1.735625*x*(x - 8.0) + 1.284375*x*(x - 4.0)

Polinomio de Lagrange: 
-0.45125*x**2 + 8.7475*x
Método del trapecio
distancia recorrida:  193.28
>>> 

El error entre los métodos es |203.2-193.28|= 9.92

Revisar el resultado usando un método con mayor precisión que el trapecio.


Algoritmo 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 sympy as sym
import matplotlib.pyplot as plt

# Literal a)
# Interpolacion de Lagrange
# divisoresL solo para mostrar valores

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

cuales = [0,2,4]

# PROCEDIMIENTO
xi = np.array(t,dtype=float)
fi = np.array(v,dtype=float)

xi = xi[cuales]
fi = fi[cuales]

# Polinomio de Lagrange
n = len(xi)
x = sym.Symbol('x')
polinomio = 0
divisorL = np.zeros(n, dtype = float)
for i in range(0,n,1):
    
    # Termino de Lagrange
    numerador = 1
    denominador = 1
    for j  in range(0,n,1):
        if (j!=i):
            numerador = numerador*(x-xi[j])
            denominador = denominador*(xi[i]-xi[j])
    terminoLi = numerador/denominador

    polinomio = polinomio + terminoLi*fi[i]
    divisorL[i] = denominador

# simplifica el polinomio
polisimple = polinomio.expand()

# para evaluación numérica
px = sym.lambdify(x,polisimple)

# Puntos para la gráfica
muestras = 51
a = np.min(xi)
b = np.max(xi)
pxi = np.linspace(a,b,muestras)
pfi = px(pxi)

# SALIDA
print('    valores de fi: ',fi)
print('divisores en L(i): ',divisorL)
print()
print('Polinomio de Lagrange, expresiones')
print(polinomio)
print()
print('Polinomio de Lagrange: ')
print(polisimple)

# Gráfica
plt.plot(t,v,'o', label = 'Puntos')
plt.plot(xi,fi,'o', label = 'Puntos en polinomio')
plt.plot(pxi,pfi, label = 'Polinomio')
plt.legend()
plt.xlabel('xi')
plt.ylabel('fi')
plt.title('Interpolación Lagrange')
plt.grid()
plt.show()

# Literal b
# INGRESO
# El ingreso es el polinomio en forma lambda
# se mantienen las muestras

# intervalo de integración
# a, b seleccionados para la gráfica anterior
tramos = muestras -1

# PROCEDIMIENTO
def integratrapecio_fi(xi,fi):
    ''' sobre muestras de fi para cada xi
        integral con método de trapecio
    '''
    n = len(xi)
    suma = 0
    for i in range(0,n-1,1):
        dx = xi[i+1]-xi[i]
        untrapecio = dx*(fi[i+1]+fi[i])/2
        suma = suma + untrapecio
    return(suma)


tramos = muestras-1
# PROCEDIMIENTO
distancia = integratrapecio_fi(xi,fi)

# SALIDA
print('Método del trapecio')
print('distancia recorrida: ', distancia)

s1Eva_IIT2018_T3 Interpolar con sistema de ecuaciones

Ejercicio: 1Eva_IIT2018_T3 Interpolar con sistema de ecuaciones

Los datos del ejercicio proporcionados son:

i 0 1 2 3 4 5
x 1.0 1.1 1.3 1.5 1.9 2.1
y(x) 1.84 1.90 2.10 2.28 2.91 3.28

Literal a

El tema es semejante al tema 1, cambiando el método de interpolación.
Se usan los puntos de las posiciones 0, 3 y 5.

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}

Matriz Aumentada

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

Pivoteo parcial por filas

Para el primer pivote no se requieren cambio de filas.
para el segundo pivote de la diagonal se deben intercambiar la fila segunda con la tercera

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

Se aplica eliminación hacia adelante:

fila = 0, columna=0  pivote = AB[0,0]=1

factor entre las filas es 1/1=1.

\begin{bmatrix}1 & 1 & 1 & 1.84 \\ 1-1 & 2.1-1 & 4.41 -1 &3.28 -1.84 \\ 1-1 & 1.5 -1 & 2.25 -1 & 2.28 - 1.84 \end{bmatrix} \begin{bmatrix} 1 & 1 & 1 & 1.84 \\ 0 & 1.1 & 3.41 &1.44 \\ 0 & 0.5 & 1.25 & 0.44 \end{bmatrix}

fila =1,  columna=1, pivote=AB[1,1] =1.1

factor entre filas es 0.5/1.1 = 1/2.2

\begin{bmatrix} 1 & 1 & 1 & 1.84 \\ 0 & 1.1 & 3.41 &1.44 \\ 0 & 0.5 -\frac{0.5}{1.1}(1.1)& 1.25 -\frac{0.5}{1.1}(3.41)& 0.44-\frac{0.5}{1.1}(1.44) \end{bmatrix} \begin{bmatrix} 1 & 1 & 1 & 1.84 \\ 0 & 1.1 & 3.41 &1.44 \\ 0 & 0 & -0.3 & -0.214545 \end{bmatrix}

aplicando sustitución hacia atrás

b2 = -0.21/(-0.3) = 0.71515 b1= \frac{1.44-3.41 b_2}{1.1} = \frac{1.44-3.41( 0.71515)}{1.1}=-0.9078 b3= \frac{1.84-b_1-b_2}{1} = \frac{1.84-(-0.9078)-(0.71515)}{1} =2.0327

con lo que el polinomio buscado es:

p_2(x) = 2.0327 -0.9078 x + 0.71515 x^2

y se obtiene el resultado de la interpolación.

E2_IIT2018_T3 Interpola Sistema Ecuaciones 01Observación: En la gráfica se muestra que el polinomio pasa por los puntos seleccionados de la tabla. En los otros puntos hay un error que se puede calcular como la resta del punto y su valor con p(x). Queda como tarea.

Usando el algoritmo del polinomio de interpolación con la matriz de Vandermonde se obtiene:

Matriz Vandermonde: 
[[1.   1.   1.  ]
 [2.25 1.5  1.  ]
 [4.41 2.1  1.  ]]
los coeficientes del polinomio: 
[ 0.71515152 -0.90787879  2.03272727]
Polinomio de interpolación: 
0.715151515151516*x**2 - 0.907878787878792*x + 2.03272727272728

 formato pprint
                   2                                         
0.715151515151516*x  - 0.907878787878792*x + 2.03272727272728
suma de columnas:  [3.   4.75 7.51]
norma D:  7.51
numero de condicion:  97.03737354737122
solucion: 
[ 0.71515152 -0.90787879  2.03272727]
>>> 


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 condición:  97.03737354737129
solución: 
[ 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

# 1Eva_IIT2018_T3 Interpolar con sistema de ecuaciones
# El polinomio de interpolación
import numpy as np
import sympy as sym
import matplotlib.pyplot as plt

# INGRESO
xj = [1.0,  1.1,  1.3,  1.5,  1.9,  2.1 ]
yj = [1.84, 1.90, 2.10, 2.28, 2.91, 3.28]
cuales = [0, 3, 5]

# muestras = tramos+1
muestras = 51

# PROCEDIMIENTO
# Convierte a arreglos numpy 
xi = np.array(xj,dtype=float)
fi = np.array(yj,dtype=float)

xi = xi[cuales]
fi  = fi[cuales]
B = fi

n = len(xi)

# Matriz Vandermonde D
D = np.zeros(shape=(n,n),dtype =float)
for i in range(0,n,1):
    for j in range(0,n,1):
        potencia = (n-1)-j # Derecha a izquierda
        D[i,j] = xi[i]**potencia

# Aplicar métodos Unidad03. Tarea
# Resuelve sistema de ecuaciones A.X=B
coeficiente = np.linalg.solve(D,B)

# Polinomio en forma simbólica
x = sym.Symbol('x')
polinomio = 0
for i in range(0,n,1):
    potencia = (n-1)-i   # Derecha a izquierda
    termino = coeficiente[i]*(x**potencia)
    polinomio = polinomio + termino

# Polinomio a forma Lambda
# para evaluación con vectores de datos xin
px = sym.lambdify(x,polinomio)

# Para graficar el polinomio en [a,b]
a = np.min(xi)
b = np.max(xi)
xin = np.linspace(a,b,muestras)
yin = px(xin)

# Usando evaluación simbólica
##yin = np.zeros(muestras,dtype=float)
##for j in range(0,muestras,1):
##    yin[j] = polinomio.subs(x,xin[j])
    
# SALIDA
print('Matriz Vandermonde: ')
print(D)
print('los coeficientes del polinomio: ')
print(coeficiente)
print('Polinomio de interpolación: ')
print(polinomio)
print('\n formato pprint')
sym.pprint(polinomio)

# Grafica
plt.plot(xj,yj,'o', label='Puntos')
plt.plot(xi,B,'o', label='[xi,fi]')
plt.plot(xin,yin, label='p(x)')
plt.xlabel('xi')
plt.ylabel('fi')
plt.legend()
plt.title(polinomio)
plt.show()


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

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

solucion = np.linalg.solve(D,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 intervalo 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 usar un método de búsqueda de raíces, se requiere encontrar el valor cuando f(x) = d’ = 0.

Un método como el de Newton Raphson requiere también f'(x) = d''

f(x) = \frac{x + (e^x - 1)e^x - 1}{\sqrt{(x - 1)^2 + (e^x - 1)^2}} f'(x)= \frac{(e^x - 1)e^x + e^{2x} + 1 - \frac{(x + (e^x - 1)e^x - 1)^2}{(x - 1)^2 + (e^x - 1)^2}} {\sqrt{(x - 1)^2 + (e^x - 1)^2}}

expresiones obtenidas usando Sympy

f(x) :
(x + (exp(x) - 1)*exp(x) - 1)/sqrt((x - 1)**2 + (exp(x) - 1)**2)
f'(x) :
((exp(x) - 1)*exp(x) + exp(2*x) + 1 - (x + (exp(x) - 1)*exp(x) - 1)**2/((x - 1)**2 + (exp(x) - 1)**2))/sqrt((x - 1)**2 + (exp(x) - 1)**2)

f(x) :
       / x    \  x        
   x + \e  - 1/*e  - 1    
--------------------------
    ______________________
   /                    2 
  /         2   / x    \  
\/   (x - 1)  + \e  - 1/  
f'(x) :
                                              2
                         /    / x    \  x    \ 
/ x    \  x    2*x       \x + \e  - 1/*e  - 1/ 
\e  - 1/*e  + e    + 1 - ----------------------
                                             2 
                                 2   / x    \  
                          (x - 1)  + \e  - 1/  
-----------------------------------------------
               ______________________          
              /                    2           
             /         2   / x    \            
           \/   (x - 1)  + \e  - 1/            


lo que permite observar la raíz de f(x) en una gráfica:
distancia mínima f(x)
con las siguientes instrucciones:

# Eva_IIT2018_T2 Distancia mínima a un punto
import numpy as np
import matplotlib.pyplot as plt
import sympy as sym

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

a = 0
b = 1
muestras = 21

# PROCEDIMIENTO
dfx = sym.diff(fx,x,1)
d2fx = sym.diff(fx,x,2)

f = sym.lambdify(x,dfx)
xi = np.linspace(a,b,muestras)
fi = f(xi)


# SALIDA
print('f(x) :')
print(dfx)
print("f'(x) :")
print(d2fx)
print()
print('f(x) :')
sym.pprint(dfx)
print("f'(x) :")
sym.pprint(d2fx)

# GRAFICA
plt.plot(xi,fi, label='f(x)')
plt.axhline(0)
plt.xlabel('x')
plt.ylabel('f(x)')
plt.legend()
plt.grid()
plt.show()

Usando el método de la bisección para el intervalo dado, se tiene:

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

itera = 0 , a = 0, b=1

c= \frac{0+1}{2} = 0.5 f(0) = \frac{0 + (e^0 - 1)e^0 - 1}{\sqrt{(0 - 1)^2 + (e^0 - 1)^2}} = -1 f(1) = \frac{1 + (e^1 - 1)e^1 - 1}{\sqrt{(1 - 1)^2 + (e^1 - 1)^2}} 2.7183 f(0.5) = \frac{(0.5) + (e^(0.5) - 1)e^(0.5) - 1}{\sqrt{((0.5) - 1)^2 + (e^(0.5) - 1)^2}} = 0.6954

cambio de signo a la izquierda,

a= 0, b=c=0.5

tramo = |0.5-0| = 0.5

itera = 1

c= \frac{0+0.5}{2} = 0.25 f(0.25) = \frac{(0.25) + (e^(0.25) - 1)e^(0.25) - 1}{\sqrt{((0.25) - 1)^2 + (e^(0.25) - 1)^2}} = -0.4804

cambio de signo a la derecha,

a=c= 0.25, b=0.5

itera = 2

c= \frac{0.25+0.5}{2} = 0.375 f(0.375) = \frac{(0.375) + (e^(0.375) - 1)e^(0.375) - 1}{\sqrt{((0.375) - 1)^2 + (e^(0.375) - 1)^2}} = 0.0479

cambio de signo a la izquierda,

a= 0.25, b=c=0.375

se continúan las iteraciones con el algoritmo, para encontrar la raíz en 0.364:

método de Bisección
i ['a', 'c', 'b'] ['f(a)', 'f(c)', 'f(b)']
   tramo
0 [0, 0.5, 1] [-1.      0.6954  2.7183]
   0.5
1 [0, 0.25, 0.5] [-1.     -0.4804  0.6954]
   0.25
2 [0.25, 0.375, 0.5] [-0.4804  0.0479  0.6954]
   0.125
3 [0.25, 0.3125, 0.375] [-0.4804 -0.2388  0.0479]
   0.0625
4 [0.3125, 0.34375, 0.375] [-0.2388 -0.1004  0.0479]
   0.03125
5 [0.34375, 0.359375, 0.375] [-0.1004 -0.0274  0.0479]
   0.015625
6 [0.359375, 0.3671875, 0.375] [-0.0274  0.01    0.0479]
   0.0078125
7 [0.359375, 0.36328125, 0.3671875] [-0.0274 -0.0088  0.01  ]
   0.00390625
8 [0.36328125, 0.365234375, 0.3671875] [-0.0088  0.0006  0.01  ]
   0.001953125
9 [0.36328125, 0.3642578125, 0.365234375] [-0.0088 -0.0041  0.0006]
   0.0009765625
raíz en:  0.3642578125

Al algoritmo anterior se complementa con las instrucciones de la función para la bisección.

# Eva_IIT2018_T2 Distancia mínima a un punto
import numpy as np
import matplotlib.pyplot as plt
import sympy as sym

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

a = 0
b = 1
muestras = 21

# PROCEDIMIENTO
dfx = sym.diff(fx,x,1)
d2fx = sym.diff(fx,x,2)

f = sym.lambdify(x,dfx)
xi = np.linspace(a,b,muestras)
fi = f(xi)


# SALIDA
print('f(x) :')
print(dfx)
print("f'(x) :")
print(d2fx)
print()
print('f(x) :')
sym.pprint(dfx)
print("f'(x) :")
sym.pprint(d2fx)

# GRAFICA
plt.plot(xi,fi, label='f(x)')
plt.axhline(0)
plt.xlabel('x')
plt.ylabel('f(x)')
plt.legend()
plt.grid()
plt.show()

# Algoritmo de Bisección
# [a,b] se escogen de la gráfica de la función
# error = tolera
import numpy as np

def biseccion(fx,a,b,tolera,iteramax = 20, vertabla=False, precision=4):
    '''
    Algoritmo de Bisección
    Los valores de [a,b] son seleccionados
    desde la gráfica de la función
    error = tolera
    '''
    fa = fx(a)
    fb = fx(b)
    tramo = np.abs(b-a)
    itera = 0
    cambia = np.sign(fa)*np.sign(fb)
    if cambia<0: # existe cambio de signo f(a) vs f(b)
        if vertabla==True:
            print('método de Bisección')
            print('i', ['a','c','b'],[ 'f(a)', 'f(c)','f(b)'])
            print('  ','tramo')
            np.set_printoptions(precision)
            
        while (tramo>=tolera and itera<=iteramax):
            c = (a+b)/2
            fc = fx(c)
            cambia = np.sign(fa)*np.sign(fc)
            if vertabla==True:
                print(itera,[a,c,b],np.array([fa,fc,fb]))
            if (cambia<0):
                b = c
                fb = fc
            else:
                a = c
                fa = fc
            tramo = np.abs(b-a)
            if vertabla==True:
                print('  ',tramo)
            itera = itera + 1
        respuesta = c
        # Valida respuesta
        if (itera>=iteramax):
            respuesta = np.nan

    else: 
        print(' No existe cambio de signo entre f(a) y f(b)')
        print(' f(a) =',fa,',  f(b) =',fb) 
        respuesta=np.nan
    return(respuesta)

# INGRESO
tolera = 0.001

# PROCEDIMIENTO
respuesta = biseccion(f,a,b,tolera,vertabla=True)
# SALIDA
print('raíz en: ', respuesta)

.

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