Etiqueta: análisis numérico

Ejercicio: 2Eva_IT2018_T1 Paracaidista wingsuit

Plantear con: [ Runge-Kutta para f''(x) ] [ Runge-Kutta para f’(x) ]

..

---

a. Planteamiento con Runge-Kutta 2do Orden para Segunda derivada

La expresión:

\frac{dv}{dt} = g - \frac{cd}{m} v^2

se puede plantear sustituir la variable con v = -\frac{dy}{dt} al considerar el sentido contrario entre la velocidad al caer y la referencia de altura hacia arriba. Ver figura.

\frac{dy^2}{dt^2} = g - \frac{cd}{m} \Big( \frac{dy}{dt} \Big) ^2

Que es una EDO de 2do orden o como 2da derivada.

La solución se propone resolver de forma simultanea para t,y,v con Runge Kutta para segunda derivada donde:

f(t,y,v) = -v g(t,y,v) = g - \frac{cd}{m} v^2

Al sustituir los valores de las constantes en la ecuación como gravedad, masa e índice de arrastre se tiene:

f(t,y,v) = -v g(t,y,v) = 9.8 - \frac{0.225}{90} v^2

con las condiciones iniciales del ejercicio  t₀ = 0 , y₀ = 1000, v₀ = 0
la velocidad se inicia con cero, si el paracaidista se deja caer desde el borde el risco, como en el video adjunto al enunciado.

Para las iteraciones, recuerde que
t se corrige con t+h (en el algoritmo era la posición para x)
y se corrige con y+K1y
v se corrige con v+K1v (en el algoritmo era la posición para z)

itera = 0

K1y = h f(t,y,v) = 2(-(0)) = 0 K1v = h g(t,y,v) = 2(9.8 - \frac{0.225}{90} (0)^2) = 19.6

..
K2y = h f(t+h, y+K1y, v + K1v) = 2(-(0 + 19.6)) = -39.2

K2v = h g(t+h, y+K1y, v + K1v) = 2(9.8 - \frac{0.225}{90} (0+19.6)^2) =17.6792

..
y_1 = y_0 + \frac{K1y+K2y}{2} = 1000 + \frac{0-39.2}{2}= 980.4

v_1 = v_0 + \frac{K1v+K2v}{2} = 0 + \frac{19.6-17.67}{2} = 18.63 t_1 =t_0 + h = 0+2 = 2

 

ti	yi	vi	K1y	K1v	K2y	K2v
0	1000	0	0	0	0	0
2	980.4	18.63	0	19.6	-39.2	17.6792

itera = 1

K1y = 2(-(18.63)) = -37.2792 K1v = 2(9.8 - \frac{0.225}{90} (18.63)^2) = 17.8628

..
K2y =2(-(18.6396+17.8628)) =-73.00

K2v = 2(9.8 - \frac{0.225}{90} (18.6396+17.8628)^2) =12.9378

..
y_2 =980.4 + \frac{ -37.2792+(-73.00)}{2}= 925.25

v_2 = 18.63 + \frac{17.8628+12.9378}{2} = 34.0399 t_2 =t_1 + h = 2+2 = 4

ti	yi	vi	K1y	K1v	K2y	K2v
0	1000	0	0	0	0	0
2	980.4	18.63	0	19.6	-39.2	17.6792
4	925.25	34.0399	-37.2792	17.8628	-73.00	12.9378

itera = 2

K1y = h f(t,y,v) = 2(-(34.0399)) = -68.0798 K1v = h g(t,y,v) = 2(9.8 - \frac{0.225}{90} (34.0399)^2) = 13.8064

..
K2y = h f(t+h, y+K1y, v + K1v) = 2(-(34.0399+13.8064)) =-95.6927

K2v = h g(t+h, y+K1y, v + K1v) = 2(9.8 - \frac{0.225}{90} (34.0399+13.8064)^2) =8.1536

..
y_2 = 925.25 + \frac{ -68.0798+(-95.6927)}{2}= 843.3716

v_2 = 34.0399 + \frac{13.8064+8.1536}{2} = 45.0199 t_2 =t_1 + h = 4+2 = 6

ti	yi	vi	K1y	K1v	K2y	K2v
0	1000	0	0	0	0	0
2	980.4	18.63	0	19.6	-39.2	17.6792
4	925.25	34.0399	-37.2792	17.8628	-73.00	12.9378
6	843.3716	45.0199	-68.0798	13.8064	-95.6927	8.1536

Algoritmo con Python

Resultados
La velocidad máxima, si no hubiese límite en la altura, se encuentra en alrededor de 62.39 m/s.
Sin embargo, luego de 20 segundos se observa que la altura alcanza el valor de cero, es decir se alcanzó tierra con velocidad de  62 m/s, que son algo mas de 223 Km/h, el objeto se estrella…!

Resultados con el algoritmo:

Runge-Kutta Segunda derivada
i  [ xi,  yi,  zi ]
   [ K1y,  K1z,  K2y,  K2z ]
0 [   0. 1000.    0.]
   [0. 0. 0. 0.]
1 [  2.     980.4     18.6396]
  [  0.      19.6    -39.2     17.6792]
2 [  4.       925.257973  34.039945]
  [-37.2792    17.862827 -73.004853  12.937864]
3 [  6.       843.371672  45.019966]
  [-68.079891  13.806411 -95.692712   8.153631]
4 [  8.       743.865726  52.131168]
  [ -90.039933    9.466013 -108.971959    4.75639 ]
5 [ 10.       633.591684  56.485537]
  [-104.262336    6.011707 -116.285749    2.697031]
6 [ 12.       516.97369   59.069216]
  [-112.971073    3.646921 -120.264915    1.520438]
7 [ 14.       396.681119  60.575537]
  [-118.138432    2.154139 -122.446709    0.858504]
8 [ 16.       274.277023  61.445121]
  [-121.151075    1.253021 -123.657117    0.486147]
9 [ 18.       150.664295  61.944336]
  [-122.890243    0.722485 -124.335213    0.275943]
10 [20.       26.361127 62.230024]
   [-123.888671    0.414496 -124.717664    0.15688 ]
11 [ 22.       -98.336041  62.393224]
   [-1.244600e+02  2.371205e-01 -1.249343e+02  8.927924e-02]
12 [  24.       -223.257917   62.486357]
   [-1.247864e+02  1.354280e-01 -1.250573e+02  5.083841e-02]
13 [  26.       -348.307907   62.539475]
   [-1.249727e+02  7.727584e-02 -1.251273e+02  2.895913e-02]
14 [  28.       -473.430927   62.56976 ]
   [-1.250789e+02  4.407055e-02 -1.251671e+02  1.649935e-02]
15 [  30.       -598.595572   62.587023]
   [-1.251395e+02  2.512592e-02 -1.251898e+02  9.401535e-03]
16 [  32.       -723.783942   62.596863]
   [-1.251740e+02  1.432256e-02 -1.252027e+02  5.357469e-03]
>>> 

Instrucciones en Python

# 2Eva_IT2018_T1 Paracaidista wingsuit
import numpy as np
import matplotlib.pyplot as plt

def rungekutta2_fg(f,g,x0,y0,z0,h,muestras,
                   vertabla=False, precision = 6):
    ''' solucion a EDO con Runge-Kutta 2do Orden Segunda derivada,
        x0,y0 son valores iniciales, h es tamano de paso,
        muestras es la cantidad de puntos a calcular.
    '''
    tamano = muestras + 1
    tabla = np.zeros(shape=(tamano,3+4),dtype=float)

    # incluye el punto [x0,y0,z0]
    tabla[0] = [x0,y0,z0,0,0,0,0]
    xi = x0
    yi = y0
    zi = z0
    i=0
    if vertabla==True:
        print('Runge-Kutta Segunda derivada')
        print('i ','[ xi,  yi,  zi',']')
        print('   [ K1y,  K1z,  K2y,  K2z ]')
        np.set_printoptions(precision)
        print(i,tabla[i,0:3])
        print('  ',tabla[i,3:])
    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
        
        tabla[i] = [xi,yi,zi,K1y,K1z,K2y,K2z]
        if vertabla==True:
            txt = ' '
            if i>=10:
                txt='  '
            print(str(i)+'',tabla[i,0:3])
            print(txt,tabla[i,3:])
    return(tabla)


# PROGRAMA PRUEBA
# INGRESO
f = lambda t,y,v: -v # el signo, revisar diagrama cuerpo libre
g = lambda t,y,v: 9.8 - (0.225/90)*(v**2)
t0 = 0
y0 = 1000
v0 = 0
h  = 2
muestras = 15+1

# PROCEDIMIENTO
tabla = rungekutta2_fg(f,g,t0,y0,v0,h,muestras, vertabla=True)
ti = tabla[:,0]
yi = tabla[:,1]
vi = tabla[:,2]

# SALIDA
# print('tabla de resultados')
# print(tabla)

# GRAFICA
plt.subplot(211)
plt.plot(ti,vi,label='velocidad v(t)', color='green')
plt.plot(ti,vi,'o')
plt.ylabel('velocidad (m/s)')
plt.title('paracaidista Wingsuit con Runge-Kutta')
plt.legend()
plt.grid()

plt.subplot(212)
plt.plot(ti,yi,label='Altura y(t)',)
plt.plot(ti,yi,'o',)
plt.axhline(0, color='red')
plt.xlabel('tiempo (s)')
plt.ylabel('Altura (m)')
plt.legend()
plt.grid()

plt.show()

Plantear con: [ Runge-Kutta para f''(x) ] [ Runge-Kutta para f’(x) ]

..

---

b. Usando Runge-Kutta 2do Orden para Primera Derivada o velocidad,

El problema para un tiempo de observación t>0, se puede dividir en dos partes: velocidad y altura.

1. Determinar velocidad: Se aplica Runge-Kutta a la expresión con primera derivada o velocidad. Use  los valores iniciales dados, descarte calcular las alturas.
2. Determinar las altura:  Con los valores de velocidades y la altura inicial de 1 km = 1000 m puede integrar con trapecio para obtener la tabla de alturas. Se integra tramo a tramo.

Observe las unidades de medida y que la  velocidad es contraria  al eje de altura dy/dt = -v.

 

La expresión:

\frac{dv}{dt} = g - \frac{cd}{m} v^2 f(t,v) = g - \frac{cd}{m} v^2

itera = 0

K1v = h f(t,v) = 2(9.8)- \frac{0.225}{90}(0)^2) = 19.6 K2v = h f(t+h , v + K1v) = 2( 9.8 - \frac{0.225}{90}(0+19.6)^2) = 17.6792 v_1 = v_0 + \frac{K1y+K2y}{2} = 0 + \frac{19.6+17.6792}{2}= 18.6396 t_1 =t_0 + h = 0+2 = 2

itera = 1

K1v = 2(9.8 - \frac{0.225}{90} (18.6396)^2) = 17.8628 K2v = 2( 9.8 - \frac{0.225}{90} (18.6396+17.8628)^2) = 12.9379 v_1 = 18.6396 + \frac{17.8628+12.9379}{2}= 34.0399 t_1 =2+2 = 4

itera = 2

K1v = 2(9.8 - \frac{0.225}{90} (34.0399)^2) = 13.8064 K2v = 2( 9.8 - \frac{0.225}{90} (34.0399+13.8064)^2) = 8.1536 v_1 = 34.0399 + \frac{13.8064+8.1536}{2}= 45.02 t_1 = 4+2 = 6

las siguientes iteraciones se completan con el algoritmo.

Resultados

La velocidad máxima, si no hubiese límite en la altura, se encuentra en alrededor de 62.39 m/s.
Sin embargo, luego de 20 segundos se observa que la altura alcanza el valor de cero, es decir se alcanzó tierra con velocidad de  62 m/s, que son algo mas de 223 Km/h, el objeto se estrella…!

velocidad con Runge-Kutta primera derivada
 [tiempo, velocidad, K1,K2]
[[ 0.      0.      0.      0.    ]
 [ 2.     18.6396 19.6    17.6792]
 [ 4.     34.0399 17.8628 12.9379]
 [ 6.     45.02   13.8064  8.1536]
 [ 8.     52.1312  9.466   4.7564]
 [10.     56.4855  6.0117  2.697 ]
 [12.     59.0692  3.6469  1.5204]
 [14.     60.5755  2.1541  0.8585]
 [16.     61.4451  1.253   0.4861]
 [18.     61.9443  0.7225  0.2759]
 [20.     62.23    0.4145  0.1569]
 [22.     62.3932  0.2371  0.0893]]
velocidad con RK2 y altura con trapecio
 [tiempo, velocidad, altura]
[[   0.      0.   1000.  ]
 [   2.     18.64  981.36]
 [   4.     34.04  928.68]
 [   6.     45.02  849.62]
 [   8.     52.13  752.47]
 [  10.     56.49  643.85]
 [  12.     59.07  528.3 ]
 [  14.     60.58  408.65]
 [  16.     61.45  286.63]
 [  18.     61.94  163.24]
 [  20.     62.23   39.07]
 [  22.     62.39  -85.55]]
>>> 

Los cálculos se realizaron usando las instrucciones en Python:

# 2da Evaluación I Término 2018
# Tema 1. Paracaidista wingsuit
import numpy as np

def rungekutta2(d1y,x0,y0,h,muestras):
    # Runge Kutta de 2do orden
    tamano = muestras + 1
    tabla = np.zeros(shape=(tamano,2+2),dtype=float)
    
    # incluye el punto [x0,y0]
    tabla[0] = [x0,y0,0,0]
    xi = x0
    yi = y0
    for i in range(1,tamano,1):
        K1 = h * d1y(xi,yi)
        K2 = h * d1y(xi+h, yi + K1)

        yi = yi + (1/2)*(K1+K2)
        xi = xi + h
        
        tabla[i] = [xi,yi,K1,K2]
    return(tabla)

def integratrapecio_fi_tabla(xi,fi,y0):
    tamano = len(xi)
    yi = np.zeros(tamano,dtype = float)
    yi[0] = y0
    for i in range(1,tamano,1):
        h = xi[i]-xi[i-1]
        trapecio = h*(fi[i]+fi[i-1])/2
        yi[i]= yi[i-1] + trapecio
    return(yi)

# PROGRAMA -------------------------

# INGRESO
g = 9.8
cd = 0.225
m = 90
d1v = lambda t,v: g - (cd/m)*(v**2)

t0 = 0
v0 = 0
h = 2
y0 = 1000
muestras = 11

# PROCEDIMIENTO
velocidad = rungekutta2(d1v,t0,v0,h,muestras)
ti = velocidad[:,0]
vi = velocidad[:,1]

# Altura, velocidad es contraria altura,
# integrar en cada tramo por trapecios o Cuadratura de Gauss
altura = integratrapecio_fi_tabla(ti,-vi,y0)

# Tabla de resultados de tiempo, velocidad, altura
altura = np.transpose([altura])
tabla = np.concatenate((velocidad[:,0:2],altura), axis = 1)

# SALIDA
np.set_printoptions(precision=4)
print('velocidad con Runge-Kutta primera derivada')
print(' [tiempo, velocidad, K1,K2]')
print(velocidad)
np.set_printoptions(precision=2)
print('velocidad con RK2 y altura con trapecio')
print(' [tiempo, velocidad, altura]')
print(tabla)

Plantear con: [ Runge-Kutta para f''(x) ] [ Runge-Kutta para f’(x) ]

Ejercicio: 1Eva_IT2018_T1 Tanque esférico canchas deportivas

a) Para seleccionar el rango para h=[a,b], se observa que el tanque puede estar vacio, a=0 o lleno al máximo, b=2R = 2(3)=6, con lo que se obtiene:

h =[0.0, 6.0]

conociendo la proporción con el valor máximo, se tiene un valor inicial para h₀ para las iteraciones.

Vmax = \frac{\pi}{3} (2R)^2 (3R-2R) Vmax = \frac{4\pi }{3}R^{3}= 113.09 h_0 = (6)*30/113.09 = 1.59

b) Usar Newton-Raphson con tolerancia 1e-6

f(h) = \frac{\pi }{3}h^2 (3(3)-h)-30 f(h) = \frac{\pi }{3}(9h^2 -h^3-28.647) f'(h) = \frac{\pi }{3}(18h-3h^2)

el punto siguiente de iteración es:

h_{i+1} = h_{i} -\frac{f(h)}{f'(h)} = h_{i}-\frac{ \frac{\pi }{3}(9h^2 -h^3-28.647)}{ \frac{\pi }{3}(18h-3h^2)} h_{i+1} = h_{i} -\frac{(9h^2 -h^3-28.647)}{(18h-3h^2)}

con lo que se realiza la tabla de iteraciones:

hi	hi+1	error	orden
1.590	2.061	0.47	10-1
2.061	2.027	-0.034	10-2
2.027	2.02686	-0.00014	10-4
2.02686	2.0268689	-2.32E-09	10-9

En valor práctico 2.028 m usando flexómetro, a menos que use medidor laser con precisión 10^-6 usará más dígitos con un tanque de 6 metros de altura ganará una precisión de una gota de agua para usar en duchas o regar el césped .

c) El orden de convergencia del error observando las tres primeras iteraciones es cuadrático

Tarea: Realizar los cálculos con Python, luego aplique otro método. Añada sus observaciones y conclusiones.

Ejercicio: 1Eva_IT2018_T3 Temperatura en nodos de placa

a) Plantear el sistema de ecuaciones. Usando el promedio para cada nodo interior:

a=\frac{50+c+100+b}{4} b=\frac{a+30+50+d}{4} c=\frac{a+60+100+d}{4} d=\frac{b+60+c+30}{4}

que reordenando se convierte en:

4a=150+c+b 4b=a+80+d 4c=a+160+d 4d=b+c+90

simplificando:

4a-b-c= 150 a-4b+d = -80 a-4c+d = -160 b+c-4d = -90

que a forma matricial se convierte en:

A = [[ 4, -1, -1, 0.0],
     [ 1, -4,  0, 1.0],
     [ 1,  0, -4, 1.0],
     [ 0,  1,  1,-4.0]]
B = [[ 150.0],
     [ -80.0],
     [-160.0],
     [ -90.0]]

Observación: la matriz A ya es diagonal dominante, no requiere pivotear por filas.  Se aumentó el punto decimal a los valores de la matriz A y el vector B  para que sean considerados como números reales.

El número de condición es: np.linalg.cond(A) = 3.0

que es cercano a 1 en un orden de magnitud, por lo que la solución matricial es «estable» y los cambios en los coeficientes afectan proporcionalmente a los resultados. Se puede aplicar métodos iterativos sin mayores inconvenientes.

b y c) método de Jacobi para sistemas de ecuaciones, con vector inicial

 
X(0) = [[60.0],
        [40], 
        [70],
        [50]] 

reemplazando los valores iniciales en cada ecuación sin cambios.

iteración 1
a=\frac{50+70+100+40}{4} = 65

b=\frac{60+30+50+50}{4} = 47.5 c=\frac{60+60+100+50}{4} = 67.5 d=\frac{40+60+70+30}{4} = 50

X(1) = [[65],
        [47.5],
        [67.5],
        [50]]

vector de error = 
     [|65-60|,
      |47.5-40|,
      |67.5-70|,
      |50-50|]
  = [|5|,
     |7.5|,
     |-2.5|,
     |0|]
errormax = 7.5

iteración 2
a=\frac{50+67.5+100+47.5}{4} = 66.25

b=\frac{65+30+50+50}{4} = 48.75 c=\frac{65+60+100+50}{4} = 68.75 d=\frac{47.5+60+67.7+30}{4} = 51.3

X(2) = [[66.25],
        [48.75],
        [68.75],
        [51.3]]

vector de error = 
       [|66.25-65|,
        |48.75-47.5|,
        |68.75-67.5|,
           |51.3-50|] 
  = [|1.25|,
     |1.25|,
     |1.25|,
     |1.3|]
errormax = 1.3

iteración 3
a=\frac{50+68.75+100+48.75}{4} = 66.875

b=\frac{66.25+30+50+51.3}{4} = 49.38 c=\frac{66.25+60+100+51.3}{4} = 69.3875 d=\frac{48.75+60+68.75+30}{4} = 51.875

X(2) = [[66.875],
        [49.38],
        [69.3875],
        [51.875]]

vector de error = 
      [|66.875-66.25|,
       |49.38-48.75|,
       |69.3875-68.75|,
       |51.875-51.3|]
    = [|0.655|,
       |0,63|,
       |0.6375|,
        |0.575|]
errormax = 0.655
con error relativo de:
100*0.655/66.875 = 0.97%

siguiendo las iteraciones se debería llegar a:

>>> np.linalg.solve(A,B)
array([[ 67.5],
       [ 50. ],
       [ 70. ],
       [ 52.5]])

Ejercicio: 1Eva_IT2018_T4 El gol imposible

Tabla de datos:

ti = [0.00, 0.15, 0.30, 0.45, 0.60, 0.75, 0.90, 1.05, 1.20]
xi = [0.00, 0.50, 1.00, 1.50, 1.80, 2.00, 1.90, 1.10, 0.30]
yi = [0.00, 4.44, 8.88,13.31,17.75,22.19,26.63,31.06,35.50]
zi = [0.00, 0.81, 1.40, 1.77, 1.91, 1.84, 1.55, 1.03, 0.30]

Observe que, un gol simplificado con física básica es un tiro parabólico cuya trayectoria se compone de movimientos en los ejes, Y y Z. Sin embargo, lo «imposible» del gol mostrado es añadir el movimiento en X. (Referencia de la imagen en el enunciado)

El movimiento de «profundidad» o dirección hacia el arco y(t) es semejante a un polinomio de primer grado, y el movimiento de «altura» z(t) es un polinomio de segundo grado. El movimiento «imposible» en el eje X, podría ser descrito por un polinomio de segundo o mayor grado.

a) Encontrar t para altura máxima, que se encuentra al igualar la derivada dz/dt a cero. Para interpolar el polinomio z(t), de segundo grado, se puede usar tres puntos de los sugeridos: 0, 0.3 y 0.6, que además son equidistantes en t (facilita usar cualquier método de interpolación).

Por ejemplo, con diferencias finitas avanzadas:

t	z(t)	d1	d2	d3
0.00	0.00	1.40	-0.89
0.30	1.40	0.51
0.60	1.91

z(t) = 0 + 1.40\frac{(t-0)}{1!(0.3)} - 0.89 \frac{(t-0)(t-0.3)}{2!(0.3)^2} = 0 + 4.66 t - 4.94(t^2-0.3t) = 4.66 t + 1.48 t - 4.94 t^2 z(t) = 6.42 t - 4.94 t^2

para encontrar el valor máximo se encuentra \frac{dz}{dt} = 0

\frac{dz}{dt} = 6.42 - 2(4.94) t 6.42 - 2(4.94) t = 0 t = \frac{6.42}{2(4.94)} t = 0.649

Observe que el resultado tiene sentido, pues según la tabla, el máximo debería estar entre 0.60 y 0.75

b) El cruce por la «barrera», corresponde al desplazamiento del balón en el eje Y a 9 metros: y(t) = 9.
El polinomio modelo de primer grado usa dos puntos para la interpolación, de los sugeridos se escogen dos, por ejemplo: 0.0 y 0.3.

Usando diferencias finitas avanzadas :

d1 = (8.88-0) = 8.88 y(t) = 0 + 8.88\frac{(t-0)}{1!(0.3)} y(t) = 29.6 t

usando y(t) = 9

29.6 t = 9
t = 0.30
z(0.30) = 1.40
(de la tabla de datos)

cuya respuesta con solo dos dígitos decimales es coherente, al estar cercano el valor a una distancia y=8.88 o aproximado a 9.
La respuesta puede ser mejorada usando más digitos significativos en las operaciones.

c)  La desviación máxima en eje X.
Determine un polinomio para la trayectoria en el eje X y obtenga el máximo igualando la derivada del polinomio x(t) a cero.

Por simplicidad, se usa un polinomio de segundo grado, alrededor de los valores máximos en el eje X

t	x(t)	d1	d2	d3
0.60	1.80	0.20	-0.30
0.75	2.00	-0.10
0.90	1.90

x(t) = 1.80 + 0.20 \frac{(t-0.60)}{1!(0.15)} -0.30 \frac{(t-0.60)(t-0.75)}{2!(0.15)^2} x(t) = 1.80 + 1.33 (t-0.60) - 6.67(t-0.60)(t-0.75)

como se busca el máximo, se usa la derivada \frac{dx}{dt} = 0

\frac{dx}{dt} = 1.33 - 6.67(2t +(-0.60-0.75)) 1.33 - 13.34t + 9.00 = 0 -13.34t + 10.33 = 0

t = 0.77

x(0.77) = 1.80 + 1.33(0.77-0.60) - 6.67(0.77-0.60)(0.77-0.75) x(0.77) = 2.003

lo que es coherente con la tabla para el eje x, pues el máximo registrado es 2, y el próximo valor es menor, la curva será decreciente.

---

Desarrollo extra para observar y verificar resultados:

Usando los puntos y un graficador 3D se puede observar la trayectoria:

Tarea: Realice el ejercicio usando los algoritmos en Python, muestre los polinomios obtenidos y compare.

Nota: La gráfica 3D presentada usa interpolación de Lagrange con todos los puntos. Realice las observaciones y recomendaciones necesarias y presente su propuesta como tarea.

Ejercicio: 1Eva_IIT2017_T1 Aproximar a polinomio usando puntos

Se dispone de tres puntos para la gráfica.

x	f(x)
0	1
0.2	1.6
0.4	2.0

Si el polinomio de Taylor fuera de grado 0, sería una constante, que si se evalúa en x₀ = 0 para eliminar los otros términos, se encuentra que sería igual a 1

Como se pide el polinomio de grado 2, se tiene la forma:

p(x) = a + bx + c x ^2 p(x) = 1 + bx + c x^2

Se disponen de dos puntos adicionales que se pueden usar para determinar b y c:

p(0.2) = 1 + 0.2 b + (0.2)^2 c = 1.6 p(0.4) = 1 + 0.4 b + (0.4)^2 c = 2.0

simplificando:

0.2 b + (0.2)^2 c = 1.6 - 1 = 0.6 0.4 b + (0.4)^2 c = 2.0 - 1 = 1

multiplicando la primera ecuación por 2 y restando la segunda ecuación:

0 - 0.08 c = 1.2-1 = 0.2 c = - 0.2/0.08 = -2.5

sustituyendo el valor de c obtenido en la primera ecuación

0.2 b + 0.04(-2.5) = 0.6 0.2 b = 0.6 - 0.04(-2.5) = 0.6 + 0.1 = 0.7 b = 0.7/0.2 = 3.5

con lo que el polinomio queda:
p(x) = 1 + 3.5 x - 2.5 x^2

validando con python:
tomando los puntos de prueba:

xi = [ 0, 0.2, 0.4]
fi = [ 1, 1.6, 2 ]
>>>

se obtiene la gráfica:

se adjunta las instrucciones usadas para validar que el polinomio pase por los puntos requeridos.

# 1Eva_IIT2017_T1 Aproximar a polinomio usando puntos
import numpy as np
import matplotlib.pyplot as plt

# INGRESO
xi = [ 0, 0.2, 0.4]
fi = [ 1, 1.6, 2 ]

px = lambda x: 1 + 3.5*x - 2.5*(x**2)
a = -0.5
b = 1
muestras = 21

# PROCEDIMIENTO
xj = np.linspace(a,b,muestras)
pxj = px(xj)

# SALIDA
print(xj)
print(pxj)

# Gráfica
plt.plot(xj,pxj,label='p(x)')
plt.plot(xi,fi,'o', label='datos')

plt.xlabel('x')
plt.ylabel('f(x)')
plt.grid()
plt.legend()
plt.show()

Nota: Se puede intentar realizar los polinomios aumentando el grado, sin embargo cada término agrega un componente adicional a los términos anteriores por la forma (x – x₀)^k

literal b

se requiere el integral aproximado de f(x) en el intervalo limitado por los 3 puntos de la tabla:

\int_{0}^{0.4}f(x) dx

Esta aproximación con un polinomio es el concepto de integración numérica con la regla de Simpson de 1/3, tema desarrollado en la unidad 5

I_{S13} = \frac{0.2}{3} \Big(1+4(1.6)+ 2 \Big) = 0.62666

Ejercicio: 2Eva_IIT2017_T4 Problema con valor de frontera

\frac{d^2T}{dx^2} + \frac{1}{x}\frac{dT}{dx} +S =0 0 \leq x \leq 1

Las diferencias finitas a usar son:

\frac{dT}{dx} =\frac{T_{i+1} - T_{i-1}}{2h}+O(h^2) \frac{d^2T}{dx^2}=\frac{T_{i+1} - 2T_i + T_{i-1}}{h^2}+O(h^2)

que al reemplazar el la ecuación:

\frac{T_{i+1} - 2T_i + T_{i-1}}{h^2} + \frac{1}{x_i}\frac{T_{i+1} -T_{i-1}}{2h}+S=0 2x_i (T_{i+1} - 4h x_i T_i + T_{i-1}) + h (T_{i+1} - T_{i-1}) = -2h^2 S x_i T_{i+1}(2x_i + h) - 4x_i T_i + T_{i-1}(2x_i - h) = -2h^2 S x_i T_{i-1}(2x_i - h) - 4x_i T_i + T_{i+1}(2x_i + h)= -2h^2 S x_i

con lo que se puede crear un sistema de ecuaciones para cada valor xi con T₀=2 y T₄ =1 que son parte del enunciado del problema.

Siendo h = 0.25 = 1/4,  y se indica al final que S=1, se crea un sistema de ecuaciones a resolver,

x = [0, 1/4, 1/4, 3/4, 1]

T_{i-1}\Big[2x_i - \frac{1}{4} \Big] - 4x_i T_i + T_{i+1}\Big[2x_i + \frac{1}{4} \Big] = -2\Big(\frac{1}{4}\Big)^2 (1) x_i T_{i-1}\Big[2x_i -\frac{1}{4}\Big] - 4x_i T_i + T_{i+1}\Big[2x_i + \frac{1}{4}\Big] = -\frac{1}{8} x_i

se sustituye con los valores conocidos para cada i:

i=1: 
T0[2(1/4) - 1/4] - 4(1/4)T1 + T2[2(1/4) + 1/4] = -(1/8)(1/4)

     - T1 + (3/4)T2 = -1/32 - (1/4)(2)
     - T1 + (3/4)T2 = -17/32

i=2: 
T1[2(1/2) - 1/4] - 4(1/2)T2 + T3[2(1/2) + 1/4] = -(1/8)(1/2)

     (3/4)T1 - 2T2 + (5/4)T3 = -1/16

i=3: 
T2[2(3/4) - 1/4] - 4(3/4)T3 + T4[2(3/4) + 1/4] = -(1/8)(3/4)

     (5/4)T2 - 3T3 = -3/32 - (7/4)(1)
     (5/4)T2 - 3T3 = -59/32

se ponen las ecuaciones en matrices para resolver con algun metodo numérico:

A = [[ -1, 3/4,   0],
     [3/4,  -2, 5/4],
     [  0, 5/4,  -3]]
B = [-17/32, -1/16, -59/32]
np.linalg.solve(A,B)
array([ 1.54,  1.34,  1.17])

Ejercicio: 2Eva_IIT2017_T1 EDO Runge Kutta 2do Orden d2y/dx2

Tema 1

Runge kutta de 2do Orden
f: y' = z
g: z' = .....
K1y = h f(xi, yi, zi)
K1z = h g(xi, y1, zi)

K2y = h f(xi+h, yi+K1y, zi+K1z)
K2z = h g(xi+h, yi+K1y, zi+K1z)

y(i+1) = yi + (1/2)(K1y + K2y)
z(i+1) = zi + (1/2)(K1z + K2z)

x(i+1) = xi + h
E = O(h3) 
xi ≤ z ≤ x(i+1)

f: z = Θ’
g: z’ = (-gr/L) sin(Θ)

Θ(0) = π/6
z(0) = 0

h=0.1

i=0, t0 = 0, Θ0 = π/6, z0 = 0
    K1y = 0.1(0) = 0
    K1z = 0.1(-9.8/2)sin(π/6) = -0.245

    K2y = 0.1(0+(-0.245)) = -0.0245
    K2z = 0.1(-9.8/2)sin(π/6+0) = -0.245

    Θ1 = π/6 + (1/2)(0+(-0.0245)) = 0.51139
    z1 = 0 + (1/2)(-0.245-0.245) = -0.245
    t1 = 0 + 0.1 = 0.1

i=1, t1 = 0.1, Θ1 = 0.51139, z1 = -0.245
    K1y = 0.1(-0.245) = -0.0245
    K1z = 0.1(-9.8/2)sin(0.51139) = -0.23978

    K2y = 0.1(-0.245+(-0.0245)) = -0.049
    K2z = 0.1(-9.8/2)sin(0.51139+(-0.0245)) = -0.22924

    Θ2 = 0.51139 + (1/2)(-0.0245+(-0.049)) = 0.47509
    z2 = -0.245 + (1/2)(-0.23978+(-0.22924)) = -0.245
    t2 = 0.1 + 0.1 = 0.2

   t         theta     z
[[ 0.        0.523599  0.      ]
 [ 0.1       0.511349 -0.245   ]
 [ 0.2       0.47486  -0.479513]
 [ 0.3       0.415707 -0.692975]
 [ 0.4       0.336515 -0.875098]
 [ 0.5       0.240915 -1.016375]
 [ 0.6       0.133432 -1.108842]
 [ 0.7       0.019289 -1.14696 ]
 [ 0.8      -0.09588  -1.128346]
 [ 0.9      -0.206369 -1.054127]
 [ 1.       -0.306761 -0.92877 ]
 [ 1.1      -0.39224  -0.759461]
 [ 1.2      -0.458821 -0.555246]
 [ 1.3      -0.503495 -0.326207]
 [ 1.4      -0.524294 -0.082851]
 [ 1.5      -0.520315  0.164197]
 [ 1.6      -0.491715  0.404296]
 [ 1.7      -0.439718  0.62682 ]
 [ 1.8      -0.366606  0.821313]
 [ 1.9      -0.275693  0.977893]
 [ 2.       -0.171235  1.087942]]

Literal b), con h= 0.25, con t = 1 ángulo= -0.352484

   t         theta     z
[[ 0.        0.523599  0.      ]
 [ 0.25      0.447036 -0.6125  ]
 [ 0.5       0.227716 -1.054721]
 [ 0.75     -0.070533 -1.170971]
 [ 1.       -0.352484 -0.910162]
 [ 1.25     -0.527161 -0.363031]
 [ 1.5      -0.540884  0.299952]
 [ 1.75     -0.387053  0.890475]
 [ 2.       -0.106636  1.221932]]

El error de del orden h³

---

Instruccciones en python, usando el algoritmo desarrollado en clase

# Runge Kutta de 2do
# EDO de 2do orden con condiciones de inicio
import numpy as np
import matplotlib.pyplot as plt

def rungekutta2_fg(f,g,v0,h,m):
    tabla = [v0]
    xi = v0[0]
    yi = v0[1]
    zi = v0[2]
    for i in range(0,m,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)

        yi1 = yi + (1/2)*(K1y+K2y)
        zi1 = zi + (1/2)*(K1z+K2z)
        xi1 = xi + h
        vector = [xi1,yi1,zi1]
        tabla.append(vector)

        xi = xi1
        yi = yi1
        zi = zi1
    tabla = np.array(tabla)
    return(tabla)

# Programa Prueba
# Funciones
f = lambda x,y,z : z
g = lambda x,y,z : (-gr/L)*np.sin(y)

gr = 9.8
L = 2

x0 = 0
y0 = np.pi/6
z0 = 0

v0 = [x0,y0,z0]

h = 0.1
xn = 2
m = int((xn-x0)/h)

# PROCEDIMIENTO
tabla = rungekutta2_fg(f,g,v0,h,m)

xi = tabla[:,0]
yi = tabla[:,1]
zi = tabla[:,2]

# SALIDA
np.set_printoptions(precision=6)
print('x, y, z')
print(tabla)
plt.plot(xi,yi, label='y')
plt.plot(xi,zi, label='z')
plt.legend()
plt.grid()
plt.show()

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

Ejercicio: 2Eva_IIT2017_T3 EDP parabólica con diferencias regresivas

\frac{dU}{dt} - \frac{1}{16} \frac{d^2U}{dx^2} = 0

Las diferencias finitas requidas en el enunciado son:

U'(x_i,t_j) = \frac{U(x_{i},t_j)-U(x_{i},t_{j-1})}{\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)

La indicación de regresiva es para la primera derivada, dependiente de tiempo t.

que al reemplazar en la ecuación sin el término de error, se convierte a.

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

Se reordenan los términos de forma semejante al modelo planteado en el método básico:

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

Se simplifica haciendo que haciendo

\lambda = \frac{\Delta t}{16\Delta x^2}

Cambiando la nomenclatura con solo los índices para las variables x y t, ordenando en forma ascendente los índices previo a crear el algoritmo.

\lambda[U(i+1,j)-2U(i,j)+U(i-1,j)] = U(i,j)-U(i,j-1)

Se reordena la ecuación como modelo para el sistema de ecuaciones.

\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) + Q U(i,j) + R U(i+1,j) = -U(i,j-1)

Se calculan los valores constantes:

λ = dt/(16*dx2) = 0.05/[16*(1/3)2] = 0.028125

P = λ = 0.028125
Q = (-1-2λ) = (1-2*0.028125) = -1.05625
R = λ = 0.028125

Usando las condiciones del problema:

U(0,t) = U(1,t) = 0, entonces, Ta = 0, Tb = 0

Para los valores de la barra iniciales se debe usar un vector calculado como 2sin(π x) en cada valor de x_i espaciados por hx = 1/3, x entre [0,1]

xi  = [0,1/3, 2/3, 1]
U[xi,0] = [2sin (0*π), 2sin(π/3), 2sin(2π/3), 2sin(π)]
U[xi,0] = [0, 2sin(π/3), 2sin(2π/3), 0]
U[xi,0] = [0, 1.732050,  1.732050, 0]

Con lo que se puede plantear las ecuaciones:

j=1: i=1
0.028125 U(0,1) + (-1.05625) U(1,1) + 0.028125 U(2,1) = -U(1,0)

j=1: i=2
0.028125 U(1,1) + (-1.05625) U(2,1) + 0.028125 U(3,1) = -U(2,0)

y reemplazando los valores de la malla conocidos:

0.028125 (0) – 1.05625 U(1,1) + 0.028125 U(2,1) = -1.732050
0.028125 U(1,1) – 1.05625 U(2,1) + 0.028125 (0) = -1.732050

hay que resolver el sistema de ecuaciones:

-1.05625  U(1,1) + 0.028125 U(2,1) = -1.732050
 0.028125 U(1,1) - 1.05625  U(2,1) = -1.732050

A = [[-1.05625 ,  0.028125],
     [ 0.028125, -1.05625 ]]
B = [-1.732050,-1.732050]
que resuelta con un método numérico:
[ 1.68,  1.68]

Por lo que la solución para una gráfica, con los índices de (fila,columna) como (t,x):

U = [[0, 1.732050,  1.732050, 0],
     [0, 1.680000,  1,680000, 0]]

El error del procedimiento, tal como fué planteado es del orden de O(Δt) y O(Δx²), o error de truncamiento E = O(Δx²) + O(Δt). Δt debe ser menor que Δx en aproximadamente un orden de magnitud

---

Usando algoritmo en python.

Usando lo resuelto en clase y laboratorio, se comprueba la solución con el algoritmo, con hx y ht mas pequeños y más iteraciones:

# EDP parabólicas d2u/dx2  = K du/dt
# método implícito
# Referencia: Chapra 30.3 p.895 pdf.917
#       Rodriguez 10.2.5 p.417
import numpy as np
import matplotlib.pyplot as plt

# INGRESO
# Valores de frontera
Ta = 0
Tb = 0
# longitud en x
a = 0
b = 1
# Constante K
K = 16
# Tamaño de paso
dx = 0.1
dt = 0.01
# temperatura en barra
tempbarra = lambda x: 2*np.sin(np.pi*x)
# iteraciones
n = 100

# PROCEDIMIENTO
# Valores de x
xi = np.arange(a,b+dx,dx)
m = len(xi)

# Resultados en tabla de u
u = np.zeros(shape=(m,n), dtype=float)
# valores iniciales de u
j=0
u[0,j] = Ta
for i in range(1,m-1,1):
    u[i,j] = tempbarra(xi[i])
u[m-1,j] = Tb

# factores P,Q,R
lamb = dt/(K*dx**2)
P = lamb
Q = -1 -2*lamb
R = lamb
vector = np.array([P,Q,R])
tvector = len(vector)

# Calcula U para cada tiempo + dt
j=1
while not(j>=n):
    u[0,j] = Ta
    u[m-1,j] = Tb
    # Matriz de ecuaciones
    tamano = m-2
    A = np.zeros(shape=(tamano,tamano), dtype = float)
    B = np.zeros(tamano, dtype = float)
    for f in range(0,tamano,1):
        for c in range(0,tvector,1):
            c1 = f+c-1
            if(c1>=0 and c1<tamano):
                A[f,c1] = vector[c]
        B[f] = -u[f+1,j-1]
    B[0] = B[0]-P*u[0,j]
    B[tamano-1] = B[tamano-1]-R*u[m-1,j]
    # Resuelve sistema de ecuaciones
    C = np.linalg.solve(A, B) 
    # copia resultados a u[i,j]
    for f in range(0,tamano,1):
        u[f+1,j] = C[f]
    j=j+1 # siguiente iteración
        
# 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, '.m')
plt.xlabel('x[i]')
plt.ylabel('t[j]')
plt.title('Solución EDP parabólica')
plt.show()

Ejercicio: 2Eva_IIT2017_T2 Volumen de isla

isla = np.array([[0,1,0,0,0],
                 [1,3,1,1,0],
                 [5,4,3,2,0],
                 [0,0,1,1,0]])

xi = np.array([0,100,200,300,400])
yi = np.array([0, 50,100,150])

Tamaño de la matriz: n=4, m=5

cantidad de elementos por fila impar, aplica Simpson 1/3
hx = (200-0)/2 =100
fila=0
    vector = [0,1,0,0,0]
    deltaA = (100/3)(0+4(1)+0) = 4(100/3)
    deltaA = (100/3)(0+4(0)+0) = 0
    area0 = 4(100/3) + 0 = 4(100/3)
fila=1
    vector = [1,3,1,1,0]
    deltaA = (100/3)(1+4(3)+1) = 14(100/3)
    deltaA = (100/3)(1+4(1)+0) = 5(100/3)
    area1 = 14(100/3) + 5(100/3) = 19(100/3)
fila=2
    vector = [5,4,3,2,0]
    deltaA = (100/3)(5+4(4)+3) = 24(100/3)
    deltaA = (100/3)(3+4(2)+0) = 11(100/3)
    area2 = 24(100/3) + 11(100/3) = 35(100/3)
fila=3
    vector = [0,0,1,1,0]
    deltaA = (100/3)(0+4(0)+1) = (100/3)
    deltaA = (100/3)(1+4(1)+0) = 5(100/3)
    area3 = (100/3) + 5(100/3) = 6(100/3)

areas = [ 4(100/3), 19(100/3), 35(100/3), 6(100/3)]
areas = (100/3)[ 4, 19, 35, 6 ]

areas tiene cantidad de elementos par, aplica Simpson 3/8
    hy = (150-0)/3 = 50
    deltaV = (3/8)(50)(100/3)(4+3(19) + 3(35)+ 6)
           = (25*25)(168)
    Volumen = 107500

tramos:  4 5
areas:  [  133.33333333   633.33333333  1166.66666667    66.66666667]
Volumen:  107500.0

las instrucciones en python para encontrar el valor son:

# 2da Eval II T 2017. Tema 2
# Formula de simpson
# Integración: Regla Simpson 1/3 y 3/8
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import axes3d

def simpson13(xi,yi):
    '''
    las muestras deben ser impares
    '''
    area = 0
    muestras = len(xi)
    impar = muestras%2
    if impar == 1:
        for i in range(0,muestras-2,2):
            h = (xi[i+2] - xi[i])/2
            deltaA = (h/3)*(yi[i]+4*yi[i+1]+yi[i+2])
            area = area + deltaA
    return(area)

def simpson38(xi,yi):
    '''
    las muestras deben ser pares
    '''
    area = 0
    muestras = len(xi)
    impar = muestras%2
    if impar == 0:
        for i in range(0,muestras-3,3):
            h = (xi[i+3] - xi[i])/3
            deltaA = (3*h/8)*(yi[i]+3*yi[i+1]+3*yi[i+2]+yi[i+3])
            area = area + deltaA
    return(area)

def simpson(xi,yi):
    '''
    Selecciona el tipo de algoritmo Simpson
    '''
    muestras = len(xi)
    impar = muestras%2
    if impar == 1:
        area = simpson13(xi,yi)
    else:
        area = simpson38(xi,yi)
    return(area)
    

# INGRESO
isla = np.array([[0,1,0,0,0],
                 [1,3,1,1,0],
                 [5,4,3,2,0],
                 [0,0,1,1,0]])

xi = np.array([0,100,200,300,400])
yi = np.array([0, 50,100,150])

# PROCEDIMIENTO
tamano = np.shape(isla)
n = tamano[0]
m = tamano[1]

areas = np.zeros(n,dtype = float)
for fila in range(0,n,1):
    unafranja = isla[fila,:]
    areas[fila] = simpson(xi,unafranja)
volumen = simpson(yi,areas)

# SALIDA
print('tramos: ', n,m)
print('areas: ', areas)
print('Volumen: ', volumen)

# Gráfica
X, Y = np.meshgrid(xi, yi)
fig = plt.figure()
ax = fig.add_subplot(111, projection = '3d')
ax.plot_wireframe(X,Y,isla)
plt.show()