s3Eva_IT2019_T1 Ecuaciones simultáneas

Ejercicio: 3Eva_IT2019_T1 Ecuaciones simultáneas

Para plantear la intersección de las ecuaciones se pueden simplificar como:

y_1 = -x^2 +x + 0.75 y+5xy=x^3 y(1+5x)=x^3 y_2=\frac{x^3}{1+5x}

Quedando dos ecuaciones simplificadas:

y_1 = -x^2 +x + 0.75 y_2 = \frac{x^3}{1+5x}

cuyas gráficas son:

dónde se puede observar la intersección alrededor de 1.3

Restando ambas ecuaciones, se tiene que encontrar el valor de x para que el resultado sea cero.

y_1(x)-y_2(x)= -x^2 +x + 0.75 -\frac{x^3}{1+5x} f(x) = -x^2 +x + 0.75 -\frac{x^3}{1+5x} = 0

Para encontrarla derivada se procesa la expresión:

(1+5x)(-x^2 +x + 0.75) -x^3 = 0(1+5x) -6x^3+4x^2+4.75x+0.75 = 0 f'(x)= -18x^2 +8x + 4.75

Se usa el punto inicial x0=1 definido en el enunciado y se realizan las iteraciones siguiendo el algoritmo.

Se tiene que la raiz es:

raiz en:  1.3310736382369661
 [  xi, 	 xnuevo,	 fxi,	 dfxi, 	 tramo]
[[ 1.000e+00  1.111e+00  5.833e-01 -5.250e+00  1.111e-01]
 [ 1.111e+00  1.160e+00  4.173e-01 -8.583e+00  4.862e-02]
 [ 1.160e+00  1.193e+00  3.353e-01 -1.018e+01  3.293e-02]
 [ 1.193e+00  1.217e+00  2.766e-01 -1.131e+01  2.445e-02]
 [ 1.217e+00  1.236e+00  2.313e-01 -1.218e+01  1.899e-02]
 [ 1.236e+00  1.251e+00  1.951e-01 -1.286e+01  1.517e-02]
....
]

Algoritmo en Python

# 3Eva I T 2019 ecuaciones simultaneas
import numpy as np
import matplotlib.pyplot as plt

def newton_raphson(funcionx, fxderiva, xi, tolera):
    '''
    funciónx y fxderiva son de forma numérica
    xi es el punto inicial de búsqueda
    '''
    tabla = []
    tramo = abs(2*tolera)
    while (tramo>=tolera):
        fxi = funcionx(xi)
        dfxi = fxderiva(xi)
        xnuevo = xi - fxi/dfxi
        tramo = abs(xnuevo-xi)
        tabla.append([xi,xnuevo,fxi,dfxi,tramo])
        xi = xnuevo
    return(xi,tabla)

# INGRESO
y1 = lambda x: -x**2 +x +0.75
y2 = lambda x: (x**3)/(1+5*x)
a = 0.5
b = 1.5
muestras = 20

f = lambda x: -x**2+x+0.75-x**3/(1+5*x)
df = lambda x: -18*(x**2)+8*x +4.75
tolera = 1e-4
x0 = 1

# PROCEDIMIENTO
# datos para la gráfica
xi = np.linspace(a,b,muestras)
yi1 = y1(xi)
yi2 = y2(xi)
fi = f(xi)
# determina raiz
raiz, tabla = newton_raphson(f, df, x0, tolera)
tabla = np.array(tabla)

# SALIDA
np.set_printoptions(precision=3)
print('raiz en: ',raiz)
print(' [  xi, \t xnuevo,\t fxi,\t dfxi, \t tramo]')
print(tabla)

# Gráfica
plt.plot(xi,yi1, label ='yi1')
plt.plot(xi,yi2, label ='yi2')
plt.plot(xi,fi, label ='fi=yi1-yi2')
plt.axvline(raiz,linestyle='dashed')
plt.axhline(0)
plt.xlabel('x')
plt.legend()
plt.title('ecuaciones simultáneas')
plt.show()

s2Eva_IT2019_T3 EDP Elíptica Placa 6×5

Ejercicio: 2Eva_IT2019_T3 EDP Elíptica Placa 6×5

La ecuación se discretiza con diferencias divididas centradas

\frac{\partial^2 u}{\partial x^2} (x,y)+\frac{\partial ^2 u}{\partial y^2 } (x,y) = -\frac{q}{K} \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{q}{K}

Para procesar los datos se toma como referencia la malla

se agrupan los términos conocidos de las diferencias divididas

\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}=-\frac{q}{K}(\Delta y)^2

Considerando que hx =2, hy= 2.5 son diferentes, se mantiene el valor de λ para el desarrollo.

\lambda = \frac{(\Delta y)^2}{(\Delta x)^2} = \frac{(2.5)^2}{(2)^2} = 1.5625 \lambda u_{i+1,j}-2\lambda u_{i,j} +\lambda u_{i-1,j} + u_{i,j+1}-2u_{i,j}+u_{i,j-1}=-\frac{q}{K}(\Delta y)^2

Se agrupan términos de u que son iguales

\lambda u_{i+1,j}-2(1+\lambda) u_{i,j} +\lambda u_{i-1,j} + u_{i,j+1}+u_{i,j-1}=-\frac{q}{K}(\Delta y)^2

con lo que se puede generar el sistema de ecuaciones a resolver para los nodos de la malla.

i = 1 , j=1

1.5625 u_{2,1}-2(1+ 1.5625) u_{1,1} + 1.5625 u_{0,1} + +u_{1,2}+u_{1,0}=-\frac{1.5}{1.04}(2.5)^2 1.5625 u_{2,1}-5.125u_{1,1} + 1.5625 y_1(5-y_1) + + 0+x_i(6-x_i)=-9.0144 1.5625 u_{2,1}-5.125u_{1,1} + 1.5625[ 1(5-1)]+ + 0+2(6-2)=-9.0144 1.5625 u_{2,1}-5.125 u_{1,1} = =-9.0144 - 1.5625 [1(5-1)] - 0-2(6-2)

 

1.5625 u_{2,1}-5.125 u_{1,1} =-23.2644

i=2, j=1

1.5625 u_{3,1}-2(1+1.5625) u_{2,1} +1.5625 u_{1,1} + +u_{2,2}+u_{2,0}=-9.0144 1.5625 (0)-5.125 u_{2,1} +1.5625 u_{1,1} + 0 +x_2(6-x_2)=-9.0144 -5.125 u_{2,1} +1.5625 u_{1,1} + 0 +4(6-4)=-9.0144

 

-5.125 u_{2,1} +1.5625 u_{1,1} =-17.0144

las ecuaciones a resolver son:

1.5625 u_{2,1}-5.125 u_{1,1} =-23.2644 -5.125 u_{2,1} +1.5625 u_{1,1} =-17.0144

cuya solución es:

u_{2,1} = 5.1858 u_{1,1} =6.1204

Resuelto usando:

import numpy as np
A = np.array([[1.5625,-5.125],
             [-5.125, 1.5625]])
B = np.array([-23.2644,-17.0144])
x = np.linalg.solve(A,B)
print(x)

s2Eva_IT2019_T2 Péndulo vertical

Ejercicio: 2Eva_IT2019_T2 Péndulo vertical

Para resolver la ecuación usando Runge-Kutta se estandariza la ecuación a la forma:

\frac{d^2\theta }{dt^2}+\frac{g}{L}\sin (\theta)=0 \frac{d^2\theta }{dt^2}= -\frac{g}{L}\sin (\theta)

Se simplifica su forma a:

\frac{d\theta}{dt}=z = f_t(t,\theta,z) \frac{d^2\theta }{dt^2}= z' =-\frac{g}{L}\sin (\theta) = g_t(t,\theta,z)

y se usan los valores iniciales para iniciar la tabla:

\theta(0) = \frac{\pi}{6} \theta '(0) = 0
ti θ(ti) θ'(ti)=z
0 π/6 = 0.5235 0
0.2 0.3602 -1.6333
0.4 -0.0815 -2.2639

para las iteraciones se usan los valores (-g/L) = (-9.8/0.6)

Iteración 1:  ti = 0 ; yi = π/6 ; zi = 0

K1y = h * ft(ti,yi,zi) 
    = 0.2*(0) = 0
K1z = h * gt(ti,yi,zi) 
    = 0.2*(-9.8/0.6)*sin(π/6) = -1.6333
        
K2y = h * ft(ti+h, yi + K1y, zi + K1z)
    = 0.2*(0-1.6333)= -0.3266
K2z = h * gt(ti+h, yi + K1y, zi + K1z)
    = 0.2*(-9.8/0.6)*sin(π/6+0) = -1.6333

yi = yi + (K1y+K2y)/2 
   = π/6+ (0+(-0.3266))/2 = 0.3602
zi = zi + (K1z+K2z)/2 
   = 0+(-1.6333-1.6333)/2 = -1.6333
ti = ti + h = 0 + 0.2 = 0.2

estimado[i] = [0.2,0.3602,-1.6333]

Iteración 2: ti = 0.2 ; yi = 0.3602 ; zi = -1.6333

K1y = 0.2*( -1.6333) = -0.3266
K1z = 0.2*(-9.8/0.6)*sin(0.3602) = -1.1515
        
K2y = 0.2*(-1.6333-0.3266)= -0.5569
K2z = 0.2*(-9.8/0.6)*sin(0.3602-0.3266) = -0.1097

yi = 0.3602 + ( -0.3266 + (-0.3919))/2 = -0.0815
zi = -1.6333+(-1.151-0.1097)/2 = -2.2639
ti = ti + h = 0.2 + 0.2 = 0.4

estimado[i] = [0.4,-0.0815,-2.2639]

Se continúan con las iteraciones, hasta completar la tabla.

Tarea: realizar la Iteración 3

Usando el algoritmo RungeKutta_fg se obtienen los valores y luego la gráfica

 [ t, 		 y, 	 dyi/dti=z]
[[ 0.          0.52359878  0.        ]
 [ 0.2         0.36026544 -1.63333333]
 [ 0.4        -0.08155862 -2.263988  ]
 [ 0.6        -0.50774327 -1.2990876 ]
 [ 0.8        -0.60873334  0.62920692]
 [ 1.         -0.29609456  2.32161986]]

si se mejora la resolución disminuyendo h = 0.05 y muestras = 20 para cubrir el dominio [0,1] se obtiene el siguiente resultado:

Tarea: Para el literal b, se debe considerar que los errores se calculan con

Runge-Kutta 2do Orden tiene error de truncamiento O(h3)
Runge-Kutta 4do Orden tiene error de truncamiento O(h5)


Algoritmo en Python

# 3Eva_IT2019_T2 Péndulo vertical
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
g = 9.8
L  = 0.6
ft = lambda t,y,z: z
gt = lambda t,y,z: (-g/L)*np.sin(y)

t0 = 0
y0 = np.pi/6
z0 = 0
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()

Otra solución propuesta puede seguir el problema semejante:

s2Eva_IT2010_T2 Movimiento angular

s2Eva_IT2019_T1 Esfuerzo en pulso cardiaco

Ejercicio: 2Eva_IT2019_T1 Esfuerzo en pulso cardiaco

Para resolver el ejercicio, la función a integrar es el cuadrado de los valores. Para minimizar los errores se usarán TODOS los puntos muestreados, aplicando los métodos adecuados.
Con aproximación de Simpson se requiere que los tamaños de paso sean iguales en cada segmento.
Por lo que primero se revisa el tamaño de paso entre lecturas.

tamaño de paso h:
[0.04 0.04 0.02 0.01 0.01 0.01 0.03 0.04 0.03 0.02 0.  ]
tiempos:
[0.   0.04 0.08 0.1  0.11 0.12 0.13 0.16 0.2  0.23 0.25]
ft:
[ 10.  18.   7.  -8. 110. -25.   9.   8.  25.   9.   9.]

Observando los tamaños de paso se tiene que:
– entre dos tamaños de paso iguales se usa Simpson de 1/3
– entre tres tamaños de paso iguales se usa Simpson de 3/8
– para tamaños de paso variables se usa trapecio.

Se procede a obtener el valor del integral,

Intervalo [0,0.8], h = 0.04

I_{S13} = \frac{0.04}{3}(10^2+4(18^2)+7^2)

Intervalo [0.08,0.1], h = 0.02

I_{Tr1} = \frac{0.02}{2}(7^2+(-8)^2)

Intervalo [0.1,0.13], h = 0.01

I_{S38} = \frac{3}{8}(0.01)((-8)^2+3(110)^2+3(-25)^2+9^2)

Intervalo [0.13,0.25], h = variable

I_{Tr2} = \frac{0.03}{2}(9^2+8^2) I_{Tr3} = \frac{0.04}{2}(8^2+25^2) I_{Tr4} = \frac{0.03}{2}(25^2+9^2) I_{Tr5} = \frac{0.02}{2}(9^2+9^2)

El integral es la suma de los valores parciales, y con el resultado se obtiene el valor Xrms requerido.

I_{total} = \frac{1}{0.08-0}I_{S13}+\frac{1}{0.1-0.08}I_{Tr1} +\frac{1}{0.13-0.1}I_{S38} + \frac{1}{0.16-0.13}I_{Tr2} + \frac{1}{0.2-0.16}I_{Tr3} +\frac{1}{0.23-0.2}I_{Tr4} + \frac{1}{0.25-0.23}I_{Tr5} X_{rms} = \sqrt{I_{total}}

Los valores resultantes son:

Is13:  19.26666666666667
ITr1:  1.1300000000000001
Is38:  143.7
ITr2:  2.175
ITr3:  13.780000000000001
ITr4:  10.59
ITr5:  1.62
Itotal:  5938.333333333333
Xrms:  77.06058222809722

Tarea: literal b

Para Simpson 1/3

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

Para Simpson 3/8

error_{truncamiento} = -\frac{3}{80} h^5 f^{(4)} (z)

Para trapecios

error_{truncar} = -\frac{h^3}{12}f''(z)


Algoritmo en Python

# 3Eva_IT2019_T1 Esfuerzo Cardiaco
import numpy as np
import matplotlib.pyplot as plt

# INGRESO
t  = np.array([0.0,0.04,0.08,0.1,0.11,0.12,0.13,0.16,0.20,0.23,0.25])
ft = np.array([10., 18, 7, -8, 110, -25, 9, 8, 25, 9, 9])

# PROCEDIMIENTO
# Revisar tamaño de paso h
n = len(t)
dt = np.zeros(n, dtype=float)
for i in range(0,n-1,1):
    dt[i]=t[i+1]-t[i]

# Integrales
Is13 = (0.04/3)*((10)**2 + 4*((18)**2) + (7)**2)
ITr1 = (0.02/2)*((7)**2 + (-8)**2)
Is38 = (3/8)*0.01*((-8)**2 + 3*((110)**2) + 3*((-25)**2) + (9)**2)

ITr2 = (0.03/2)*((9)**2 + (8)**2)
ITr3 = (0.04/2)*((8)**2 + (25)**2)
ITr4 = (0.03/2)*((25)**2 + (9)**2)
ITr5 = (0.02/2)*((9)**2 + (9)**2)

Itotal = (1/(0.08-0.0))*Is13 + (1/(0.1-0.08))*ITr1
Itotal = Itotal + (1/(0.13-0.1))*Is38 + (1/(0.16-0.13))*ITr2
Itotal = Itotal + (1/(0.20-0.16))*ITr3 + (1/(0.23-0.20))*ITr4
Itotal = Itotal + (1/(0.25-0.23))*ITr5
Xrms = np.sqrt(Itotal)

# SALIDA
print('tamaño de paso h:')
print(dt)
print('tiempos:')
print(t)
print('ft: ')
print(ft)
print('Is13: ', Is13)
print('ITr1: ', ITr1)
print('Is38: ', Is38)
print('ITr2: ', ITr2)
print('ITr3: ', ITr3)
print('ITr4: ', ITr4)
print('ITr5: ', ITr5)
print('Itotal: ', Itotal)
print('Xrms: ', Xrms)

# Grafica
plt.plot(t,ft)
for i in range(1,n,1):
    plt.axvline(t[i], color='green', linestyle='dashed')
plt.xlabel('tiempo s')
plt.ylabel('valor sensor')
plt.title('Un pulso cardiaco con sensor')
plt.show()

s1Eva_IT2019_T3 Vector perpendicular a plano

Ejercicio: 1Eva_IT2019_T3 Vector perpendicular a plano

Literal a

Para que un vector sea perpendicular a otro, se debe cumplir que
V1.V2 =0.

\begin{bmatrix} 2 \\ -3 \\ a \end{bmatrix} . \begin{bmatrix} b \\ 1 \\ -4 \end{bmatrix} = 0

se obtiene entonces la ecuación:

(2)(b)+(-3)(1)+(a)(-4)=0
2b -3 -4a =0

procediendo de la misma forma con los siguientes pares de vectores:

\begin{bmatrix} 2 \\ -3 \\ a \end{bmatrix} . \begin{bmatrix} 3 \\ c \\ 2 \end{bmatrix} = 0

se obtiene entonces la ecuación:

(2)(3)+(-3)(c)+(a)(2)=0
6 -3c +2a = 0
\begin{bmatrix} b \\ 1 \\ -4 \end{bmatrix} . \begin{bmatrix} 3 \\ c \\ 2 \end{bmatrix} = 2

se obtiene entonces la ecuación:

(b)(3)+(1)(c)+(-4)(2)=2
3b +c -8 =2

se obtiene el sistema de ecuaciones:

\begin{cases}-4a+2b=3 \\ 2a-3c=-6 \\3b+c=10 \end{cases}

Literal b

Se convierte a la forma matricial Ax=B

\begin{bmatrix} -4 && 2 && 0 \\ 2 && 0 && -3 \\ 0 && 3 && 1\end{bmatrix}.\begin{bmatrix} a \\b\\c \end{bmatrix} = \begin{bmatrix} 3 \\ -6 \\ 10 \end{bmatrix}

se crea la matriz aumentada:

\begin{bmatrix} -4 && 2 && 0 && 3 \\ 2 && 0 && -3 &&-6 \\ 0 && 3 && 1 && 10 \end{bmatrix}

se pivotea por filas buscando hacerla diagonal dominante:

\begin{bmatrix} -4 && 2 && 0 && 3 \\ 0 && 3 && 1 && 10 \\ 2 && 0 && -3 &&-6 \end{bmatrix}

se aplica el algoritmo de eliminación hacia adelante:
1ra Iteración

la eliminación del primer término columna es necesario solo para la tercera fila:

\begin{bmatrix} -4 && 2 && 0 && 3 \\ 0 && 3 && 1 && 10 \\ 2 -(-4)\Big( \frac{2}{-4}\Big) && 0-2\Big(\frac{2}{-4}\Big) && -3 -0\Big(\frac{2}{-4}\Big) &&-6 -3\Big(\frac{2}{-4}\Big) \end{bmatrix} \begin{bmatrix} -4 && 2 && 0 && 3 \\ 0 && 3 && 1 && 10 \\ 0 && 1 && -3 && -4.5 \end{bmatrix}

2da Iteración

\begin{bmatrix} -4 && 2 && 0 && 3 \\ 0 && 3 && 1 && 10 \\ 0 && 1 -3\Big(\frac{1}{3}\Big) && -3-(1)\Big(\frac{1}{3}\Big) && -4.5-10\Big(\frac{1}{3}\Big) \end{bmatrix} \begin{bmatrix} -4 && 2 && 0 && 3 \\ 0 && 3 && 1 && 10 \\ 0 && 0 && -\frac{10}{3} && -7.8333 \end{bmatrix}

Aplicando eliminación hacia atras

(-10/3)c = -7.8333
c = -7.8333(-3/10) = 2.35

3b +c = 10
b= (10-c)/3 = (10-2.35)/3 = 2.55

-4a +2b =3
a= (3-2b)/(-4) = (3-2(2.55))/(-4) = 0.525

como resultado, los vectores buscados:

v1 = (2,-3,0.525)
v2 = (2.55,1,-4)
v3 = (3,2.35,2)

comprobando resultados:

>>> np.dot((2,-3,0.525),(2.55,1,-4))
-4.440892098500626e-16
>>> np.dot((2,-3,0.525),(3,2.35,2))
-6.661338147750939e-16
>>> np.dot((2.55,1,-4),(3,2.35,2))
2.0

Los dos primeros resultados son muy cercanos a cero, por lo que se los considera válidos

literal c

Para usar el método de Jacobi, se despeja una de las variables de cada ecuación:

\begin{cases} a = \frac{2b -3}{4} \\b = \frac{10-c}{3} \\c = \frac{2a+6}{3} \end{cases}

usando el vector x(0) = [0,0,0]

1ra iteración

a = \frac{2b -3}{4} = \frac{2(0) -3}{4} = -\frac{3}{4} b = \frac{10-c}{3} = \frac{10-0}{3} = \frac{10}{3} c = \frac{2a+6}{3 }= \frac{2(0)+6}{3} = 2

x(1) = [-3/4,10/3,2]
diferencias = [-3/4-0,10/3-0,2-0] = [-3/4,10/3,2]
error = max|diferencias| = 10/3 = 3.3333

2da iteración

a = \frac{2b -3}{4} = \frac{2(10/3) -3}{4} = \frac{11}{12} b = \frac{10-c}{3} = \frac{10-2}{3} = \frac{8}{3} c = \frac{2a+6}{3} = \frac{2(-3/4)+6}{3} = \frac{3}{2}

x(2) = [11/12,8/3,3/2]
diferencias = [11/12-(-3/4),8/3-10/3,3/2-2] = [20/12, -2/3, -1/2]
error = max|diferencias| = 5/3= 1.666

3ra iteración

a = \frac{2b -3}{4} = \frac{2(8/3)-3}{4} = \frac{7}{12} b = \frac{10-c}{3} = \frac{10-3/2}{3} = \frac{17}{6} c = \frac{2a+6}{3} = \frac{2(11/12)+6}{3} = 2.6111

x(3) = [7/12, 17/6, 2.6111]
diferencias = [7/12-11/12, 17/6-8/3, 2.6111-3/2] = [-1/3, 1/6, 1.1111]
error = max|diferencias| = 1.1111

Los errores disminuyen en cada iteración, por lo que el método converge,
si se analiza en número de condición de la matriz A es 2, lo que es muy cercano a 1, por lo tanto el método converge.


Revisión de resultados

Resultados usando algoritmos desarrollados en clase:

matriz aumentada: 
[[-4.   2.   0.   3. ]
 [ 0.   3.   1.  10. ]
 [ 2.   0.  -3.  -6. ]]
Elimina hacia adelante
[[-4.   2.   0.   3. ]
 [ 0.   3.   1.  10. ]
 [ 0.   1.  -3.  -4.5]]
Elimina hacia adelante
[[-4.   2.   0.        3.      ]
 [ 0.   3.   1.       10.      ]
 [ 0.   0.  -3.333333 -7.833333]]
Elimina hacia adelante
[[-4.   2.   0.        3.      ]
 [ 0.   3.   1.       10.      ]
 [ 0.   0.  -3.333333 -7.833333]]
Elimina hacia atras
[[ 1.  -0.  -0.   0.525]
 [ 0.   1.   0.   2.55 ]
 [-0.  -0.   1.   2.35 ]]
el vector solución X es:
[[0.525]
 [2.55 ]
 [2.35 ]]
verificar que A.X = B
[[ 3.]
 [10.]
 [-6.]]

Número de condición de A: 2.005894

los resultados para [a,b,c]:
[0.525 2.55  2.35 ]

producto punto entre vectores:
v12:  0.0
v13:  1.3322676295501878e-15
v23:  2.0

Algoritmos en Python:

# 1Eva_IT2019_T3 Vector perpendicular a plano
import numpy as np

# INGRESO
A = np.array([[-4.,2,0],
              [ 0., 3,1],
              [ 2.,0,-3]])
B = np.array([3.,10,-6])

# PROCEDIMIENTO
B = np.transpose([B])

casicero = 1e-15 # 0
AB = np.concatenate((A,B),axis=1)
tamano = np.shape(AB)
n = tamano[0]
m = tamano[1]

print('matriz aumentada: ')
print(AB)
# Gauss elimina hacia adelante
# tarea: verificar términos cero
for i in range(0,n,1):
    pivote = AB[i,i]
    adelante = i+1 
    for k in range(adelante,n,1):
        if (np.abs(pivote)>=casicero):
            factor = AB[k,i]/pivote
            AB[k,:] = AB[k,:] - factor*AB[i,:]
    print('Elimina hacia adelante')
    print(AB)

# Gauss-Jordan elimina hacia atras
ultfila = n-1
ultcolumna = m-1
for i in range(ultfila,0-1,-1):
    # Normaliza a 1 elemento diagonal
    AB[i,:] = AB[i,:]/AB[i,i]
    pivote = AB[i,i] # uno
    # arriba de la fila i
    atras = i-1 
    for k in range(atras,0-1,-1):
        if (np.abs(AB[k,i])>=casicero):
            factor = pivote/AB[k,i]
            AB[k,:] = AB[k,:]*factor - AB[i,:]
        else:
            factor= 'division para cero'
print('Elimina hacia atras')
print(AB)

X = AB[:,ultcolumna]

# Verifica resultado
verifica = np.dot(A,X)

# SALIDA
print('el vector solución X es:')
print(np.transpose([X]))

print('verificar que A.X = B')
print(np.transpose([verifica]))

numcond = np.linalg.cond(A)

# para comprobar respuestas

v1 = [2,-3,X[0]]
v2 = [X[1],1,-4]
v3 = [3,X[2],2]

v12 = np.dot(v1,v2)
v13 = np.dot(v1,v3)
v23 = np.dot(v2,v3)
      
# SALIDA
print('\n Número de condición de A: ', numcond)

print('\n los resultados para [a,b,c]:')
print(X)

print('\n productos punto entre vectores:')
print('v12: ',v12)
print('v13: ',v13)
print('v23: ',v23)

Tarea, usar el algoritmo de Jacobi usado en el taller correspondiente.

s1Eva_IT2019_T2 Catenaria cable

Ejercicio: 1Eva_IT2019_T2 Catenaria cable

Las fórmulas con las que se requiere trabajar son:

y = \frac{T_A}{w} cosh \Big( \frac{w}{T_A}x \Big) + y_0 - \frac{T_A}{w}

Donde la altura y del cable está en función de la distancia x.

Además se tiene que:

cosh(z) = \frac{e^z+ e^{-z}}{2}

que sustituyendo la segunda en la primera se convierte en:

y = \frac{T_A}{w} \frac{e^{\frac{w}{T_A}x}+ e^{-\frac{w}{T_A}x}}{2} + y_0 - \frac{T_A}{w}

y usando los valores del enunciado w=12, y0=6 , y=15, x=50 se convierte en:

15 = \frac{T_A}{12} \frac{e^{\frac{12}{T_A}50}+ e^{-\frac{12}{T_A}50}}{2} + 6 - \frac{T_A}{12}

simplificando, para usar el método de búsqueda de raices:

\frac{1}{2}\frac{T_A}{12} e^{\frac{12}{T_A}50} + \frac{1}{2}\frac{T_A}{12} e^{-\frac{12}{T_A}50} - \frac{T_A}{12} - 9 = 0

cambiando la variable \frac{12}{T_A}=x

\frac{1}{2x} e^{50x} + \frac{1}{2x} e^{-50x} - \frac{1}{x}-9=0

la función a usar para la búsqueda de raices es:

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

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

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

por lo que se determina:

f'(x)= - \frac{1}{2x^2}e^{50x} + \frac{1}{2x}(50) e^{50x} + - \frac{1}{2x^2} e^{-50x} + \frac{1}{2x}(-50)e^{-50x} + \frac{1}{x^2} f'(x)= -\frac{1}{2x^2}[e^{50x}+e^{-50x}] + + \frac{25}{x}[e^{50x}-e^{-50x}] +\frac{1}{x^2} f'(x)= \Big[\frac{25}{x} -\frac{1}{2x^2}\Big]\Big[e^{50x}+e^{-50x}\Big] +\frac{1}{x^2}

Con lo que se puede inicar las iteraciones.

Por no disponer de valor inicial para TA, considere que el cable colgado no debería tener tensión TA=0 N, pues en la forma x=12/TA se crea una indeterminación. Si no dispone de algún criterio para seleccionar el valor de TA puede iniciar un valor positivo, por ejemplo 120 con lo que el valor de x0=12/120=0.1

Iteración 1

f(0.1)=\frac{1}{2(0.1)} e^{50(0.1)} + \frac{1}{2(0.1)} e^{-50(0.1)} - \frac{1}{0.1}-9 =723.0994 f'(0.1)=\Big[\frac{25}{0.1} - \frac{1}{2(0.1)^2}\Big]\Big[e^{50(0.1)}+e^{-50(0.1)}\Big] + +\frac{1}{(0.1)^2} = 29780.61043 x_{1} = 0.1 -\frac{723.0994}{29780.61043} = 0.07571

error = | x1 – x0| = | 0.07571 – 0.1| = 0.02428

Iteración 2

f(0.07571)=\frac{1}{2(0.07571)} e^{50(0.07571)}+ + \frac{1}{2(0.07571)} e^{-50(0.07571)} - \frac{1}{0.07571}-9 = 269.0042 f'(0.07571)= \Big[\frac{25}{0.07571} -\frac{1}{2(0.07571)^2}\Big]. .\Big[e^{50(0.07571)}+e^{-50(0.07571)}\Big] + +\frac{1}{(0.07571)^2} = 10874.0462 x_{2} = 0.07571 -\frac{269.0042}{10874.0462} = 0.05098

error = | x2 – x1| = |0.05098- 0.02428| = 0.02473

Iteración 3

f(0.05098) = 97.6345 f'(0.05098) = 4144.1544 x_{3} = 0.0274

error = | x3 – x2| = |0.05098- 0.0274| = 0.0236

finalmente después de varias iteraciones, la raiz se encuentra en: 0.007124346154337298

que convitiendo

T_A = \frac{12}{x} = \frac{12}{0.0071243461} = 1684.36 N

Revisión de resultados

Usando como base los algoritmos desarrollados en clase:

['xi', 'xnuevo', 'tramo']
[0.1    0.0757 0.0243]
[0.0757 0.051  0.0247]
[0.051  0.0274 0.0236]
[0.0274 0.0111 0.0163]
[0.0111 0.0072 0.0039]
[7.2176e-03 7.1244e-03 9.3199e-05]
[7.1244e-03 7.1243e-03 3.8351e-08]
raiz en:  0.007124346154337298
TA = 12/x =  1684.365096815854

Algoritmos Python usando el procedimiento de:

http://blog.espol.edu.ec/analisisnumerico/2-3-1-newton-raphson-ejemplo01/

# 1Eva_IT2019_T2 Catenaria cable
import numpy as np
import matplotlib.pyplot as plt

# INGRESO
a = 0.001
b = 0.1
muestras = 51

x0 = 0.1
tolera = 0.00001

fx = lambda x: 0.5*(1/x)*np.exp(50*x) + 0.5*(1/x)*np.exp(-50*x)-1/x -9
dfx = lambda x: -0.5*(1/(x**2))*(np.exp(50*x)+np.exp(-50*x)) + (25/x)*(np.exp(50*x)-np.exp(-50*x)) + 1/(x**2)

# PROCEDIMIENTO
tabla = []
tramo = abs(2*tolera)
xi = x0
while (tramo>=tolera):
    xnuevo = xi - fx(xi)/dfx(xi)
    tramo = abs(xnuevo-xi)
    tabla.append([xi,xnuevo,tramo])
    xi = xnuevo

tabla = np.array(tabla)
n=len(tabla)

TA = 12/xnuevo

# para la gráfica
xp = np.linspace(a,b,muestras)
fp = fx(xp)

# SALIDA
print(['xi', 'xnuevo', 'tramo'])
np.set_printoptions(precision = 4)
for i in range(0,n,1):
    print(tabla[i])
print('raiz en: ', xi)
print('TA = 12/x = ', TA)

# Grafica
plt.plot(xp,fp)
plt.xlabel('x=12/TA')
plt.ylabel('f(x)')
plt.axhline(0, color = 'green')
plt.grid()
plt.show()

s1Eva_IT2019_T1 Oxígeno y temperatura en agua

Ejercicio: 1Eva_IT2019_T1 Oxígeno y temperatura en agua

Literal a

Se requiere un polinomio de grado 3 siendo el eje x correspondiente a temperatura. Son necesarios 4 puntos de referencia alrededor de 15 grados, dos a la izquierda y dos a la derecha.

Se observa que los datos en el eje x son equidistantes, h=8, y ordenados en forma ascendente, se cumple con los requisitos para usar diferencias finitas avanzadas. que tiene la forma de:

p_n (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) + + \frac{\Delta^3 f_0}{3!h^3} (x - x_0)(x - x_1)(x - x_2) + \text{...}

Tabla

xi f[xi] f[x1,x0] f[x2,x1,x0] f[x2,x1,x0] f[x3,x2,x1,x0]
8 11.5 9.9-11.5=
-1.6
-1.5-(-1.6) =
0.1
0.4-0.1=
0.3
16 9.9 8.4-9.9=
-1.5
-1.1-(1.5)=
0.4
24 8.4 7.3-8.4=
-1.1
32 7.3

Con lo que el polinomio buscado es:

p_3 (x) = 11.5 + \frac{-1.6}{8} (x - 8) + + \frac{0.1}{2!8^2} (x - 8)(x - 16) + \frac{0.3}{3!8^3} (x - 8)(x - 16)(x - 24)

Resolviendo y simplificando el polinomio, se puede observar que al aumentar el grado, la constante del término disminuye.

p_3(x)=12.9- 0.15 x - 0.00390625 x^2 + 0.00009765625 x^3

para el cálculo del error se puede usar un término adicional del polinomio, añadiendo un punto más a la tabla de diferencia finitas. Se evalúa éste término y se estima el error que dado que el término de grado 3 es del orden de 10-5, el error será menor. (Tarea)

p_3(15)=12.9- 0.15 (15) - 0.00390625 (15)^2 + 0.00009765625 (15)^3

Evaluando el polinomio en temperatura = 15:

p3(15) = 10.1006835937500

literal b

se deriva el polinomio del literal anterior y se evalua en 16:

p'_3(x)=0- 0.15 - 0.00390625 (2) x + 0.00009765625 (3)x^2 p'_3(16)=0- 0.15 - 0.00390625 (2)(16) + 0.00009765625 (3)(16)^2

p’3(16) = -0.20

literal c

El valor de oxígeno usado como referencia es 9, cuyos valores de temperatura se encuentran entre 16 y 24 que se toman como rango inicial de búsqueda [a,b]. Por lo que el polinomio se iguala a 9 y se crea la forma estandarizada del problema:

p_3(x)=9 9 = 12.9- 0.15 x - 0.00390625 x^2 + 0.00009765625 x^3 12.9- 0.15 x - 0.00390625 x^2 + 0.00009765625 x^3 -9 = 0 f(x) = 3.9- 0.15 x - 0.00390625 x^2 + 0.00009765625 x^3

Para mostrar el procedimiento se realizan solo tres iteraciones,

1ra Iteración
a=16 , b = 24, c = (16+24)/2 = 20
f(a) = 0.9, f(b) = -0.6, f(c) = 0.011
error = |24-16| = 8
como f(c) es positivo, se mueve el extremo f(x) del mismo signo, es decir a

2da Iteración
a=20 , b = 24, c = (20+24)/2 = 22
f(a) = 0.119, f(b) = -0.6, f(c) = -0.251
error = |24-20|= 4
como f(c) es negativo, se mueve el extremo f(x) del mismo signo, b

3ra Iteración
a=20 , b = 22, c = (20+22)/2 = 21
f(a) = 0.119, f(b) = -0.251, f(c) = -0.068
error = |22-20| = 2
como f(c) es negativo, se mueve el extremo f(x) del mismo signo, b
y así sucesivamente hasta que error< que 10-3

Usando el algoritmo en python se obtendrá la raiz en 20.632 con la tolerancia requerida.


Revisión de resultados

Usando como base los algoritmos desarrollados en clase:

Literal a
[[ 1.   8.  11.5 -1.6  0.1  0.3  0. ]
 [ 2.  16.   9.9 -1.5  0.4  0.   0. ]
 [ 3.  24.   8.4 -1.1  0.   0.   0. ]
 [ 4.  32.   7.3  0.   0.   0.   0. ]]
9.765625e-5*x**3 - 0.00390625*x**2 - 0.15*x + 12.9
p(15) =  10.1006835937500
Literal b
0.00029296875*x**2 - 0.0078125*x - 0.15
dp(16) = -0.200000000000000

literal c
[ i, a, c, b, f(a), f(c), f(b), paso]
1 16.000 20.000 24.000 0.900 0.119 -0.600 8.000 
2 20.000 22.000 24.000 0.119 -0.251 -0.600 4.000 
3 20.000 21.000 22.000 0.119 -0.068 -0.251 2.000 
4 20.000 20.500 21.000 0.119 0.025 -0.068 1.000 
5 20.500 20.750 21.000 0.025 -0.022 -0.068 0.500 
6 20.500 20.625 20.750 0.025 0.001 -0.022 0.250 
7 20.625 20.688 20.750 0.001 -0.010 -0.022 0.125 
8 20.625 20.656 20.688 0.001 -0.004 -0.010 0.062 
9 20.625 20.641 20.656 0.001 -0.002 -0.004 0.031 
10 20.625 20.633 20.641 0.001 -0.000 -0.002 0.016 
11 20.625 20.629 20.633 0.001 0.001 -0.000 0.008 
12 20.629 20.631 20.633 0.001 0.000 -0.000 0.004 
13 20.631 20.632 20.633 0.000 0.000 -0.000 0.002 
14 20.632 20.632 20.633 0.000 0.000 -0.000 0.001 
raiz:  20.63232421875

Algoritmos Python usando la funcion de interpolación y un procedimiento encontrado en:

http://blog.espol.edu.ec/analisisnumerico/5-1-1-diferencias-finitas-avanzadas-polinomio/

http://blog.espol.edu.ec/analisisnumerico/2-1-1-biseccion-ejemplo01/

# 1Eva_IT2019_T1 Oxígeno y temperatura en mar
import numpy as np
import matplotlib.pyplot as plt
import sympy as sym

def interpola_dfinitas(xi,fi):
    '''
    Interpolación de diferencias finitas avanzadas
    resultado: polinomio en forma simbólica
    '''
    # Tabla de diferencias finitas
    titulo = ['i','xi','fi']
    n = len(xi)
    # cambia a forma de columnas
    i = np.arange(1,n+1,1)
    i = np.transpose([i])
    xi = np.transpose([xi])
    fi = np.transpose([fi])
    # Añade matriz de diferencias
    dfinita = np.zeros(shape=(n,n),dtype=float)
    tabla = np.concatenate((i,xi,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 x equidistantes
    dfinita = tabla[:,3:]
    n = len(dfinita)
    x = sym.Symbol('x')
    h = xi[1,0]-xi[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*(x-xi[f])
        polinomio = polinomio + termino*factor
    # simplifica polinomio, multiplica los (x-xi)
    polinomio = polinomio.expand()
    
    return(tabla, polinomio)

# INGRESO
tm = [0.,8,16,24,32,40]
ox = [14.6,11.5,9.9,8.4,7.3,6.4]

xi = [8,16,24,32]
fi = [11.5,9.9,8.4,7.3]

# PROCEDIMIENTO
x = sym.Symbol('x')
# literal a
tabla, polinomio = interpola_dfinitas(xi,fi)
p15 = polinomio.subs(x,15)
# literal b
deriva = polinomio.diff(x,1)
dp16 = deriva.subs(x,16)

px =  sym.lambdify(x,polinomio)
xk = np.linspace(np.min(xi),np.max(xi))
pk = px(xk)

# SALIDA
print('Literal a')
print(tabla)
print(polinomio)
print('p(15) = ',p15)
print('Literal b')
print(deriva)
print('dp(16) =', dp16)

# gráfica
plt.plot(tm,ox,'ro')
plt.plot(xk,pk)
plt.axhline(9,color="green")
plt.xlabel('temperatura')
plt.ylabel('concentracion de oxigeno')
plt.grid()
plt.show()

# --------literal c ------------

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

# se convierte forma de símbolos a numéricos
buscar = polinomio-9
fx = sym.lambdify(x,buscar)

# INGRESO
a = 16
b = 24
tolera = 0.001

# PROCEDIMIENTO
tabla = []
tramo = b-a

fa = fx(a)
fb = fx(b)
i = 1
while (tramo > tolera):
    c = (a+b)/2
    fc = fx(c)
    tabla.append([i,a,c,b,fa,fc,fb,tramo])
    i = i+1
                 
    cambia = np.sign(fa)*np.sign(fc)
    if (cambia<0):
        b = c
        fb = fc
    else:
        a=c
        fa = fc
    tramo = b-a
c = (a+b)/2
fc = fx(c)
tabla.append([i,a,c,b,fa,fc,fb,tramo])
tabla = np.array(tabla)

raiz = c

# SALIDA
print('\n literal c')
np.set_printoptions(precision = 4)
print('[ i, a, c, b, f(a), f(c), f(b), paso]')
# print(tabla)

# Tabla con formato
n=len(tabla)
for i in range(0,n,1):
    unafila = tabla[i]
    formato = '{:.0f}'+' '+(len(unafila)-1)*'{:.3f} '
    unafila = formato.format(*unafila)
    print(unafila)
    
print('raiz: ',raiz)

s3Eva_IIT2018_T2 Drenar tanque cilíndrico

Ejercicio: 3Eva_IIT2018_T2 Drenar tanque cilíndrico

La ecuación a desarrollar es:

\frac{\delta y}{\delta t} = -k\sqrt{y}

con valores de k =0.5, y(0)=9


Formula de Taylor con término de error:

P_{n}(x) = \sum_{k=0}^{n} \frac{f^{(k)}(x_0)}{k!} (x-x_0)^k P_{n}(x) = f(x_0)+\frac{f'(x_0)}{1!} (x-x_0) + + \frac{f''(x_0)}{2!}(x-x_0)^2 + + \frac{f'''(x_0)}{3!}(x-x_0)^3 + \text{...}

Se requiere la 2da y 3ra derivadas:

\frac{\delta^2 y}{\delta t^2} = -k\frac{1}{2} y^{(\frac{1}{2}-1)} = -\frac{k}{2} y^{-\frac{1}{2}} \frac{\delta^3 y}{\delta t^3} = -\frac{k}{2}\Big(-\frac{1}{2}\Big) y^{(-\frac{1}{2}-1)} = \frac{k}{4} y^{-\frac{3}{2}}

con lo que inicia las iteraciones y cálculo del error, con avance de 0.5 para t.


t=0 , y(0) = 9


t = 0.5

\frac{\delta y(0)}{\delta t} = -(0.5)\sqrt{9} = -1.5 \frac{\delta^2 y(0)}{\delta t^2} = -\frac{0.5}{2} 9^{-\frac{1}{2}} = - 0.08333 \frac{\delta^3 y(0)}{\delta t^3} = \frac{0.5}{4} 9^{-\frac{3}{2}} = 0.004628 P_{2}(0.5) = 9 - 1.5 (0.5-0) + \frac{-0.08333}{2}(0.5-0)^2 P_{2}(0.5) = 8.2395

Error orden de:

Error = \frac{0.004628}{3!}(0.5-0)^3 = 9.641 . 10^{-5}

t = 1

\frac{\delta y(0.5)}{\delta t} = -(0.5)\sqrt{8.2395} = -1.4352 \frac{\delta^2 y(0.5)}{\delta t^2} = -\frac{0.5}{2} (8.2395)^{-\frac{1}{2}} = - 0.08709 \frac{\delta^3 y(0.5)}{\delta t^3} = \frac{0.5}{4} (8.2395)^{-\frac{3}{2}} = 0.005285 P_{2}(1) = 8.2395 - 1.4352(1-0.5) + \frac{-0.08709}{2}(1-0.5)^2 P_{2}(1) = 7.5110

Error orden de:

Error = \frac{0.005285}{3!}(1-0.5)^3 = 4.404 . 10^{-4}

t = 1.5

\frac{\delta y(1)}{\delta t} = -(0.5)\sqrt{7.5110} = -1.3703 \frac{\delta^2 y(1)}{\delta t^2} = -\frac{0.5}{2} (7.5110)^{-\frac{1}{2}} = - 0.09122 \frac{\delta^3 y(1)}{\delta t^3} = \frac{0.5}{4} (7.5110)^{-\frac{3}{2}} = 0.006072 P_{2}(1.5) = 7.5110 - 1.3703(1.5-1) + \frac{-0.09122}{2}(1.5-1)^2 P_{2}(1.5) = 6.8144

Error orden de:

Error = \frac{0.006072}{3!}(1.5-1)^3 = 1.4 . 10^{-4}

t = 2

\frac{\delta y(1.5)}{\delta t} = -(0.5)\sqrt{6.8144} = -1.3052 \frac{\delta^2 y(1.5)}{\delta t^2} = -\frac{0.5}{2} (6.8144)^{-\frac{1}{2}} = - 0.09576 \frac{\delta^3 y(1.5)}{\delta t^3} = \frac{0.5}{4} (6.8144)^{-\frac{3}{2}} = 0.007026 P_{2}(2) = 6.8144 - 1.3052 (2-1.5) - \frac{0.09576}{2}(2-1.5)^2 P_{2}(2) = 6.1498

Error orden de:

Error = \frac{0.007026}{3!}(2-1.5)^3 = 1.4637 . 10^{-4}

Se estima que el próximo término pasa debajo de 6 pies.
Por lo que estima esperar entre 2 y 2.5 minutos.

resultados usando el algoritmo:

ti, p_i,  error
[[0.00000000e+00 9.00000000e+00 0.00000000e+00]
 [5.00000000e-01 8.23958333e+00 9.64506173e-05]
 [1.00000000e+00 7.51107974e+00 1.10105978e-04]
 [1.50000000e+00 6.81451855e+00 1.26507192e-04]
 [2.00000000e+00 6.14993167e+00 1.46391550e-04]
 [2.50000000e+00 5.51735399e+00 1.70751033e-04]]

Algoritmo en Python

# 3Eva_IIT2018_T2 Drenar tanque cilíndrico
import numpy as np
import matplotlib.pyplot as plt

# INGRESO
y0 = 9
t0 = 0
buscar = 6
k = 0.5
h = 0.5

dy  = lambda t,y: -k*np.sqrt(y)
d2y = lambda t,y: -(k/2)*(y**(-1/2))
d3y = lambda t,y: (k/4)*(y**(-3/2))

# PROCEDIMIENTO
resultado = [[t0,y0,0]]
yi = y0
ti = t0
while not(yi<buscar):
    ti = ti+h
    dyi = dy(ti,yi)
    d2yi = d2y(ti,yi)
    d3yi = d3y(ti,yi)
    p_i = yi +dyi*(h) + (d2yi/2)*(h**2)
    errado = (d3yi/6)*(h**3)
    yi = p_i
    resultado.append([ti,p_i,errado])
resultado = np.array(resultado)

# SALIDA
print('ti, p_i,  error')
print(resultado)

# Grafica
plt.plot(resultado[:,0],resultado[:,1])
plt.ylabel('nivel de agua')
plt.xlabel('tiempo')
plt.grid()
plt.show()

s2Eva_IIT2018_T3 EDP

Ejercicio: 2Eva_IIT2018_T3 EDP

Se indica en el enunciado que b = 0

\frac{\delta u}{\delta t} = \frac{\delta ^2 u}{\delta x^2} + b\frac{\delta u}{\delta x}

simplificando la ecuación a:

\frac{\delta u}{\delta t} = \frac{\delta ^2 u}{\delta x^2}

Reordenando la ecuación a la forma estandarizada:

\frac{\delta ^2 u}{\delta x^2} = \frac{\delta u}{\delta t}

Seleccione un método: explícito o implícito.
Si el método es explícito, las diferencias finitas a usar son hacia adelante y centrada:

U'(x_i,t_j) = \frac{U(x_i,t_{j+1})-U(x_i,t_j)}{\Delta t} + O(\Delta t) U''(x_i,t_j) = \frac{U(x_{i+1},t_j)-2U(x_{i},t_j)+U(x_{i-1},t_j)}{\Delta x^2} + O(\Delta x^2)

como referencia se usa la gráfica.

Se selecciona la esquina inferior derecha como 0,  por la segunda ecuación de condiciones y facilidad de cálculo. (No hubo indicación durante el examen que muestre lo contrario)

condiciones de frontera U(0,t)=0, U(1,t)=1
condiciones de inicio U(x,0)=0, 0≤x≤1

aunque lo más recomendable sería cambiar la condición de inicio a:

condiciones de inicio U(x,0)=0, 0<x<1

Siguiendo con el tema de la ecuación, al reemplazar las diferencias finitas en la ecuación:


\frac{U(x_{i+1},t_j)-2U(x_{i},t_j)+U(x_{i-1},t_j)}{\Delta x^2} = = \frac{U(x_i,t_{j+1})-U(x_i,t_j)}{\Delta t}

se reagrupan los términos que son constantes y los términos de error se acumulan:

\frac{\Delta t}{\Delta x^2} \Big[U(x_{i+1},t_j)-2U(x_i,t_j)+U(x_{i-1},t_j) \Big] = U(x_i,t_{j+1})-U(x_i,t_j)

siendo,

\lambda= \frac{\Delta t}{\Delta x^2} error \cong O(\Delta t) + O(\Delta x^2)

continuando con la ecuación, se simplifica la escritura usando sólo los índices i,j y se reordena de izquierda a derecha como en la gráfica

\lambda \Big[U[i-1,j]-2U[i,j]+U[i+1,j] \Big] = U[i,j+1]-U]i,j] \lambda U[i-1,j]+(-2\lambda+1)U[i,j]+\lambda U[i+1,j] = U[i,j+1] U[i,j+1] = \lambda U[i-1,j]+(-2\lambda+1)U[i,j]+\lambda U[i+1,j] U[i,j+1] = P U[i-1,j]+QU[i,j]+R U[i+1,j] P=R = \lambda Q = -2\lambda+1

En las iteraciones, el valor de P,Q y R se calculan a partir de λ ≤ 1/2

iteraciones: j=0, i=1

U[1,1] = P*0+Q*0+R*0 = 0

j=0, i=2

U[2,1] = P*0+Q*0+R*0=0

j=0, i=3

U[3,1] = P*0+Q*0+R*1=R=\lambda=\frac{1}{2}

iteraciones: j=1, i=1

U[1,2] = P*0+Q*0+R*0 = 0

j=1, i=2

U[2,2] = P*0+Q*0+R*\lambda = \lambda ^2 = \frac{1}{4}

j=1, i=3

U[3,2] = P*0+Q*\frac{1}{4}+R (\lambda) U[3,2] = (-2\lambda +1) \frac{1}{4}+\lambda^2 = \Big(-2\frac{1}{2}+1\Big) \frac{1}{4}+\Big(\frac{1}{2}\Big)^2 U[3,2] =0\frac{1}{4} + \frac{1}{4} = \frac{1}{4}

Literal b. Para el cálculo del error:

\lambda \leq \frac{1}{2} \frac{\Delta t}{\Delta x^2} \leq \frac{1}{2} \Delta t \leq \frac{\Delta x^2}{2}

en el enunciado se indica h = 0.25 = ¼ = Δ x

\Delta t \leq \frac{(1/4)^2}{2} = \frac{1}{32} error \cong O(\Delta t) + O(\Delta x^2) error \cong \frac{\Delta x^2}{2}+ \Delta x^2 error \cong \frac{3}{2}\Delta x^2 error \cong \frac{3}{2}( \frac{1}{4})^2 error \cong \frac{3}{32} = 0.09375

s2Eva_IIT2018_T2 Kunge Kutta 2do Orden x»

Ejercicio: 2Eva_IIT2018_T2 Kunge Kutta 2do Orden x»

\frac{\delta ^2 x}{\delta t^2} + 5t\frac{\delta x}{\delta t} +(t+7)\sin (\pi t) = 0 x'' + 5tx' +(t+7)\sin (\pi t) = 0 x'' = -5tx' +(t+7)\sin (\pi t) = 0

si se usa z=x’

z' = -5tz +(t+7)\sin (\pi t) = 0

se convierte en:
f(t,x,z) = x’ = z
g(t,x,z) = x» = z’ = -5tz +(t+7)sin (π t) = 0

Donde se aplica el algoritmo de Runge Kutta
http://blog.espol.edu.ec/analisisnumerico/8-2-2-runge-kutta-d2y-dx2/

   t,              x,              z
[[ 0.00000000e+00  6.00000000e+00  1.50000000e+00]
 [ 2.00000000e-01  6.30000000e+00  1.77320538e+00]
 [ 4.00000000e-01  6.70381805e+00  2.26987703e+00]
 [ 6.00000000e-01  7.20775473e+00  2.41163944e+00]
 [ 8.00000000e-01  7.68994485e+00  1.90531839e+00]
 [ 1.00000000e+00  8.01027755e+00  9.52659193e-01]
 [ 1.20000000e+00  8.10554347e+00 -5.65431040e-03]
 [ 1.40000000e+00  8.00869435e+00 -6.09147239e-01]
 [ 1.60000000e+00  7.81236802e+00 -7.16247408e-01]
 [ 1.80000000e+00  7.62013640e+00 -3.92947221e-01]
 [ 2.00000000e+00  7.50882725e+00  1.63598524e-01]]

Instrucciones en Python

# 2Eva_IIT2018_T2 Kunge Kutta 2do Orden x''
import numpy as np

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

# PROGRAMA
# INGRESO
f = lambda t,x,z: z
g = lambda t,x,z: -5*t*z+(t+7)*np.sin(np.pi*t)
t0 = 0
x0 = 6
z0 = 1.5
h = 0.2
muestras = 10

# PROCEDIMIENTO
tabla = rungekutta2_fg(f,g,t0,x0,z0,h,muestras)

# SALIDA
print(tabla)
# GRAFICA
import matplotlib.pyplot as plt
plt.plot(tabla[:,0],tabla[:,1])
plt.xlabel('t')
plt.ylabel('x(t)')
plt.show()