notebook.community

Edit and run

Taller 3 - Métodos Numéricos

Sergio Agudelo Bernal
Código: 285675

1. Interpolación por método de Newton:

En estudios de polimerización inducida por radiación, se emplea una fuente de rayos $\gamma$ $ $ para obtener dosis medidas en radiación. Sin embargo, la dosis varía con la posición del aparato, según los datos que se dan a continuación:

</tr> </table> </center>

(a) ¿Cuál es la estimación para el nivel de dosis en 2.5 pulgadas?

Posición
(en pulgadas)

</td>

Dosis
($\times 10^5$ rads/h)

1.0

2.71

1.5

2.98

2.0

3.20

3.0

3.20

 





In [1]:



    


%load_ext sympy.interactive.ipythonprinting
from IPython.display import display, Math, Latex
from sympy import *
from numpy import *
from __future__ import division
import sympy as sym
import numpy as np

Pos = array([1.00,1.50,2.00,3.00])
Dosis = array([2.71,2.98,3.20,3.20])

print (Pos)
print (Dosis)
display(Math(r'Dosis = {} \times 10^5 rads/h'.format(interp(2.5,Pos,Dosis))))



    


















    






[ 1.   1.5  2.   3. ]
[ 2.71  2.98  3.2   3.2 ]










    







$$Dosis = 3.2 \times 10^5 rads/h$$

(b) Si se efectúa una nueva medición que indica que a 3.5 pulgadas el nivel de dosis correspondiente es de 2.98, ¿cuál será ahora la estimación para el nivel de dosis en 2.5 pulgadas?

Se generan de nuevo los arreglos, añadiendo los nuevos datos:





In [2]:



    


Pos = np.append(Pos,3.5)
Dosis = np.append(Dosis,2.98)

print (Pos)
print (Dosis)
display(Math(r'Dosis \approx {} \times 10^5 rads/h'.format(interp(2.5,Pos,Dosis))))



    


















    






[ 1.   1.5  2.   3.   3.5]
[ 2.71  2.98  3.2   3.2   2.98]










    







$$Dosis \approx 3.2 \times 10^5 rads/h$$

2. Interpolación por método de Newton:

Para los valores siguientes

</tr> </table> </center>

Donde $E$ son los voltios y $P$ los kilowatts en una curva de pérdida en el núcleo para un motor eléctrico.

(a) Elaborar una tabla de diferencias divididas

$E$

</td>

$P$

$40$

$0.63$

$60$

$1.36$

$80$

$2.18$

$100$

$3.00$

$120$

$3.93$

 





In [3]:



    


E = array([40, 60, 80, 100, 120])
P = array([0.63, 1.36, 2.18, 3.00, 3.93])
Diff1 = np.zeros(4)
Diff2 = np.zeros(3)
Diff3 = np.zeros(2)

#Primera Iteración
for i in range(1, len(E)):
    Diff1[i-1]  = ((P[i] - P[i-1])/(E[i] - E[i-1]))
    
#Segunda Iteración
for i in range(2,len(E)):
    Diff2[i-2] = ((Diff1[i-1] - Diff1[i-2])/(E[i] - E[i-2]))

#Tercera Iteración
for i in range(3,len(E)):
    Diff3[i-3] = ((Diff2[i-2] - Diff2[i-3])/(E[i] - E[i-3]))
    
#Cuarta y Última Iteración
Diff4 = ((Diff3[1] - Diff3[0])/(E[4] - E[0]))

print (Diff1)
print (Diff2)
print (Diff3)
print (Diff4)



    


















    






[ 0.0365  0.041   0.041   0.0465]
[  1.12500000e-04  -1.73472348e-19   1.37500000e-04]
[ -1.87500000e-06   2.29166667e-06]
5.20833333333e-08

Se alcanza a observar que el segundo valor de la segunda columna de las diferencias divididas debería ser cero, pero por manejo numérico en Python, no alcanza tal valor. Pero de acá en adelante se toma su valor como nulo:

</tr> </table> </center>

(b) Calcular el polinomio de interpolación de Newton de segundo grado para $E = 80,\ 100,\ 120$. Utilizarlo para estimar el valor de $P$ correspondiente a $E = 90$ voltios.

Los tres coeficientes ubicados diagonalmente a $y = 2.18$ en la tabla de diferencias divididas son los necesarios para armar el polinomio de segundo grado, de la siguiente manera:

$$P_2(x) = f[x_2,x_3,x_4](x-x_2)(x-x_3)+f[x_2,x_3](x-x_2)+f[x_2]$$

$x_k$

</td>

$f[x_k]$

$f[x_{k-1},x_k]$

$f[x_{k-2},x_{k-1},x_k]$

$f[x_{k-3},...,x_k]$

$f[x_{k-4},...,x_k]$

$40$

$0.63$

$—$

$—$

$—$

$—$

$60$

$1.36$

$0.0365$

$—$

$—$

$—$

$80$

$2.18$

$0.041$

$1.125 \times 10^{-4}$

$—$

$—$

$100$

$3.00$

$0.041$

$0$

$-1.875 \times 10^{-6}$

$—$

$120$

$3.93$

$0.0465$

$1.375 \times 10^{-4}$

$2.292 \times 10^{-6}$

$5.208 \times 10^{-8}$

 





In [4]:



    


x = Symbol('x')
pol = 1.375e-4*(x-80)*(x-100) + 0.041*(x-80) + 2.18
pol



    


















    
Out[4]:







$$0.041 x + \left(0.0001375 x - 0.011\right) \left(x - 100\right) - 1.1$$
















In [5]:



    


expand(pol)



    


















    
Out[5]:







$$0.0001375 x^{2} + 0.01625 x$$
















In [6]:



    


display(Math(r'P_2(90) \approx {}\ kW'.format(pol.subs(x,90).evalf(5))))



    


















    







$$P_2(90) \approx 2.5762\ kW$$

3. Ajustar los siguientes datos a una línea recta:

</tr> </table> </center>


Grafique la nube de puntos y estime el valor de la función en 1.

El ajuste para una recta tiene los siguientes parámetros:

$$ \left[ \begin{array}{ccc} N & \sum x_k \\ \sum x_k & \sum x^2_k \end{array} \right] \left[ \begin{array}{ccc} C_1 \\ C_2 \end{array} \right] = \left[ \begin{array}{ccc} \sum y_k \\ \sum x_k y_k \end{array} \right] \\ $$
$$ y = C_2 x + C_1 $$

$x$

</td>

$f(x)$

$-1$

$2$

$0$

$0$

$3$

$4$

$7$

$7$

 





In [7]:



    


%pylab inline

x = array([-1, 0, 3, 7])
y = array([2, 0, 4, 7])
xsq = x**2
xy = x*y

plot(x,y,'k-')

A = array([[len(x), sum(x)],[sum(x), sum(xsq)]])
B = array([sum(y), sum(xy)])
print (A)
print (B)



    


















    






Welcome to pylab, a matplotlib-based Python environment [backend: module://IPython.zmq.pylab.backend_inline].
For more information, type 'help(pylab)'.
[[ 4  9]
 [ 9 59]]
[13 59]










    
























In [8]:



    


C = linalg.solve(A,B)
C



    


















    
Out[8]:








array([ 1.52258065,  0.76774194])

La ecuación es:

$$y \approx 0.768 x + 1.523$$

Ahora, se compara la aproximación encontrada con los datos, y se encuentra el valor estimado en $x = 1$:





In [9]:



    


def f(X):
    return X*C[1]+C[0]

p1, = plot(x, y, 'ko-')
p2, = plot(x, f(x), 'r')
xlim([-1,7])
ylim([-1,8])
legend([p1, p2], ["Datos", "Ajuste"], loc=2)
show()
display(Math(r'f(1) \approx {}'.format(f(1))))



    


















    

















    







$$f(1) \approx 2.2903225806451615$$

4. Ajustar los siguientes datos a una línea recta:

</tr> </table> </center>



Grafique la nube de puntos y estime el valor de la función en 10.

De manera idéntica al ejercicio anterior, se encuentra la matriz:

$x$

</td>

$f(x)$

$1$

$0.20$

$2$

$0.25$

$3$

$0.20$

$4$

$0.35$

$5$

$0.45$

$6$

$0.40$

 





In [10]:



    


%pylab inline

x2 = array([1, 2, 3, 4, 5, 6])
y2 = array([0.20, 0.25, 0.20, 0.35, 0.45, 0.40])
x2sq = x2**2
x2y2 = x2*y2

plot(x2,y2,'k-')

A2 = array([[len(x2), sum(x2)],[sum(x2), sum(x2sq)]])
B2 = array([sum(y2), sum(x2y2)])
print (A2)
print (B2)



    


















    






Welcome to pylab, a matplotlib-based Python environment [backend: module://IPython.zmq.pylab.backend_inline].
For more information, type 'help(pylab)'.
[[ 6 21]
 [21 91]]
[ 1.85  7.35]










    
























In [11]:



    


C2 = linalg.solve(A2,B2)
C2



    


















    
Out[11]:








array([ 0.13333333,  0.05      ])

La ecuación es:

$$y \approx 0.05 x + 0.1\overline{3}$$

Ahora, se compara la aproximación encontrada con los datos, y se encuentra el valor estimado en $x = 10$:





In [12]:



    


def f2(X):
    return X*C2[1]+C2[0]

p3, = plot(x2, y2, 'ko-')
p4, = plot(x2, f2(x2), 'r')
xlim([1,6])
ylim([0.1,0.5])
legend([p3, p4], ["Datos", "Ajuste"], loc=2)
show()
display(Math(r'f(10) \approx {}'.format(f2(10))))



    


















    

















    







$$f(10) \approx 0.6333333333333333$$

5. Determine a qué función se ajusta mejor a los siguientes puntos:

</tr> </table> </center>



Grafique la nube de puntos y estime el valor de la función en 0 y en 2.

Inicialmente se grafican los puntos para saber cuál función es la más adecuada:

$x$

</td>

$f(x)$

$-3$

$14$

$-1$

$4$

$1$

$2$

$3$

$8$

$5$

$22$

$7$

$44$

 





In [13]:



    


%pylab inline

x3 = array([-3, -1, 1, 3, 5, 7])
y3 = array([14, 4, 2, 8, 22, 44])

plot(x3,y3,'k-')



    


















    






Welcome to pylab, a matplotlib-based Python environment [backend: module://IPython.zmq.pylab.backend_inline].
For more information, type 'help(pylab)'.










    
Out[13]:








[<matplotlib.lines.Line2D at 0xb09d1c8c>]

Parece asemejarse a una parábola. Entonces, se toma el ajuste para una ecuación de orden 2:

$$ \left[ \begin{array}{ccc} N & \sum x_k & \sum x^2_k \\ \sum x_k & \sum x^2_k & \sum x^3_k \\ \sum x^2_k & \sum x^3_k & \sum x^4_k \end{array} \right] \left[ \begin{array}{ccc} C_1 \\ C_2 \\ C_3 \end{array} \right] = \left[ \begin{array}{ccc} \sum y_k \\ \sum x_k y_k \\ \sum x^2_k y_k \end{array} \right] \\ $$
$$ y = C_3 x^2 + C_2 x + C_1 $$ 





In [14]:



    


x3sq = x3**2
x3cub = x3**3
x3fth = x3**4
x3y3 = x3*y3
x3sqy3 = (x3**2)*y3

A3 = array([[len(x3), sum(x3), sum(x3sq)],[sum(x3), sum(x3sq), sum(x3cub)],[sum(x3sq), sum(x3cub), sum(x3fth)]])
B3 = array([sum(y3), sum(x3y3), sum(x3sqy3)])
print (A3)
print (B3)



    


















    






[[   6   12   94]
 [  12   94  468]
 [  94  468 3190]]
[  94  398 2910]

















In [15]:



    


C3 = linalg.solve(A3,B3)
C3



    


















    
Out[15]:








array([ 2., -1.,  1.])

La ecuación es:

$$y = x^2 - x + 2$$

Ahora, se compara la aproximación encontrada con los datos, y se encuentran los valores estimado en $x = 0,\ x = 2$:





In [16]:



    


def f3(X):
    return (X**2)*C3[2]+X*C3[1]+C[0]

xx = linspace(-3, 7, num=1000)
p5, = plot(x3, y3, 'ko-')
p6, = plot(xx, f3(xx), 'r')
xlim([-3,7])
ylim([1,45])
legend([p5, p6], ["Datos", "Ajuste"], loc=2)
show()
display(Math(r'f(0) = {}'.format(round(f3(0), 2)))) # Se sabe que Python no redondeó las cifras, 
display(Math(r'f(2) = {}'.format(round(f3(2), 2)))) # entonces se redondea forzadamente aquí
                                                    # únicamente para presentar bien los datos.



    


















    

















    







$$f(0) = 1.52$$









    







$$f(2) = 3.52$$

6. Proyectar la oferta de un cierto producto tomando en cuenta los datos obtenidos en el estudio de mercado, ver cuál de los métodos o curvas de proyección se ajusta mejor a la nube de puntos y determinar la oferta para los próximos diez años.

</tr> </table> </center>

Lo primero que se hace es graficar los puntos, y así estimar la ecuación necesaria para el ajuste:

Año

</td>

Tiempo$(x)$

Oferta$(y)$

$1989$

$1$

$100\,000$

$1990$

$2$

$120\,000$

$1991$

$3$

$140\,000$

$1992$

$4$

$110\,000$

$1993$

$5$

$170\,000$

$1994$

$6$

$150\,000$

$1995$

$7$

$180\,000$

$1996$

$8$

$200\,000$

$1997$

$9$

$210\,000$

$1998$

$10$

$200\,000$

 





In [17]:



    


t = array([1, 2, 3, 4, 5, 6, 7, 8, 9, 10])
O = array([100000, 120000, 140000, 110000, 170000, 150000, 180000, 200000, 210000, 200000])
p7, = plot(t, O, 'ko-')
xlim([1,10])
ylim([0,220000])
show()

La oferta parece tener tendencia lineal, pero también podría ser exponencial. El propósito es realizar ambas estimaciones y compararlas.

(a) Ajuste lineal:





In [18]:



    


tsq = t**2
tO = t*O

M1 = array([[len(t), sum(t)],[sum(t), sum(tsq)]])
M2 = array([sum(O), sum(tO)])
print (M1)
print (M2)



    


















    






[[ 10  55]
 [ 55 385]]
[1580000 9700000]

















In [19]:



    


Coeff = linalg.solve(M1,M2)
Coeff



    


















    
Out[19]:








array([ 90666.66666667,  12242.42424242])

La ecuación lineal resultante es:

$$y = 12242.\overline{42} x + 90666.\overline{6}$$

(b) Aproximación exponencial:





In [20]:



    


lnO = log(O)
tlnO = t*lnO
M3 = array([sum(lnO), sum(tlnO)])
print (M3)



    


















    






[ 119.39547035  663.32072503]

















In [21]:



    


CoeffExp = linalg.solve(M1,M3)
CoeffExp



    


















    
Out[21]:








array([ 11.4965045 ,   0.08055319])

La ecuación exponencial resultante es:

$$y \approx e^{11.4965}\, e^{0.08055 x}$$

Finalmente, se grafican tanto los datos como las aproximaciones obtenidas:





In [22]:



    


def flin(X):
    return X*Coeff[1]+Coeff[0]

def fe(X):
    return exp(CoeffExp[0])*exp(CoeffExp[1]*X)

tt = linspace(1, 20, num=1000)
p8, = plot(t, O, 'ko-')
p9, = plot(tt, flin(tt), 'r')
p0, = plot(tt, fe(tt), 'g')
xlim([1,20])
ylim([0,500000])
legend([p8, p9, p0], ["Datos", "Ajuste lineal", "Ajuste exponencial"], loc=2)
show()

Se observa que dentro del rango la estimación exponencial se comporta muy bien, pero proyectando hacia los 10 años siguientes se aleja demasiado de la estimación lineal (eso asumiendo que la lineal sea la más adecuada). Ignorando esto, se muestran los valores de oferta para los diez años siguientes:

• Lineal





In [23]:



    


for i in range(11,20): display(Math(r'f({}) = {}'.format(i, int(around(flin(i), 0)))))



    


















    







$$f(11) = 225333$$









    







$$f(12) = 237576$$









    







$$f(13) = 249818$$









    







$$f(14) = 262061$$









    







$$f(15) = 274303$$









    







$$f(16) = 286545$$









    







$$f(17) = 298788$$









    







$$f(18) = 311030$$









    







$$f(19) = 323273$$

• Exponencial





In [24]:



    


for i in range(11,20): display(Math(r'f({}) = {}'.format(i, int(around(fe(i), 0)))))



    


















    







$$f(11) = 238611$$









    







$$f(12) = 258627$$









    







$$f(13) = 280322$$









    







$$f(14) = 303838$$









    







$$f(15) = 329326$$









    







$$f(16) = 356952$$









    







$$f(17) = 386895$$









    







$$f(18) = 419350$$









    







$$f(19) = 454528$$

7. Ajustar a la función:

$$Y = Ax^M$$

Los siguientes datos:

</tr> </table> </center>

En este caso, las variables usadas en la matriz son $X = log x$ y $Y = log y$, y la equivalencia es:

$$M\overline{X} + B \longrightarrow y = 10^B X^M$$

$x$

</td>

$y$

$1$

$7$

$2$

$30$

$3$

$90$

$4$

$170$

$5$

$290$

$6$

$450$

$7$

$650$

 





In [25]:



    


x4 = array([1, 2, 3, 4, 5, 6, 7])
y4 = array([7, 30, 90, 170, 290, 450, 650])
logx4 = log10(x4)
logx4sq = logx4**2
logy4 = log10(y4)
logx4logy4 = logx4*logy4

plot(x4,y4,'k-')

A4 = array([[len(logx4), sum(logx4)],[sum(logx4), sum(logx4sq)]])
B4 = array([sum(logy4), sum(logx4logy4)])
print (A4)
print (B4)



    


















    






[[ 7.          3.70243054]
 [ 3.70243054  2.48900912]]
[ 14.43543459   8.88286304]










    
























In [26]:



    


C4 = linalg.solve(A4,B4)
C4



    


















    
Out[26]:








array([ 0.81875802,  2.35092283])

La función encontrada es:

$$y \approx 10^{0.8188}\,x^{2.3509} \approx 6.59\,x^{2.3509}$$





In [27]:



    


def f4(X):
    return 10**(C4[0])*(X**(C4[1]))

xx4 = linspace(1, 7, num=1000)
pa, = plot(x4, y4, 'ko-')
pb, = plot(xx4, f4(xx4), 'r')
xlim([1,7])
ylim([0,700])
legend([pa, pb], ["Datos", "Ajuste"], loc=2)
show()

8. Ajustar a la función:

$$y = Ce^{Ax}$$

Los siguientes datos:

</tr> </table> </center

$x$

</td>

$y$

$1$

$1.2379$

$1.5$

$0.7598$

$2$

$0.7612$

$3.2$

$0.4342$

$5$

$0.2123$

 





In [28]:



    


x5 = array([1, 1.5, 2, 3.2, 5])
y5 = array([1.2379, 0.7598, 0.7612, 0.4342, 0.2123])
x5sq = x5**2
lny5 = log(y5)
lny5x5 = lny5*x5

plot(x5,y5,'k-')

A5 = array([[len(x5), sum(x5)],[sum(x5), sum(x5sq)]])
B5 = array([sum(lny5), sum(lny5x5)])
print (A5)
print (B5)



    


















    






[[  5.    12.7 ]
 [ 12.7   42.49]]
[ -2.71814772 -11.16272657]










    
























In [29]:



    


C5 = linalg.solve(A5,B5)
C5



    


















    
Out[29]:








array([ 0.51353657, -0.41620713])

La función encontrada es:

$$y \approx e^{0.5135}\,e^{-0.4162x}$$





In [30]:



    


def f5(X):
    return exp(C5[0])*exp(C5[1]*X)

xx5 = linspace(1, 5, num=1000)
pc, = plot(x5, y5, 'ko-')
pd, = plot(xx5, f5(xx5), 'r')
xlim([1,5])
ylim([0,1.5])
legend([pc, pd], ["Datos", "Ajuste"], loc=1)
show()



    


















    
























In [ ]:

Content source: queq/just-stuff

Similar notebooks:

notebook.community | gallery | about