Universidad Nacional Autónoma de Honduras

Facultad de Ciencias / Depto. Matemática Aplicada

MM412 Análisis Numérico / Kerin Cardona

Ejercicios de Laboratorio I

Nombre y cuenta:Carlos Eduardo Lopez Alvarado 20151002214

Codigo de honor:
Me comprometo a responder de forma individual, sin ayuda, sin comunicación con mis
compañeros o alguien más. Comprendo que yo soy la única persona que sabrá si estoy
respondiendo éticamente los ejercicios y, por lo tanto, decido ser ético al presentar mi trabajo.

Instrucciones
Resolver de forma clara y ordenada los siguientes ejercicios. Debe agregar la cantidad de
codigo y programas necesarios.

Fecha de entrega : Miercoles 27 de febrero

Valor : 14%

import sys

import numpy

from math import factorial

import sys

from numpy import exp

from sympy import*
import numpy as np

import sympy as sp

1. La definición matemática de la función exponencial es

2 n ∞ n
x x x
x
exp(x) = e = 1 + x + +. . . + +. . . = ∑ .
2! n! n!
n=0

a. Para x 2.9, calcular ex utilizando la función exp del módulo numpy.

=

b. Escriba un programa que utilice la definición anterior para calcular el mismo valor
anterior.

c. Encuentre el valor de n para que ambos métodos den una respuesta con los mismos 16
dígitos.

inciso a

y=exp(2.9)

18.1741453694431

inciso b

a=0
for i in range (0,30,1):
 a=a+(2.9**i)/factorial(i)
 

print(a)

18.1741453694431

y , a

(18.1741453694431, 18.1741453694431)

round(y,14)

18.1741453694431
z=0
for i in range (0,26,1):
  z=z+(2.9**i)/factorial(i)
  print(z),14

1.00000000000000
3.90000000000000
8.10500000000000
12.1698333333333
15.1168375000000
16.8260999166667
17.6522434180556
17.9945028686309
18.1185719194645
18.1585497247331
18.1701432882610
18.1731997731911
18.1739384237159
18.1741031996022
18.1741373317500
18.1741439306320
18.1741451266793
18.1741453307109
18.1741453635827
18.1741453685999
18.1741453693275
18.1741453694279
18.1741453694412
18.1741453694428
18.1741453694430
18.1741453694431

inciso c

round (y,14) ,z

(18.1741453694431, 18.1741453694431)

2. La operación de Resistence Temperature Detector (RTD) está basada en el hecho de que la

resistencia eléctrica del material detector cambia con la temperatura.

Para el Nickel, que usualmente se utiliza en tales detectores, la resistencia, RT ,a

temperatura T (∘ C ) como una función de la temperatura está dada por
2 4 6
RT = R0 (1 + AT + BT + CT + DT )

donde R0 es la resistencia del detector a 0∘ C y A = 5.485 × 10

−3
,B = 6.65 × 10
−6
,
yD son constantes.
−11 −17
C = 2.805 × 10 = −2 × 10

Considere un detector con R0 = 100Ω y determine la temperatura cuando la resistencia

es 300Ω. Utilice el método de bisección con un nivel de tolerancia de 1 × 10 −6 .

codigo de biseccion
import numpy as np

def bisection(a,b,fun,tol,kmax):
    
    x=np.array([ a , (a+b)/2 , b ])
    fx = fun(x)
    
    if fx[0]*fx[2] > 0 : 
        print (' El signo de f debe ser diferente en los extremos del intervalo') 
    
    elif fx[0] == 0 :
        zero , res , niter  = a , 0, 0 
        
    elif fx[2] == 0 :
        zero , res , niter  = b , 0, 0
        
    else : 
        
        niter = 0
        I = ( b - a )*0.5
        
        while I >= tol and niter < kmax :
            niter = niter + 1
            
            if  fx[0]*fx[1] < 0 :
                x[2] = x[1]
                x[1] = x[0]+(x[2] - x[0])*0.5
                
                fx = fun(x)
                I = (x[2]-x[0])*0.5
                
            elif fx[1]*fx[2] < 0 :
                x[0] = x[1]
                x[1] = x[0]+(x[2] - x[0])*0.5
                
                fx = fun(x)
                I = (x[2]-x[0])*0.5
                
            else :
                x[1]=np.find(fx==0)
        
        if  niter == kmax and I > tol :
            print ('Se alcanzó el max de iteraciones sin convergencia ')
            
        zero = x[1]
        x = x[1]
        res = fun(x)
        
    return zero, res, niter

A=5.485*10**(-3)
B= 6.65*10**(-6) 
C=2.805*10**(-11)
D=-2*10**(-17)
R0=100
Rt=300

f = lambda T:R0*(1+A*T+B*T**2+C*T**4+D*T**6)-Rt

bisection(200,400,f,1e-6,100)

(260.21868363022804, 7.873809408920351e-07, 27)

3. La ecuación de van der Waals establece la relación entre la presión P (en atm.), volumen
V (en L.), y temperatura T (en K) para un gas real.
2
nRT n a
P = −
2
V − nb V

donde n es el número de moles, R = 0.08206 (L atm)(mol K) es la constante del gas y a

(en L2 atm/mole2) son las constantes del material.

Considere 1.5 moles de nitrógeno (a=1.39 L2 atm/mol2 , b=0.03913 L/mol) a 25 ∘ C

almacenado en un recipiente a presión. Determinar el volumen del recipiente si la presión
es de 13.5 atm, utilice el método de Newton con un nivel de tolerancia de 1 × 10 .
−12

Codigo para metodo de Newton

def newton(fun,dfun,x0,tol,kmax):
    
    x = x0
    fx , dfx = fun(x) , dfun(x)
    niter = 0
    diff = tol + 1
    
    while diff >= tol and niter < kmax :
        niter = niter + 1;
        diff = - fx / dfx
        x = x + diff
        diff = abs ( diff )
        
        fx , dfx = fun(x) , dfun(x)
        
    if  niter == kmax and diff > tol :
        print  (' Newton alcanzo max de iteraciones ')
        
    zero , res = x , fx
    
    return zero, res,niter
R=0.08206
n=1.5
a=1.39
b=0.03913
t=25+273
P=13.5

V=Symbol('V')

fa= lambda V: (n*R*t)/(V-n*b)-(n**2*a)/V**2-P

fa(3)

-1.3765661798079432

fa1=diff(fa(V),V,1)

fa1

36.68082 6.255
− 2 +
V3
​ ​

(V − 0.058695)

fa11= lambda V: -36.68082/(V-0.058695)**2+6.255/V**3

newton(fa,fa11,5,1e-12,100)

(2.691599427783978, 0.0, 9)

4. Cuando se calcula el pago de una hipoteca, la relación entre el monto del préstamo, Loan,
el pago mensual, Mpay , la duración del préstamo en meses Months y la tasa anual de
interés Rate , está dada por la ecuación
Loan ⋅ Rate
M P ay =

1
12 (1 − )
M onths
Rate
(1+ )
12

Determinar la tasa de un préstamo de 300,000 a 20 años si el pago mensual es de

1,684.57.

Realice las primeras 10 iteraciones del método de bisección y tome ese valor como
aproximación inicial del método de Newton para alcanzar una presición del orden
.
−12
1 × 10

Mpay=1684.57
Months=12*20
Loan=300000

Tasa = lambda t:(Loan*t)/(12*(1-1/(1+t/12)**Months))-Mpay

bisection(0.03,2,Tasa,1e-12,10)

Se alcanzó el max de iteraciones sin convergencia

(0.0309619140625, -6.294266783107105, 10)

t=Symbol('t')

Tasa1=diff(Tasa(t),t,1)

Tasa1

72000000t 300000
− 2 +
12 − t 12 240
​ ​

(12 − ) ( 12t + 1)
241
12 ( 12 +1)
​

240
​

( 12t +1)
​ ​

Tasa11=lambda t:-((72000000)*t)/((12-12/(t/12+1)**240)**2*(t/12+1)**241)+300000/(12-12/(t/

Tasa11(0.08)

18670.66662169849

newton (Tasa,Tasa11,0.0309619140625,1e-12,100)

(0.03137843640854579, -3.2514435588382185e-11, 3)

5. La resistividad ρ del silicio dopado se basa en la carga q en un electrón, la densidad

electrónica n y la movilidad electrónica μ .

La densidad de electrones se da en términos de la densidad de dopaje N y la densidad

intrínseca de portadores ni .

La movilidad de los electrones es descrito por la temperatura T , la temperatura de

referencia T0 y la movilidad de referencia μ0 .

Las ecuaciones requeridas para calcular la resistividad son

1
ρ =
qnμ

donde
−−−−−−−−
1 2 2
n = (N + √ N + 4n )
2 i
 −2.42
T
μ = μ0 ( )
T0

Determinar N, dado T0 ,
= 300K T = 1000K μ0 = 1300cm (V s)
2 −1
,
, , y un ρ objetivo de 6 × 10 .
−19 9 −3 6
q = 1.6 × 10 C n i = 6.21 × 10 cm Vscm\C

Utilice el método de Newton para alcanzar una presición del orden 1 × 10 −8 .

ejercicio # 5

T0=300
u0=1300
q=1.6*10**(-19)
ni=6.21*10**9
p=6*10**6
T=1000

u=u0*(T/T0)**(-2.42)

70.56295972119374

n=1/(p*q*u)

14762230365.371134

Den=lambda N:1/2*(N+sqrt(N**2+4*(ni)**2))-n

N=Symbol('N')

Den1=diff(Den(N),N,1)

Den1

4.02576489533011 ⋅ 10−11 N
+ 0.5
6.48271319698891 ⋅ 10−21 N 2 + 1
​

Den11=lambda N:(4.02576489533011*10**-11*N)/(sqrt(6.48271319698891*10**-21*N**2+1))+0.5

Den11(1)

0.500000000040258
newton (Den,Den11,12149881211.78,1e-6,10000)

Newton alcanzo max de iteraciones

(12149881211.7864, -1.90734863281250e-6, 10000)

6. Dos fluidos a diferentes temperaturas entran en un mezclador y salen a la misma

temperatura. La capacidad calorífica del líquido A esta dada por
−2 −6 2
cp (A) = 3.381 + 1.804 × 10 T − 4.300 × 10 T

y la capacidad de calor del fluido B esta dada por

−1 −5 2
cp (B) = 8.592 + 1.290 × 10 T − 4.078 × 10 T

donde cp tiene unidades cal/mol K, y T en K. Notar que

T2

ΔH = ∫ cp dT
T1

El fluido A entra al mezclador a 400 ∘ C y el fluido B a 600 ∘ C . Entra dos veces la cantidad
del fluido A con respecto al fluido B.

Encontrar las temperaturas de salida de ambos fluidos. Utilice un nivel de tolerancia de

.
−6
1 × 10

T=sp.Symbol('T')

cpA= 3.381+1.804*10**-2*T-4.300*10**-6*T**2

cpB=8.592+1.290*10**-1*T-4.078*10**-5*T**2

T=sp.Symbol('T')

cpA1=sp.integrate(cpA,(T,673,T))

cpA1

−1.43333333333333 ⋅ 10−6 T 3 + 0.00902T 2 + 3.381T − 5923.92216896667

cpB1=sp.integrate(cpB,(T,873,T))

cpB1

−1.35933333333333 ⋅ 10−5 T 3 + 0.0645T 2 + 8.592T − 47613.96689958

FU=lambda T: -1.43333333333333*10**-6*T*3+0.00902*T**2+3.381*T-5923.92216896667-2*(-1.3593
FU(700)

30176.735286859974

T=Symbol('T')

FU1=diff(FU(T),T,1)

FU1

8.15599999999998 ⋅ 10−5 T 2 − 0.23996T − 13.8030043

FU11=lambda T:8.15599999999998*10**-5**T*2-0.23996*T-13.8030043

bisection(672,900,FU,1e-6,100)

(897.0445343703032, 1.225577580044046e-05, 27)

7. Determinar la velocidad de sedimentación de las partículas a través de los líquidos es de

gran importancia de muchas áreas de la ingeniería y la ciencia.

Dichos cálculos dependen del régimen bajo representado por el número de Reynolds
adimensional,
ρdv
Re =
μ

donde ρ = densidad del fluido (kg/m3 ), d = diametro de la particula (m), v = velocidad

de la particula (m/s) y μ = viscosidad dinamica del fluido (N s/m2 )

En condiciones laminares (Re < 0.1 ),la velocidad de sedimentación de una partícula
esférica se puede calcular con la siguiente fórmula basada en la ley de Stokes,
g ρs − ρ
2
v = ( )d
18 μ

donde g es la constante gravitacional (9.81 m/s2 ) y ρs es la densidad de la particula (kg/

m
3
).

Para condiciones de turbulencia ( número de Reynolds mas altos),un enfoque alternativo

se basa en la siguiente formula
− −−−− −− − −
4g(ρs − ρ)
v = √
3C D ⋅ ρ

donde CD es el coeficiente de arrastre,el cual depende del número de Reynolds de la

forma
 24 3
CD = + + 0.34
−
−−
Re √Re

a. Encuentre una expresion para f (v) = 0 .

b. Utilice el metodo de Newton con ϵ , para determinar v para una particula de

−6
= 1 × 10

hierro esferica asentada en agua, esto es d = 200μ m, ρ = 1 g/cm3 , ρs = 7.874g/cm

3

yμ = 0.014 g/(cmxs).

Para el valor inicial utilizar la primera ecuacion dada para v.

c. Para el resultado encontrado en (b.), calcular el numero de Reynols y el coeficiente de

arrastre y utilizar lo anterior para confirmar que el regimen de flujo no es la minar.

g=9.81
E=1*10**(-6)
d=200*1e-6
p=1000
ps=7874
u=0.0014

a. Encuentre una expresion para f(v)=0

Re=lambda v: (p*d*v)/u

(p*d*v)/u

142.857142857143v

v0=(g/18)*((ps-p)/u)*d**2

v0

0.10703799999999998

Cd= lambda v:24/(Re(v))+3/sqrt(Re(v))+0.34

x= lambda v:v- sqrt((4*g)*(ps-p)/(3*Cd(v)*p))

v- sqrt((4*g)*(ps-p)/(3*Cd(v)*p))

1
v − 519.360914971468 504.0 752.994023880668
1020.0 + +
​

v v
​ ​

v=Symbol('v')
vel1=diff(x(v),v,1)

vel1

0.000249596748832885 ( 504.0
v2
​ + 376.497011940334
3 ​ ) (1020.0 + 504.0
v
​ + 752.994023880668
v
​) 1020.0+
v2
​

−
​

2
(1 + 0.494117647058824
v
​ + 0.738229435177125
v
) ​

a=0.000249596748832885
b=lambda v: (504/v**2)+(376.497011940334)/(v**1.5)
c=lambda v:1020+504/v+(752.994023880668)/(sqrt(v))
d=lambda v:1+(0.494117647058824/v)+(0.738229435177125/sqrt(v))

veli11= lambda v: 1-(a*b(v)*c(v)*(sqrt(1/c(v))))/(d(v)**2)

newton (x,veli11,v0,10**-6,100)

por el metodo de newton no es posibre resolver la funcion esto porque en x1 queda un numero
negativo y en la funcion hay raices

Metodo de biseccion

bisection(13,18,x,1e-6,100)

(14.684988141059875, 4.910598754293005e-07, 22)

Re=(p*d*14.6849)/u

Re

2097.842857142857

el regimen de flujo no es la minar.