3Eva_IT2019_T3 Difusión en sólidos

3ra Evaluación I Término 2019-2020. 10/Septiembre/2019. MATG1013

Tema 3. (30 Puntos).  En el año 1855, los experimentos de Adolf Fick tratan sobre la medición de concentraciones y sus flujos, también ahora aplicados a la difusión en sólidos que en ese tiempo no se consideraba posible.

La gráfica muestra los cambios en el tiempo de concentración Φ de un gas en un sólido (estado no-estacionario) para un sólido semi infinito (eje y).

La segunda ley de Fick predice la forma en que la difusión causa que la concentración cambie con el tiempo. Se trata de una ecuación diferencial parcial que en una dimensión se escribe:

\frac{\partial \phi}{\partial t} = D\frac{\partial ^2 \phi}{\partial x^2}

Φ(0, t) = 5; Φ(L, t) = 0; Φ(x,0) = 0; D = 0.16; L =0.1

a. Plantee las ecuaciones, la malla, desarrolle y obtenga el modelo Φ(xi,tj)
b. Aproxime la solución con Δx = 0.02, Δt = Δx/100. Realice al menos tres iteraciones en el eje tiempo.
c. Estime el error de Φ(xi,tj)

Rúbrica: Construir la malla (5 puntos), plantear la ecuación en el nodo i,j (5 puntos), modelo de ecuación (5 puntos), literal b (10 puntos), literal c (5 puntos).

Referencia: https://es.wikipedia.org/wiki/Leyes_de_Fick;
Difusión 2ª Ley de Fick|7/22|UPV (2011) https://www.youtube.com/watch?v=HHBvZDNvTic

s3Eva_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: Presentar el orden de error de la ecuación basado en las fórmulas de diferenciación


Algorirmo en Python

# 3EIT2019T4.Difusión en sólidos. 2da Ley de Fick
# EDP parabólicas. método explícito, usando diferencias finitas
# http://blog.espol.edu.ec/matg1013/3eva_it2018_t3-edp-parabolica-temperatura-en-varilla/
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.

3Eva_IT2019_T2 Integral con interpolación

3ra Evaluación I Término 2019-2020. 10/Septiembre/2019. MATG1013

Tema 2. (40 Puntos) Construya un polinomio que aproxime a f(x) = sin(\pi x) usando los puntos x=0, π/4, π/2 y aproxime la integral de 0 a π/2.

a. Realice la interpolación mediante el método de trazador cúbico fijo
b. Integre usando el método de cuadratura de Gauss
c. Estime el error para el ejercicio.

Rúbrica: Bosquejo de gráficas (5 puntos), literal a, planteo de fórmulas (5 puntos), calcula los parámetros (10 puntos), literal b (15 puntos), literal c (5 puntos).

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

3Eva_IT2019_T1 Ecuaciones simultáneas

3ra Evaluación I Término 2019-2020. 10/Septiembre/2019. MATG1013

Tema 1. (30 Puntos).  Deteremine las raíces de las ecuaciones simultáneas siguientes:

y = -x^2 +x + 0.75 y+5xy=x^3

a. Realice un bosquejo para cada ecuación
b. Use el método de Newton-Raphson con x0=1 , y0=0.75, realice 3 iteraciones
c. Estime el orden del error

Rúbrica: literal a (5 puntos), literal b planteo (5 puntos), iteraciones (15 puntos), literal c (5 puntos)

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

2Eva_IT2019_T3 EDP Elíptica Placa 6×5

2da Evaluación I Término 2019-2020. 27/Agosto/2019. MATG1013

Tema 3. (30 Puntos) Una placa rectangular de plata de 6×5 cm tiene calor que se genera uniformemente en todos los puntos, con una rapidez q = 1.5 cal/cm3 s.
Al representar con x la distancia a lo largo del borde de longitud 6 cm y con y la de 5 cm.

Suponga que la temperatura en los bordes se mantiene como se indica:

u(x,0) = x(6-x) u(x,5)=0 0≤x≤6
u(0,y) = y(5-y) u(6,y)=0 0≤y≤5

Donde el origen se encuentra en una esquina de la placa y los bordes se hayan a lo largo de los ejes positivos x, y.

La temperatura de estado estable u(x,y) satisface la ecuación de Poisson:

\frac{\partial^2 u}{\partial x^2} (x,y)+\frac{\partial ^2 u}{\partial y^2 } (x,y) = -\frac{q}{K}

0≤x≤6
0≤y≤5

Donde K, la conductividad térmica es 1.04 cal/cm deg s.

a. Aproxime la temperatura u(x,y) en los nodos de la malla con hx =2, hy= 2.5

b. Exprese el término del error

Rúbrica: literal a expresiones (10 puntos), valor (5 puntos), literal b (5 puntos)


Referencia: Ejercicio 12.1.8, Burden 9Ed, p724.

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

2Eva_IT2019_T2 Péndulo vertical

2da Evaluación I Término 2019-2020. 27/Agosto/2019. MATG1013

Tema 2. (40 Puntos) Suponga que un péndulo tiene 0.6 m de Longitud, se desplaza θ desde la posición vertical de equilibrio.

\frac{d^2\theta }{dt^2}+\frac{g}{L}\sin (\theta)=0 0\lt t \lt 1 g = 9.81 \frac{m}{s^2} \theta(0) = \frac{\pi}{6} \theta '(0) = 0

a. Aproxime la solución de la ecuación para t = [0,1] con pasos de h=0.2
b. Aproxime el valor del error

Rúbrica: literal a, expresiones (20 puntos), valor (10 puntos), literal b (10 puntos)


Referencia: Ejercicio 5.9.8, Burden 9Ed, p338.
2Eva_IT2010_T2 Movimiento angular

Professor of Physics Emeritus Walter Lewin.  Lec 11 | 8.01 Physics I: Classical Mechanics, Fall 1999.

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