Categoría: Sol_3Eva 2007-2010

Ejercicio: 3Eva_IIT2010_T4 EDO con Taylor

$\frac{\delta y}{\delta x} = \frac{y^3}{1-2xy^2}$ $y(0) = 1, 0 \leq x \leq 1$

escrita en forma simplificada

$y' = \frac{y^3}{1-2xy^2}$

tomando como referencia Taylor de 2 términos más el término de error O(h²)

$y_{i+1} = y_i +\frac{h}{1!}y'_i + \frac{h^2}{2!}y''_i$

Se usa hasta el segundo término para el algoritmo.

$y_{i+1} = y_i +\frac{h}{1!}$

Con lo que se puede realizar el algoritmo

estimado[xi,yi]
[[0.  1. ]
 [0.2 1.2 ]
 [0.4 2.01509434]
 [0.6 1.28727044]
 [0.8 0.85567954]
 [1.  0.12504631]]

Algoritmo en Python

# 3Eva_IIT2010_T4 EDO con Taylor
# estima la solución para muestras espaciadas h en eje x
# valores iniciales x0,y0
# entrega arreglo [[x,y]]
import numpy as np

def edo_taylor2t(d1y,x0,y0,h,muestras):
    tamano = muestras + 1
    estimado = np.zeros(shape=(tamano,2),dtype=float)
    # incluye el punto [x0,y0]
    estimado[0] = [x0,y0]
    x = x0
    y = y0
    for i in range(1,tamano,1):
        y = y + h*d1y(x,y)
        x = x+h
        estimado[i] = [x,y]
    return(estimado)

# PROGRAMA PRUEBA
# 3Eva_IIT2010_T4 EDO con Taylor

# INGRESO.
# d1y = y' = f, d2y = y'' = f'
d1y = lambda x,y: (y**3)/(1-2*x*(y**2))
x0 = 0
y0 = 1
h = 0.2
muestras = 5

# PROCEDIMIENTO
puntos = edo_taylor2t(d1y,x0,y0,h,muestras)
xi = puntos[:,0]
yi = puntos[:,1]

# SALIDA
print('estimado[xi,yi]')
print(puntos)

# Gráfica
import matplotlib.pyplot as plt
plt.plot(xi[0],yi[0],'o', color='r', label ='[x0,y0]')
plt.plot(xi[1:],yi[1:],'o', color='g', label ='y estimada')
plt.title('EDO: Solución con Taylor 2 términos')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.grid()
plt.show()

Nota: Revisar los resultados no lineales con los valores de h=0.02

Ejercicio: 3Eva_IT2010_T2 EDO problema con valor inicial dy/dx

Ecuación del problema:

$(2y^2 + 4x^2)\delta x -xy \delta y =0$

se despeja dy/dx:

$(2y^2 + 4x^2)\delta x = xy \delta y$ $\frac{\delta y}{\delta x} =\frac{2y^2 + 4x^2}{xy}$

con valores iniciales de x₀ = 1, y₀=-2 , h=0.2 e intervalo [1,2]

Usando Runge-Kutta de 2do Orden

Iteración 1

$K_1= h\frac{\delta y}{\delta x}(1,-2)$ $=(0.2)\frac{2(-2)^2 + 4(1)^2}{(1)(-2)}= -1.2$ $K_2= h\frac{\delta y}{\delta x}(1+0.2,-2+(-1.2))$ $=(0.2)\frac{2(-2-1.2)^2 + 4(1+0.2)^2}{(1+0.2)(-2-1.2)}= -1.3667$ $y_1 = -2 + \frac{-1.2-1.3667}{2} = -3.2833$ $x_1 = x_0 + h = 1 + 0.2 = 1.2$ $error = O(h^3) = O(0.2^3) = 0.008$

Iteración 2

$K_1= h\frac{\delta y}{\delta x}(1.2,-3.2833)$ $=(0.2)\frac{2(-3.2833)^2 + 4(1.2)^2}{(1.2)(-3.2833)}= -1.3868$ $K_2= h\frac{\delta y}{\delta x}(1.2+0.2,-3.2833+(-1.3868))$ $=(0.2)\frac{2(-3.2833+(-1.3868))^2 + 4(1.2+0.2)^2}{(1.2+0.2)(-3.2833+(-1.3868))}= -1.5742$ $y_2 = -3.2833 + \frac{-1.3868-1.5742}{2} = -4.7638$ $x_2 = x_1 + h = 1.2 + 0.2 = 1.4$

Iteración 3

$K_1= h\frac{\delta y}{\delta x}(1.4,-4.7638)$ $=(0.2)\frac{2(-4.76383)^2 + 4(1.4)^2}{(1.4)(-4.7638)}= -1.5962$ $K_2= h\frac{\delta y}{\delta x}(1.4+0.2,-4.7638+(-1.5962))$ $=(0.2)\frac{2(-4.7638+(-1.5962))^2 + 4(1.4+0.2)^2}{(1.4+0.2)(-4.7638+(-1.5962))}= -1.7913$ $y_3 = -4.7638 + \frac{-1.5962-1.7913}{2} = -6.4576$ $x_3 = x_2 + h = 1.4 + 0.2 = 1.6$

con lo que usando el algoritmo se obtiene la tabla y gráfica:

 [xi,       yi,       K1,     K2     ]
[[  1.      -2.       0.       0.    ]
 [  1.2     -3.2833  -1.2     -1.3667]
 [  1.4     -4.7638  -1.3868  -1.5742]
 [  1.6     -6.4576  -1.5962  -1.7913]
 [  1.8     -8.3698  -1.8126  -2.0119]
 [  2.     -10.5029  -2.032   -2.2342]

Algoritmo en Python

# 3Eva_IT2010_T2 EDO
import numpy as np
import matplotlib.pyplot as plt

# INGRESO
d1y = lambda x,y : (2*(y**2)+4*(x**2))/(x*y)
x0 = 1
y0 = -2
h = 0.2
a = 1
b = 2

# PROCEDIMIENTO
muestras = int((b -a)/h)+1
tabla = np.zeros(shape=(muestras,4),dtype=float)

i = 0
xi = x0
yi = y0
tabla[i,:] = [xi,yi,0,0]

i = i+1
while not(i>=muestras):
    K1 = h* d1y(xi,yi)
    K2 = h* d1y(xi+h,yi+K1)
    yi = yi + (K1+K2)/2
    xi = xi +h
    tabla[i,:] = [xi,yi,K1,K2]
    i = i+1
# vector para gráfica
xg = tabla[:,0]
yg = tabla[:,1]

# SALIDA
# muestra 4 decimales
np.set_printoptions(precision=4)
print(' [xi, yi, K1, K2]')
print(tabla)

# Gráfica
plt.plot(xg,yg)
plt.xlabel('xi')
plt.ylabel('yi')
plt.grid()
plt.show()

Tarea: Realizar con Runge Kutta de 4to Orden

Ejercicio: 3Eva_IT2010_T1 Envase cilíndrico

Se conoce que el volumen del cilindro es 1000 cm³
que se calcula como:

$Volumen = area_{base} . altura = \pi r^2 h$ $\pi r^2 h = 1000$ $h = \frac{1000}{\pi r^2}$

con lo que la altura queda en función del radio, pues el volumen es una constante.

Para conocer el área de la hoja de material para el envase se conoce que las tapas cilíndricas deben ser 0.25 mas que el radio para sellar el recipiente

Area de una tapa circular:

$Area_{tapa} = \pi (r+ 0.25)^2$

para la hoja lateral:

$Area_{lateral} = base * altura = (2 \pi r + 0.25) h$

que en función del radio, usando la fórmula para h(r) se convierte en:

$Area_{lateral} = (2 \pi r + 0.25) \frac{1000}{\pi r^2}$

por lo que el total de material de hoja a usar corresponde a dos tapas circulares y una hoja lateral.

$Area total = 2 \pi (r+ 0.25)^2 + (2 \pi r + 0.25) \frac{1000}{\pi r^2}$

---

la gráfica permite observar el comportamiento entre el área de la hoja necesaria para fabricar el envase y el radio de las tapas circulares.

El objetivo es encontrar el área mínima para consumir menos material. Con la fórmula mostrada se puede aplicar uno de los métodos para encontrar raíces a partir de la derivada de la fórmula.

Tarea: Encontrar el valor requerido para r con la tolerancia indicada para el examen.

---

Instrucciones para la gráfica:

# 3Eva_IT2010_T1 Envase cilíndrico
import numpy as np
import matplotlib.pyplot as plt

# INGRESO
a = 0.5
b = 20
muestras = 51

Area = lambda r: (2*np.pi)*(r+0.25)**2 +(2*np.pi*r+0.25)*1000/(np.pi*(r**2))

# PROCEDIMIENTO
ri = np.linspace(a,b,muestras)
Areai = Area(ri)

# SALIDA
plt.plot(ri,Areai)
plt.xlabel('r (cm)')
plt.ylabel('Area (cm2)')
plt.title('Area hoja para material cilindro')
plt.show()

Ejercicio: 3Eva_IIT2009_T2 Valor inicial Runge-Kutta 4to orden dy/dx

La ecuación del problema:

$(1-x^2)y' - xy = x (1-x^2)$

se despeja:

$(1-x^2)y' = x (1-x^2) - xy$ $y' = x - \frac{xy}{(1-x^2)}$

con valores iniciales x₀ = 0 , y₀ = 2, h = 0.1 en el intervalo [0, 0.5]

Usando Runge-Kutta de 2do Orden

Iteración 1

$K_1 = h y'(0,2)$ $= (0.1)\Big[0 - \frac{0 (2)}{(1-(0^2))}\Big] = 0$ $K_2 = h y'(0+0.1, 2 + 0)$ $= (0.1)\Big[0.1 - \frac{0.1 (2+0)}{(1-(0.1^2))}\Big] = 0.0302$ $y_1= 2 + \frac{0+0.0302}{2} = 2.0151$ $x_1 = x_0 + h = 0 + 0.1 = 0.1$

Iteración 2

$K_1 = h y'(0.1,2.0151)$ $= (0.1)\Big[0.1 - \frac{0.1 (2.0151)}{(1-(0.1^2))}\Big] = 0.0304$ $K_2 = h y'(0.1+0.1, 2.0151 + 0.0304)$ $= (0.1)\Big[0.2 - \frac{0.2 (2.0151 + 0.0304)}{(1-(0.2^2))}\Big] = 0.0626$ $y_2 = 2.0151 + \frac{0.0304+0.0626}{2} = 2.0616$ $x_2 = x_1 + h = 0.1 + 0.1 = 0.2$

Iteración 3

$K_1 = h y'(0.2,2.0616)$ $= (0.1)\Big[0.2 - \frac{0.2 (2.0616)}{(1-(0.2^2))}\Big] = 0.0629$ $K_2 = h y'(0.2+0.1, 2.0616 + 0.0629)$ $= (0.1)\Big[0.3 - \frac{0.3 (2.0616 + 0.0629)}{(1-(0.3^2))}\Big] = 0.1$ $y_3 = 2.0151 + \frac{0.0629+0.1}{2} = 2.1431$ $x_3 = x_2 + h = 0.2 + 0.1 = 0.3$

siguiendo el algoritmo se completa la tabla:

 [xi,    yi,    K1,    K2    ]
[[0.     2.     0.     0.    ]
 [0.1    2.0151 0.     0.0302]
 [0.2    2.0616 0.0304 0.0626]
 [0.3    2.1431 0.0629 0.1   ]
 [0.4    2.2668 0.1007 0.1468]
 [0.5    2.4463 0.1479 0.211 ]]

Algoritmo en Python

# 3Eva_IIT2009_T2 Valor inicial Runge-Kutta 4to orden dy/dx
import numpy as np
import matplotlib.pyplot as plt

# INGRESO
d1y = lambda x,y : x + x*y/(1-x**2)
x0 = 0
y0 = 2
h = 0.1
a = 0
b = 1/2

# PROCEDIMIENTO
muestras = int((b -a)/h)+1
tabla = np.zeros(shape=(muestras,4),dtype=float)

i = 0
xi = x0
yi = y0
tabla[i,:] = [xi,yi,0,0]

i = i+1
while not(i>=muestras):
    K1 = h* d1y(xi,yi)
    K2 = h* d1y(xi+h,yi+K1)
    yi = yi + (K1+K2)/2
    xi = xi +h
    tabla[i,:] = [xi,yi,K1,K2]
    i = i+1
# vector para gráfica
xg = tabla[:,0]
yg = tabla[:,1]

# SALIDA
# muestra 4 decimales
np.set_printoptions(precision=4)
print(' [xi, yi, K1, K2]')
print(tabla)

# Gráfica
plt.plot(xg,yg)
plt.xlabel('xi')
plt.ylabel('yi')
plt.grid()
plt.show()

Tarea: Realizar con Runge-Kutta de 4to Orden

Ejercicio: 3Eva_IT2009_T2 EDO Taylor Seno(x)

La ecuación del problema es:

$xy'+ 2y = \sin (x)$ $\frac{\pi}{2} \leq x \leq \frac{3\pi}{2}$ $y\Big(\frac{\pi}{2} \Big) = 1$

Para el algoritmo se escribe separando y’:

$y' = \frac{\sin (x)}{x} - 2\frac{y}{x}$

tomando como referencia Taylor de 3 términos más el término de error O(h³)

$y_{i+1} = y_i +\frac{h}{1!}y'_i + \frac{h^2}{2!}y''_i + \frac{h^3}{3!}y'''_i$

Se usa hasta el tercer término para el algoritmo.

$y_{i+1} = y_i +\frac{h}{1!}y'_i + \frac{h^2}{2!}y''_i$

Se determina que se requiere la segunda derivada para completar la aproximación. A partir de la ecuación del problema se aplica en cada término:

$\Big(\frac{u}{v}\Big)' = \frac{u'v-uv' }{v^2}$ $y'' = \frac{\cos (x) x - sin(x)}{x^2} - 2\frac{y'x-y}{x^2}$ $y'' = \frac{\cos (x) x - sin(x)}{x^2} - 2\frac{\Big(\frac{\sin (x)}{x} - 2\frac{y}{x} \Big)x-y}{x^2}$

Con lo que se realizan las iteraciones para llenar la tabla

iteración 1

x₀ = π/2= 1.57079633

y₀= 1

$y' = \frac{\sin (π/2)}{π/2} - 2\frac{1}{π/2} = -0.40793719$ $y'' = \frac{\cos (π/2) π/2 - sin(π/2)}{(π/2)^2}+$ $- 2\frac{(-0.40793719)(π/2)-1}{(π/2)^2}= 0.48531359$ $y_{1} = 1 +\frac{\pi/10}{1!}(-0.40793719) + \frac{(\pi/10)^2}{2!}(0.48531359) = 0.86$

x₁ = x₀+h = 1.57079633 + π/10 =1.88495559

iteración 2

x₁ = 1.88495559

y₁ = 0.86

$y' = \frac{\sin (1.88495559)}{1.88495559} - 2\frac{0.86}{1.88495559} = -0.31947719$ $y'' = \frac{\cos (1.88495559)(1.88495559) - sin(1.88495559)}{(1.88495559)^2}$ $- 2\frac{(-0.31947719) (1.88495559)-y}{(1.88495559)^2} = 0.16854341$ $y_{2} = 0.86 +\frac{\pi /10}{1!}(-0.31947719) + \frac{(\pi /10)^2}{2!}(0.16854341) = 0.75579202$

x₂ = x₁+ h = 1.88495559 + π/10 = 2.19911486

iteración 3

x₂ = 2.19911486

y₂ = 0.75579202

$y' = \frac{\sin (2.19911486)}{2.19911486} - 2\frac{y}{2.19911486} = -0.29431724$ $y'' = \frac{\cos (2.19911486)(2.19911486) - sin(2.19911486)}{(2.19911486)^2} +$ $- 2\frac{(-0.29431724)(2.19911486)-0.75579202}{(2.19911486)^2} = 0.0294177$ $y_{3} = 0.75579202 +\frac{\pi/10}{1!}(-0.29431724) + \frac{(\pi /10)^2}{2!}(0.0294177)$

x₃ = x₃+h = 2.19911486 + π/10 = 2.19911486

Con lo que se puede realizar el algoritmo, obteniendo la siguiente tabla:

estimado
 [xi,          yi,         d1y,        d2y       ]
[[ 1.57079633  1.         -0.40793719  0.48531359]
 [ 1.88495559  0.86       -0.31947719  0.16854341]
 [ 2.19911486  0.75579202 -0.29431724  0.0294177 ]
 [ 2.51327412  0.66374258 -0.29583247 -0.02247944]
 [ 2.82743339  0.5727318  -0.30473968 -0.02730493]
 [ 3.14159265  0.47868397 -0.31027009 -0.00585871]
 [ 3.45575192  0.38159973 -0.30649476  0.02930236]
 [ 3.76991118  0.28383639 -0.29064275  0.06957348]
 [ 4.08407045  0.18899424 -0.26221809  0.10859762]
 [ 4.39822972  0.10111944 -0.22243506  0.14160656]
 [ 4.71238898  0.02410027  0.          0.        ]]

Algoritmo en Python

# 3Eva_IT2009_T2 EDO Taylor Seno(x)
# EDO. Método de Taylor 3 términos 
# estima solucion para muestras separadas h en eje x
# valores iniciales x0,y0
# entrega arreglo [[x,y]]
import numpy as np

def edo_taylor3t(d1y,d2y,x0,y0,h,muestras):
    tamano = muestras + 1
    estimado = np.zeros(shape=(tamano,4),dtype=float)
    # incluye el punto [x0,y0]
    estimado[0] = [x0,y0,0,0]
    x = x0
    y = y0
    for i in range(1,tamano,1):
        y = y + h*d1y(x,y) + ((h**2)/2)*d2y(x,y)
        x = x+h
        estimado[i-1,2:]= [d1y(x,y),d2y(x,y)]
        estimado[i,0:2] = [x,y]
    return(estimado)

# PROGRAMA PRUEBA
# 3Eva_IT2009_T2 EDO Taylor Seno(x).

# INGRESO.
# d1y = y', d2y = y''
d1y = lambda x,y: np.sin(x)/x - 2*y/x
d2y = lambda x,y: (x*np.cos(x)-np.sin(x))/(x**2)-2*(x*(np.sin(x)/x - 2*y/x)-y)/(x**2)
x0 = np.pi/2
y0 = 1
h = np.pi/10
muestras = 10

# PROCEDIMIENTO
puntos = edo_taylor3t(d1y,d2y,x0,y0,h,muestras)
xi = puntos[:,0]
yi = puntos[:,1]

# SALIDA
print('estimado[xi, yi, d1yi, d2yi]')
print(puntos)

# Gráfica
import matplotlib.pyplot as plt
plt.plot(xi[0],yi[0],'o', color='r', label ='[x0,y0]')
plt.plot(xi[1:],yi[1:],'o', color='g', label ='y estimada')
plt.title('EDO: Solución con Taylor 3 términos')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.grid()
plt.show()

Ejercicio: 3Eva_IIT2007_T3 EDO Taylor orden 2

La ecuación del problema es:

$y'= 1 +\frac{y}{t} + \Big(\frac{y}{t}\Big) ^2$ $1\leq t\leq 2$ $y(1)=0, h=0.2$

Tomando como referencia Taylor de 3 términos más el término de error O(h³)

$y_{i+1} = y_i +\frac{h}{1!}y'_i + \frac{h^2}{2!}y''_i + \frac{h^3}{3!}y'''_i$

Se usa hasta el tercer término para el algoritmo.

$y_{i+1} = y_i +\frac{h}{1!}y'_i + \frac{h^2}{2!}y''_i$

Se determina que se requiere la segunda derivada para completar la aproximación. A partir de la ecuación del problema se aplica en cada término:

$\Big(\frac{u}{v}\Big)' = \frac{u'v-uv' }{v^2}$ $y'= 1 +\frac{y}{t} + \frac{y^2}{t^2}$ $y''= 0+\frac{y't-y}{t^2} + \frac{2yy't^2-2y^2t}{t^4}$ $y''= \frac{y't-y}{t^2} + \frac{2y\Big(1 +\frac{y}{t} + \frac{y^2}{t^2}\Big) t^2-2y^2t}{t^4}$

Con lo que se realizan las iteraciones para llenar la tabla

iteración 1

t₀ = 1

y₀= 0

$y'= 1 +\frac{0}{1} + \Big(\frac{0}{1}\Big)^2 = 1$ $y''= \frac{(1)(1)-0}{(1)^2} + \frac{2(0)(1)(1)^2-2(0)^2(1)}{(1)^4} = 1$ $y_{1} = 0 +\frac{0.2}{1!}(1) + \frac{0.2^2}{2!}(1) = 0.22$

t₁ = t₀ + h = 1 + 0.2 = 1.2

iteración 2

t₁ = 1.2

y₁= 0.22

$y'= 1 +\frac{0.22}{1.2} + \Big(\frac{0.22}{1.2}\Big) ^2 = 1.21694444$ $y''= \frac{y'(1.2)-0.22}{1.2^2} + \frac{2(0.22)y'(1.2)^2-2(0.22)^2(1.2)}{(1.2)^4} = 1.17716821$ $y_{2} = 0.22 +\frac{0.2}{1!}(1.21694444) + \frac{0.2^2}{2!}(1.17716821) = 0.48693225$

t₂ = t₁ + h = 1.2 + 0.2 = 1.4

iteración 3

t₂ = 1.4

y₂= 0.48693225

$y'= 1 +\frac{0.48693225}{1.4} + \Big(\frac{0.48693225}{1.4}\Big) ^2 = 1.46877968$ $y''= \frac{(1.46877968)(1.4)-0.48693225}{1.4^2} +$ $\frac{2(0.48693225)(1.46877968)(1.4)^2-2(0.48693225)^2(1.4)}{1.4^4} = 1.35766995$ $y_{3} = 0.48693225 +\frac{0.2}{1!}(1.46877968) + \frac{0.2^2}{2!}(1.35766995) = 0.80784159$

t₃ = t₂ + h = 1.4 + 0.2 = 1.6

Con lo que se puede realizar el algoritmo, obteniendo la siguiente tabla:

estimado
 [xi,        yi,        d1yi,      d2yi]
[[1.         0.         1.         1.        ]
 [1.2        0.22       1.21694444 1.17716821]
 [1.4        0.48693225 1.46877968 1.35766995]
 [1.6        0.80784159 1.759826   1.57634424]
 [1.8        1.19133367 2.09990017 1.8564435 ]
 [2.         1.64844258 0.         0.        ]]

Ejercicio: 3Eva_IIT2007_T1 EDP Eliptica, problema de frontera

Con los datos del ejercicio se plantean de la siguiente forma en los bordes:

La ecuación a resolver es:

$\frac{\delta^2u}{\delta x^2} +\frac{\delta^2 u}{\delta y^2} = 4$

Que en su forma discreta, con diferencias divididas centradas es:

$\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} = 4$

Al agrupar constantes se convierte en:

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

Siendo Δx= 1/3 y Δy =2/3, se mantendrá la relación λ = (Δx/Δy)²

$u[i-1,j]-2u[i,j]+u[i+1,j] +$ $+ \lambda u[i,j-1]-2\lambda u[i,j]+\lambda u[i,j+1] =$ $= 4(\Delta x)^2$

---

$u[i-1,j]-2(1+\lambda)u[i,j]+u[i+1,j] +$ $+ \lambda u[i,j-1]+\lambda u[i,j+1] = 4(\Delta x)^2$

---

Se pueden realizar las iteraciones para los nodos y reemplaza los valores de frontera:

i=1 y j=1

$u[0,1]-2(1+\lambda)u[1,1]+u[2,1] +$ $+ \lambda u[1,0]+\lambda u[1,2] = 4(\Delta x)^2$ $\Delta y^2-2(1+\lambda)u[1,1]+u[2,1] +$ $+ \Delta x^2+\lambda u[1,2] = 4(\Delta x)^2$ $-2(1+\lambda)u[1,1]+u[2,1] +\lambda u[1,2]$ $= 4(\Delta x)^2 - \Delta y^2 - \Delta x^2$ $-2(1+0.25)u[1,1]+u[2,1] +0.25 u[1,2]$ $= 4(1/3)^2 - (2/3)^2- (1/3)^2$ $-2.5u[1,1]+u[2,1] +0.25 u[1,2] = -0.1111$

i=2 , j=1

$u[1,1]-2(1+\lambda)u[2,1]+u[3,1] +$ $+ \lambda u[2,0]+\lambda u[2,2] = 4(\Delta x)^2$ $u[1,1]-2(1+\lambda)u[2,1]+(\delta y-1)^2 +$ $+ \Delta x^2 +\lambda u[2,2] = 4(\Delta x)^2$ $u[1,1]-2(1+\lambda)u[2,1] + \lambda u[2,2]$ $= 4(\Delta x)^2 - (\Delta y-1)^2 - \Delta x^2$ $u[1,1]-2(1+0.25)u[2,1] + 0.25 u[2,2]$ $= 4(1/3)^2 - (2/3-1)^2 - (1/3)^2$ $u[1,1]-2.5u[2,1] + 0.25 u[2,2] = 0.2222$

i=1 , j=2

$u[0,2]-2(1+\lambda)u[1,2]+u[2,2] +$ $+ \lambda u[1,1]+\lambda u[1,3] = 4(\Delta x)^2$ $(2\Delta y^2)-2(1+\lambda)u[1,2]+u[2,2] +$ $+ \lambda u[1,1]+\lambda (\Delta x-1)^2 = 4(\Delta x)^2$ $-2(1+\lambda)u[1,2]+u[2,2] + \lambda u[1,1]$ $= 4(\Delta x)^2 - (2\Delta y^2) -\lambda (\Delta x-1)^2$ $-2(1+0.25)u[1,2]+u[2,2] + 0.25 u[1,1]$ $= 4(1/3)^2 - (2(2/3)^2] -0.25 (1/3-1)^2$ $-2.5u[1,2]+u[2,2] + 0.25 u[1,1] = -0.5555$

i=2, j=2

$u[1,2]-2(1+\lambda)u[2,2]+u[3,2] +$ $+ \lambda u[2,1]+\lambda u[2,3] = 4(\Delta x)^2$ $u[1,2]-2(1+\lambda)u[2,2]+(\Delta y-1)^2 +$ $+ \lambda u[2,1]+\lambda (\Delta x-1)^2 = 4(\Delta x)^2$ $u[1,2]-2(1+\lambda)u[2,2] + \lambda u[2,1]$ $=4(\Delta x)^2- (\Delta y-1)^2 -\lambda (\Delta x-1)^2$ $u[1,2]-2(1+0.25)u[2,2] + 0.25 u[2,1]$ $=4(1/3)^2- (2/3-1)^2 -0.25(1/3-1)^2$ $u[1,2]-2.5u[2,2] + 0.25 u[2,1] =-0.2222$

Obtiene el sistema de ecuaciones a resolver:

$-2.5u[1,1]+u[2,1] +0.25 u[1,2] = -0.1111$ $u[1,1]-2.5u[2,1] + 0.25 u[2,2] = 0.2222$ $-2.5u[1,2]+u[2,2] + 0.25 u[1,1] = -0.5555$ $u[1,2]-2.5u[2,2] + 0.25 u[2,1] =-0.2222$

Se escribe en forma matricial Ax=B

$\begin {bmatrix} -2.5 && 1&&0.25&&0 \\1&&-2.5&&0&&0.25\\0.25&&0&&-2.5&&1\\0&&0.25&&1&&-2.5\end{bmatrix} \begin {bmatrix} u[1,1] \\ u[2,1]\\ u[1,2]\\u[2,2] \end{bmatrix} = \begin {bmatrix}-0.1111\\0.2222\\-0.5555\\-0.2222 \end{bmatrix}$

que se resuelve con un algoritmo:

import numpy as np

A = np.array([[-2.5,1,0.25,0],
              [1,-2.5,0,0.25],
              [0.25,0,-2.5,1],
              [0,0.25,1,-2.5]])
B = np.array([-0.1111,0.2222,-0.5555,-0.2222])

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

print(x)

y los resultados son:

[ 0.05762549 -0.04492835  0.31156835  0.20901451]

Ejercicio: 3Eva_IT2009_T3 Integrar Simpson compuesta

La fórmula a integrar es:

$\int_0^1 \frac{\cos (2x)}{x^{1/3}} \delta x$

Que tiene la forma:

da como resultado:

el área bajo la curva es:  0.879822622256

Algoritmo en Python

# 3Eva_IT2009_T3 Integrar Simpson compuesta
import numpy as np
import matplotlib.pyplot as plt

funcionx = lambda x: np.cos(2*x)/(x**(1/3))

# INGRESO
a = 0.00001
b = 1
tramos = 10000

# PROCEDIMIENTO
h = (b-a)/tramos
x = a
area = 0
for i in range(0,tramos,2):
    deltaA = (h/3)*(funcionx(x)+4*funcionx(x+h)+funcionx(x+2*h))
    area = area + deltaA
    x = x + 2*h
    
# para la gráfica
xi = np.linspace(a,b,tramos+1)
fxi = funcionx(xi)

print('el área bajo la curva es: ', area)

# Gráfica
plt.plot(xi,fxi)
plt.title('Funcion a integrar')
plt.grid()
plt.xlabel('x')
plt.show()

Tarea: convertir a función el cálculo de Simpson