s3Eva_IIT2019_T1 Lanzamiento de Cohete

Ejercicio: 3Eva_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:

# 3Eva_IIT2019_T1 Lanzamiento de Cohete
# 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()

s3Eva_IT2019_T3 Difusión en sólidos

Ejercicio: 3Eva_IT2019_T3 Difusión en sólidos

Siguiendo el procedimiento planteado en la sección EDP parabólicas, se plantea la malla del ejercicio:

Para plantear la ecuación en forma discreta:

\frac{\phi_{i,j+1}-\phi_{i,j}}{\Delta t}=D\frac{\phi_{i+1,j}-2\phi_{i,j}+\phi_{i-1,j}}{(\Delta x)^2}

y resolver usando el método explícito para ecuaciones parabólicas, obteniendo el siguiente resultado:

\phi_{i,j+1}-\phi_{i,j}=D\frac{\Delta t }{\Delta x^2}(\phi_{i+1,j}-2\phi_{i,j}+\phi_{i-1,j}) \lambda = D\frac{\Delta t }{\Delta x^2} \phi_{i,j+1}-\phi_{i,j}=\lambda (\phi_{i+1,j}-2\phi_{i,j}+\phi_{i-1,j}) \phi_{i,j+1} =\lambda \phi_{i+1,j}-2\lambda\phi_{i,j}+\lambda\phi_{i-1,j}+\phi_{i,j} \phi_{i,j+1} =\lambda \phi_{i+1,j}(1-2\lambda)\phi_{i,j}+\lambda\phi_{i-1,j} \phi_{i,j+1} =P \phi_{i+1,j}+Q\phi_{i,j}+R\phi_{i-1,j}

siendo:
P = λ = 0.16 (Δx/100)/Δx2 = 0.0016/Δx = 0.0016/0.02=0.08
Q = 1-2λ = 1-2*(0.08) = 0.84
R = λ =0.08

\phi_{i,j+1} =0.08 \phi_{i+1,j}+ 0.84\phi_{i,j}+0.08\phi_{i-1,j}

Iteración 1 en tiempo:
i=1, j=0

\phi_{1,1} =0.08 \phi_{2,0}+ 0.84\phi_{1,0}+0.08\phi_{0,0} \phi_{1,1} =0.08 (0)+ 0.84(0)+0.08(5)=0.4

i=2,j=0

\phi_{2,1} =0.08 \phi_{3,0}+ 0.84\phi_{2,0}+0.08\phi_{1,0} = 0

Para los proximos valores i>2, todos los resultados son 0

Iteración 2 en tiempo
i=1, j=1

\phi_{1,2} =0.08 \phi_{2,0}+ 0.84\phi_{1,0}+0.08\phi_{0,0}

\phi_{1,2} =0.08 (0)+ 0.84(0.4)+0.08(5)=0.736
i=2, j=1

\phi_{2,2} =0.08 \phi_{3,1}+ 0.84\phi_{2,1}+0.08\phi_{1,1} \phi_{2,2} =0.08(0)+ 0.84(0)+0.08(0.4) = 0.032

i=3, j=1

\phi_{3,2} =0.08\phi_{4,1}+ 0.84\phi_{3,1}+0.08\phi_{2,1}=0

Para los proximos valores i>3, todos los resultados son 0

Tarea: Desarrollar la iteración 3 en el tiempo.

siguiendo las iteraciones se tiene la siguiente tabla:

[[5.0, 0.000, 0.000, 0.00000, 0.00000 , 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
 [5.0, 0.400, 0.000, 0.00000, 0.00000 , 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
 [5.0, 0.736, 0.032, 0.00000, 0.00000 , 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
 [5.0, 1.021, 0.085, 0.00256, 0.00000 , 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
 [5.0, 1.264, 0.153, 0.00901, 0.00020, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
...
]

Con lo que se obtiene la siguiente gráfica.

El resultado se interpreta mejor con una animación: (tarea)

Tarea: Presentar el orden de error de la ecuación basado en las fórmulas de diferenciación


Algorirmo en Python

# 3Eva_IT2019_T3 Difusión en sólidos
# EDP parabólicas. método explícito,usando diferencias finitas
import numpy as np
import matplotlib.pyplot as plt

# INGRESO
# Valores de frontera
Ta = 5
Tb = 0
T0 = 0
# longitud en x
a = 0
b = 0.1
# Constante K
K = 1/(1.6e-1)
# Tamaño de paso
dx = 0.02
dt = dx/100
# iteraciones en tiempo
n = 50

# PROCEDIMIENTO
# iteraciones en longitud
xi = np.arange(a,b+dx,dx)
m = len(xi)
ultimox = m-1

# Resultados en tabla u[x,t]
u = np.zeros(shape=(m,n), dtype=float)

# valores iniciales de u[:,j]
j=0
ultimot = n-1
u[0,j]= Ta
u[1:ultimox,j] = T0
u[ultimox,j] = Tb

# factores P,Q,R
lamb = dt/(K*dx**2)
P = lamb
Q = 1 - 2*lamb
R = lamb

# Calcula U para cada tiempo + dt
j = 0
while not(j>=ultimot):
    u[0,j+1] = Ta
    for i in range(1,ultimox,1):
        u[i,j+1] = P*u[i-1,j] + Q*u[i,j] + R*u[i+1,j]
    u[m-1,j+1] = Tb
    j=j+1

# SALIDA
print('Tabla de resultados')
np.set_printoptions(precision=2)
print(u)

# Gráfica
salto = int(n/10)
if (salto == 0):
    salto = 1
for j in range(0,n,salto):
    vector = u[:,j]
    plt.plot(xi,vector)
    plt.plot(xi,vector, '.r')
plt.xlabel('x[i]')
plt.ylabel('phi[i,j]')
plt.title('Solución EDP parabólica')
plt.show()

La animación se complementa con lo mostrado en la sección de Unidades.

s3Eva_IT2019_T2 Integral con interpolación

Ejercicio: 3Eva_IT2019_T2 Integral con interpolación

El ejercicio considera dos partes: interpolación e integración

a. Interpolación

Se requiere aproximar la función usando tres puntos. Para comprender la razón del método solicitado, se compara la función con dos interpolaciones:

a.1 Lagrange
a.2 Trazador cúbico sujeto

Observando la gráfica se aclara que en éste caso, una mejor aproximación se obtiene con el método  de trazador cúbico sujeto. Motivo por lo que el tema tiene un peso de 40/100 puntos

Los valores a considerar para la evaluación son:

puntos referencia xi,yi: 
[0.         0.78539816 1.57079633]
[ 0.          0.62426595 -0.97536797]
derivadas en los extremos:  
    3.141592653589793 
    0.6929852019184021
Polinomio de Lagrange
-1.80262534301178*x**2 + 2.21061873102778*x
Trazadores cúbicos sujetos
[0.         0.78539816]
-0.548171611756137*x**3 - 2.55744517923506*x**2 + 3.14159265358979*x

[0.78539816 1.57079633]
4.66299098804068*x**3 - 14.8359577843727*x**2 + 12.7851139029174*x - 2.52466795930204

------------------
Valores calculados para Trazadores cúbicos sujetos:
Matriz A: 
[[-0.26179939 -0.13089969  0.        ]
 [ 0.78539816  3.14159265  0.78539816]
 [ 0.          0.13089969  0.26179939]]
Vector B: 
[  2.34675256 -16.9893436    2.72970237]
coeficientes S: 
[-5.11489036 -7.69808822 14.27573913]
coeficientes a,b,c,d
[-0.54817161  4.66299099]
[-2.55744518 -3.84904411]
[ 3.14159265 -1.89005227]
[0.         0.62426595]

b. Integración

Como forma de comparacíon de resultados, se requiere integrar con varios métodos para comparar resultados y errores.

b.1 Integración con Cuadratura de Gauss, usando el resultado de trazador cúbico.

Se integra en cada tramo de cada polinomio:

Trazadores cúbicos sujetos
[0.         0.78539816]
-0.548171611756137*x**3 - 2.55744517923506*x**2 + 3.14159265358979*x

Se obtienen los puntos del método de cuadratura desplazados en el rango:

xa:  0.16597416116944688
xb:  0.6194240022280014
area:  0.5037962958529855

Para el segundo tramo:

[0.78539816 1.57079633]
4.66299098804068*x**3 - 14.8359577843727*x**2 + 12.7851139029174*x - 2.52466795930204
xa:  0.9513723245668951
xb:  1.4048221656254496
area:  -0.2706563884589365

Con lo que el integral total es:

Integral total:  0.23313990739404894

b.2 Integración analítica

\int_0^{\pi /2}sin(\pi x) dx

u = πx
du/dx = π
dx = du/π

se convierte en:

\frac{1}{\pi}\int sin(u) du \frac{1}{\pi}(-cos(u))

volviendo a la variable x:

\frac{1}{\pi}(-cos(\pi x)) \Big\rvert_{0}^{\frac{\pi}{2}} -\frac{1}{\pi}(cos(\pi \frac{\pi}{2})-cos(\pi(0))) = 0.24809580527879377

c. Estimación del error

Se restan los resultados de las secciones b.1 y b.2

error = |0.24809580527879377 – 0.23313990739404894 |

error = 0.014955897884744829


Algoritmo en Python

separado por literales

# 3Eva I T 2019 Interpola e Integra
import numpy as np
import sympy as sym
import matplotlib.pyplot as plt

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

def traza3sujeto(xi,yi,u,v):
    '''
    Trazador cúbico sujeto, splines
    resultado: polinomio en forma simbólica
    '''
    n = len(xi)
    # Valores h
    h = np.zeros(n-1, dtype=float)
    # Sistema de ecuaciones
    A = np.zeros(shape=(n,n), dtype=float)
    B = np.zeros(n, dtype=float)
    S = np.zeros(n-1, dtype=float)
    # 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)
    
    polinomios=[]
    
    if (n>=3):
        for i in range(0,n-1,1):
            h[i]=xi[i+1]-xi[i]
        A[0,0] = -h[0]/3
        A[0,1] = -h[0]/6
        B[0] = u-(yi[1]-yi[0])/h[0]
        for i in range(1,n-1,1):
            A[i,i-1] = h[i-1]
            A[i,i] = 2*(h[i-1]+h[i])
            A[i,i+1] = h[i]
            B[i] = 6*((yi[i+1]-yi[i])/h[i] - (yi[i]-yi[i-1])/h[i-1])
        A[n-1,n-2] = h[n-2]/6
        A[n-1,n-1] = h[n-2]/3
        B[n-1] = v-(yi[n-1]-yi[n-2])/h[n-2]

        # Resolver sistema de ecuaciones
        S = np.linalg.solve(A,B)

        # Coeficientes
        for i in range(0,n-1,1):
            a[i]=(S[i+1]-S[i])/(6*h[i])
            b[i]=S[i]/2
            c[i]=(yi[i+1]-yi[i])/h[i]-(2*h[i]*S[i]+h[i]*S[i+1])/6
            d[i]=yi[i]
      
        # polinomio en forma simbólica
        x=sym.Symbol('x')
        polinomios=[]
        for i in range(0,n-1,1):
            ptramo = a[i]*(x-xi[i])**3 + b[i]*(x-xi[i])**2 + c[i]*(x-xi[i])+ d[i]
            ptramo = ptramo.expand()
            polinomios.append(ptramo)
        parametros = [A,B,S,a,b,c,d]                                                           
    return(polinomios, parametros)

# INGRESO
f = lambda x: np.sin(np.pi*x)
muestrasf = 20
a = 0
b = np.pi/2
# Derivadas en los extremos
u = np.pi*np.cos(np.pi*a)
v = np.pi*np.cos(np.pi*b)
muestras = 3

# literal a
# PROCEDIMIENTO
xif = np.linspace(a,b,muestrasf)
yif = f(xif)

xi = np.linspace(a,b,muestras)
yi = f(xi)

# Usando Lagrange
x = sym.Symbol('x')
pL = interpola_lagrange(xi,yi)
pxL = sym.lambdify(x,pL)
pxiL =  pxL(xif)

# Trazador cúbico sujeto
pS, parametros = traza3sujeto(xi,yi,u,v)
pxiS = np.zeros(muestrasf,dtype=float)

# Evalua trazadores cúbicos sujetos
i=0
ap = xi[i]
bp = xi[i+1]
poli = sym.lambdify(x, pS[i])
for j in range(0,muestrasf,1):
    punto = xif[j]
    if (punto>bp):
        i = i+1
        ap = xi[i]
        bp = xi[i+1]
        poli = sym.lambdify(x,pS[i])
    pxiS[j] = poli(punto)

# SALIDA
print('puntos referencia xi,yi: ')
print(xi)
print(yi)
print('derivadas en los extremos: ',u,v)
print('Polinomio de Lagrange')
print(pL)
print('Trazadores cúbicos sujetos')
n = len(xi)
for i in range(0,n-1,1):
    print(xi[i:i+2])
    print(pS[i])
# Parametros de Trazadores cúbicos sujetos
print('Matriz A: ')
print(parametros[0])
print('Vector B: ')
print(parametros[1])
print('coeficientes S: ')
print(parametros[2])
print('coeficienetes a,b,c,d')
print(parametros[3])
print(parametros[4])
print(parametros[5])
print(parametros[6])

# Gráficas
plt.plot(xif,yif, label='funcion')
plt.plot(xi,yi,'o', label='muestras')
plt.plot(xif,pxiL, label='p(x)_Lagrange')
plt.plot(xif,pxiS, label='p(x)_Traza3Sujeto')
plt.legend()
plt.xlabel('x')
plt.show()

# literal b
# cuadratura de Gauss de dos puntos
def integraCuadGauss2p(funcionx,a,b):
    x0 = -1/np.sqrt(3)
    x1 = -x0
    xa = (b+a)/2 + (b-a)/2*(x0)
    xb = (b+a)/2 + (b-a)/2*(x1)
    area = ((b-a)/2)*(funcionx(xa) + funcionx(xb))
    print('xa: ',xa)
    print('xb: ',xb)
    print('area: ', area)
    return(area)

# INGRESO
f0 = sym.lambdify(x,pS[0])
f1 = sym.lambdify(x,pS[1])
# Procedimiento
I0 = integraCuadGauss2p(f0,xi[0],xi[1])
I1 = integraCuadGauss2p(f1,xi[1],xi[2])
It = I0+I1

# SALIDA
print('Integral 1: ', I0)
print('Integral 2: ', I1)
print('Integral total: ',It)

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

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

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

s3Eva_IT2017_T3 Sustancia en lago

Ejercicio: 3Eva_IT2017_T3 Sustancia en lago

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

Sección transversal

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

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

repidiendo la fórmula entre cada par de puntos consecutivos

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

Concentración promedio

Para los integrales usamos la regla del trapecio:

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

Aplicando los algoritmos en Python para todos los puntos:

# 3Eva_IT2017_T3 Sustancia en lago
import numpy as np
import matplotlib.pyplot as plt

# INGRESO
zi = np.array([0.  , 4   , 8   , 12    , 16])
vi = np.array([9.82, 5.11, 1.96,  0.393,  0.])
ci = np.array([10.2, 8.5 , 7.4 ,  5.2  ,  4.1])

# PROCEDIMIENTO
n = len(zi)
# primera derivada hacia adelante con error O(h)
dv = np.zeros(n,dtype=float)
for i in range(0,n-1,1):
    h = zi[i+1]-zi[i]
    dv[i]=(vi[i+1]-vi[i])/h

As = -dv*zi

# integrales por rectángulo
numerador = 0
for i in range(0,n-1,1):
    altura = (ci[i]*As[i]+ci[i+1]*As[i+1])/2
    numerador = numerador +altura*(zi[i+1]-zi[i])

denominador = 0
for i in range(0,n-1,1):
    altura = (As[i]+As[i+1])/2
    denominador = denominador +altura*(zi[i+1]-zi[i])

cpromedio = numerador/denominador

# SALIDA
print('dv/dz: ')
print(dv)
print('numerador: ',numerador)
print('denominador: ',denominador)
print('concentracion promedio: ',cpromedio)

# Grafica
plt.subplot(121)
plt.plot(zi,vi)
plt.plot(zi,vi,'bo')
plt.xlabel('profundidad z')
plt.ylabel('Volumen')
plt.grid()
plt.subplot(122)
plt.plot(zi,ci, color = 'orange')
plt.plot(zi,ci,'ro')
plt.xlabel('profundidad z')
plt.ylabel('concentración')
plt.grid()
plt.show()

s3Eva_IT2017_T4 EDP elíptica, placa desplazada

Ejercicio: 3Eva_IT2017_T4 EDP elíptica, placa desplazada

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

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

1 <  x < 2
1 <  y < 2

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


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

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

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

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


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

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


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

Iteraciones:

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

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

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

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

i=1, j=1


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

i = 2, j =1


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

Método iterativo

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

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

Algoritmo en Python

# 3Eva_IT2017_T4 EDP elíptica, placa desplazada
# método iterativo
import numpy as np

# INGRESO
# longitud en x
a = 1
b = 2
# longitud en y
c = 1
d = 2
# tamaño de paso
dx = 0.25
dy = 0.25
# funciones en los bordes de la placa
abajo     = lambda x,y: x*np.log(x)
arriba    = lambda x,y: x*np.log(4*(x**2))
izquierda = lambda x,y: y*np.log(y)
derecha   = lambda x,y: 2*y*np.log(2*y)
# función de la ecuación
fxy = lambda x,y: x/y + y/x

# control de iteraciones
maxitera = 100
tolera = 0.0001

# PROCEDIMIENTO
# tamaño de la matriz
n = int((b-a)/dx)+1
m = int((d-c)/dy)+1
# vectores con valore de ejes
xi = np.linspace(a,b,n)
yj = np.linspace(c,d,m)
# matriz de puntos muestra
u = np.zeros(shape=(n,m),dtype=float)

# valores en los bordes
u[:,0]   = abajo(xi,yj[0])
u[:,m-1] = arriba(xi,yj[m-1])
u[0,:]   = izquierda(xi[0],yj)
u[n-1,:] = derecha(xi[n-1],yj)

# valores interiores
# para menos iteraciones
mitadx = int(n/2)
mitady = int(m/2)
promedio = (u[mitadx,0]+u[mitadx,m-1]+u[0,mitady]+u[n-1,mitady])/4
u[1:n-1,1:m-1] = promedio

# método iterativo
itera = 0
converge = 0
while not(itera>=maxitera or converge==1):
    itera = itera +1
    # copia u para calcular errores entre iteraciones
    nueva = np.copy(u)
    for i in range(1,n-1):
        for j in range(1,m-1):
            # usar fórmula desarrollada para algoritmo
            u[i,j] = (u[i-1,j]+u[i+1,j]+u[i,j-1]+u[i,j+1]-(dy**2)*fxy(xi[i],yj[j]))/4 
    diferencia = nueva-u
    erroru = np.linalg.norm(np.abs(diferencia))
    if (erroru<tolera):
        converge=1

# SALIDA
print('iteraciones: ',itera)
print('error entre iteraciones: ',erroru)
print('solución para u: ')
print(u)

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

# matrices de ejes para la gráfica 3D
X, Y = np.meshgrid(xi, yj)
U = np.transpose(u) # ajuste de índices fila es x
figura = plt.figure()

grafica = Axes3D(figura)
grafica.plot_surface(X, Y, U, rstride=1, cstride=1, cmap=cm.Reds)

plt.title('EDP elíptica')
plt.xlabel('x')
plt.ylabel('y')
plt.show()