GUIA Programacion Aplicada
GUIA Programacion Aplicada
GUIA Programacion Aplicada
FACULTAD DE INGENIERIA
PROGRAMACION APLICADA
PROGRAMA ANALITICO
I. JUSTIFICACION
La computación es una ciencia en constante evolución y crecimiento. Su aplicación
prácticamente no tiene límites, se emplea en todas las áreas, por supuesto que en la
Ingeniería Petrolera tiene un amplio desarrollo. Si bien en la actualidad existe
software desarrollado por empresas especializadas, la Universidad en general y el
estudiante en particular no tiene acceso a ellos, por su alto costo y limitaciones
presupuestarias; además este software es cerrado y solo permite el acceso a nivel de
usuario, sin poder ver las técnicas y el desarrollo efectuado. De ahí la necesidad de
formar a los futuros profesionales en el manejo de una herramienta de programación
como el Visual Basic para resolver problemas específicos relacionados con la
Industria Petrolera
II. OBJETIVOS
Una parte de la materia se dicta mediante clases magistrales, otra mediante clases
participativas con planteamientos de problemas y resolución individual y por grupos.
En ambos casos se utiliza proyección de métodos utilizados y solución de problemas
paso a paso. Se complementa la enseñanza con el uso de paquetes de aplicación en
la Industria Petrolera y la asignación y ejecución de proyectos individuales por
alumno con la correspondiente exposición, prueba y defensa de los mismos.
V. CRITERIOS DE EVALUACION
Primer Parcial 15
Segundo Parcial 15
Prácticas 15
Ayudantía 15
Laboratorio 20
Asistencia 5
Examen Final 15
Total 100
VI. CRONOGRAMA
PROGRAMA DE TRABAJO
Mes 1 Mes 2 Mes 3 Mes 4 Mes 5
Ca-pí-
Tema Horas Semanas
tulo
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
1 Consideraciones Generales 10
2 Ecuaciones no Lineales 10
3 Sistemas de Ecuaciones 12
1a3 Primer Examen Parcial 2
4 Aproximación Funcional e Interpolación 10
5 Integración Numérica 10
6 Diferenciación Numérica 4
4a6 Segundo Examen Parcial 2
7a 9 Exposición y Defensa de Proyectos 14
1a9 Examen Final 2
1a9 Examen Segundo Turno 2
Registro Final de Notas 2
Total 80
VII. BIBLIOGRAFIA
4.1 Introducción 1 de 21
4.2 Aproximación Polinomial Simple e Interpolación 2 de 21
4.3 Polinomios de Lagrange 5 de 21
4.4 Diferencias Divididas 9 de 21
4.5 Aproximación Polinomial de Newton 12 de 21
4.6 Aproximación Polinomial con Mínimos Cuadrados 15 de 21
4.7 Aproximación Multilineal con Mínimos Cuadrados 20 de 21
5.1 Introducción 1 de 17
5.2 Métodos de Newton – Cotes 1 de 17
5.2.1 Método Trapezoidal 2 de 17
5.2.2 Método de Simpson 4 de 17
5.3 Caso General 6 de 17
5.4 Métodos Compuestos de Integración 7 de 17
Ing. Hermas Herrera Callejas Página: 2 de 3
PET230 – Programación Aplicada Contenido de la Asignatura
6.1 Introducción 1 de 8
6.2 Derivación con Polinomios de Lagrange 2 de 8
Existe una gama de paquetes de software desarrollados para las distintas áreas de la
industria petrolera, se dará un vistazo a los más utilizados en las distintas funciones sobre la
base del acuerdo mínimo de estandarización logrado entre un grupo importante de empresas
con la finalidad de buscar un lenguaje común.
1.1.2 Exploración
1.1.2.3 PETCOM.- Aplicación que corre en una PC bajo el sistema operativo DOS o
Windows. Permite la captura de información desde las cintas magnéticas de sísmica y se
utiliza en el análisis petrofísico y de registros de pozos, generando los reportes y gráficos
relacionados a este campo.
1.1.3 Perforación
1.1.3.1 DES-II (Drilling Expert System II).- Permite analizar y optimizar las tareas
operativas diarias de perforación de pozos. Corre en PC’s con Windows. Tiene siete módulos
independientes pero relacionados entre sí tales como el de perforación direccional, diseño de
perforación, optimización de la perforación, sistema experto de perforación, hidráulica y
cementación, ingeniería de lodos, control del pozo y detección de presiones anormales.
1.1.4 Producción
1.1.4.2 AUTOMATE.- Programa que permite efectuar un análisis de los distintos tipos de
presiones en los pozos productores y los reservorios petrolíferos.
1.1.5.4 SimBest II.- Para simulación de reservorios de la Scientific Software Inc, emplea
el paquete ESPIDO (Equation Solution Program based on an Incomplete Direct Method
acelerated via Orthomin), que usa la eliminacion de Gauss para problemas pequeños y
SOR (successive over relaxation) para los grandes.
1.1.5.6 CHEMCAD.- Aplicación que corre en una PC bajo el sistema operativo DOS o
Windows. Tiene como fin la preparación de reportes y gráficas relacionadas a la composición
de hidrocarburos de gas y petróleo producido en campos petrolíferos.
Ing. Hermas Herrera Callejas Página: 2 de 20
Programación Aplicada Capítulo 1 – Consideraciones Generales
1.1.6.2 FOAS.- Sistema considerado como el cerebro central que se encarga del proceso
contable de la compañía petrolera. Permite captar la información de los asientos contables
de todo el movimiento económico generado en la empresa y una vez procesada muestra los
resultados de la operación, emite los estados financieros de la compañía y los libros legales
correspondientes. Permite generar reportes financieros detallados por proyectos o
reservorios así como consolidados de la compañía. Esta característica permite el proceso
contable de multicompañías. Como funciones adicionales se tienen los siguientes módulos
que complementan el FOAS
1.1.6.2.3 CUENTAS POR COBRAR.- Aplicación en línea que permite captar la información
de las facturas emitidas por la compañía por conceptos de venta de gas y petróleo. También
permite introducir las cobranzas o los pagos parciales y llevar un control de las deudas,
calculando intereses por montos vencidos así como clasificar las facturas vencidas por
períodos de 30 días, 60 días, 90 días, 120 días, 180 días, o más de 180 días.
1.1.6.2.4 CUENTAS POR PAGAR.- Aplicación en línea que permite hacer un control y
seguimiento a las facturas pendientes de pago y las canceladas. Ayuda en la planificación de
los pagos y a controlar la no-duplicación de estos últimos, especialmente aquellos casos de
contratos que implican el pago mensual por servicios.
1.1.6.3 OPICS.- Aplicación interactiva en línea para el control de inventarios. Permite llevar
un control de los materiales de la compañía en los distintos almacenes, así como hacer un
seguimiento a las ordenes de compra y generar los asientos contables en forma automática
para cada transacción que debe ser contabilizada.
1.1.6.4 HRIAS (Human Resources Information Application System).- Permite contar con
información de recursos humanos en línea. Incluye los siguientes módulos:
1.2 Algoritmos
Término genérico para nombrar las instrucciones del programa, utilizadas en dos
sentidos generales derivados del diagrama de flujo.
Constituida por todos los documentos que se elaboran en cada una de las etapas
del análisis, diseño y desarrollo de la aplicación, es muy importante para facilitar su
mantenimiento y obtener un mayor rendimiento.
Denominamos documentación interna al contenido del propio programa fuente.
Debe incluir los comentarios explicativos suficientes que posibiliten su comprensión y
actualización.
Asimismo, se debe utilizar un código autodocumentado; es decir, debe ser escrito
de una forma clara y legible.
La documentación externa la forman el resto de documentos que se acompañan
con el programa sin formar parte de él. Dentro de ellos están los manuales internos del
sistema que incluyen detalles de técnicas y diseños de bases de datos, programas, etc,
que constituyen la aplicación; los manuales del usuario que describen la manera en que el
usuario puede obtener mejor provecho de la aplicación así como una explicación de los
reportes y la información que proporciona. También forma parte de este tipo de
documentación los manuales en línea de las aplicaciones así como los textos de ayuda a
los que el usuario puede acudir
1.- Diagrama de Flujo para calcular el área de 2.- D.F. para hallar el cociente y el residuo
un triángulo de A\B enteros
Inicio
Inicio
Def A, B, C, D
Def b, h
Leer A, B
Leer b, h
C = A Mod B
A = b*h D = A\B
Imprimir A Imprimir C, D
Fin Fin
3.- D.F. para hallar la longitud de una circunferencia y el área del círculo
4.- D.F. para convertir metros en Km y cm
5.- D.F. para convertir Kb a Gb, Mb y bytes
6.- Hallar el mayor de 3 números diferentes
7.- Hallar el mayor y el menor de 3 números diferentes
8.- Hallar el mayor y el menor de 3 números cualesquiera
9.- Determinar si un número es par o impar
10.- Desplegar los números enteros de N hasta M
11.- Imprimir la tabla del 4
12.- Hallar la suma de los primeros 10 números pares
13.- Hallar la suma de los primeros 10 números impares
14.- Hallar los cuadrados de los primeros 10 números pares
15.- Determinar si el número introducido es positivo o negativo
16.- Hallar el factorial de un número entero positivo
17.- Crear el vector I = 1, 2, 3, …10
Inicio
Def I, V(I)
I = 1 … 10
V(I) = I
Imprimir V
Fin
Inicio
A
Fin ? Si Fin B
No
A No Ejecutar ? DIVE>2 ? Si
Si No
Leer N K=K+1
A P(K) = J
I
J Mod I = 0 ?
No
Si A
DIVE = DIVE + 1
Inicio
Def I, V(I), N
Inicio
Leer N
Def I, V(I), N
N>0? Leer N
I=1…N I=1…N
V(I) = 2 ^ I V(I) = 0
I I
Imprimir V Imprimir V
Fin Fin
Inicio Inicio
Leer N Leer N
I=1…N I=1…N
V(I) = N - I V(I) = I * I
I I
Imprimir V Imprimir V
Fin Fin
23.- Crear el vector de N elementos donde 24.- Sea N un Nro entero. Hacer un D.F.
c/elemento a partir del 3ro sea la suma para invertir sus dígitos (Ej, 3457 a 7543)
de los dos anteriores y V(1)=1 V(2)=2
Inicio Inicio
Leer N Leer N
A=N
N>2?
N1 = 0
V(1) = 1 V(2) = 2
A>0?
I=3…N
Dig = A Mod 10
V(I) = V(I-1) + V(I-2) N1 = N1 * 10 + Dig
A = A Div 10
I
Imprimir V Imprimir N, N1
Fin Fin
25.- Generar la serie de Fibonacci para 26.- Crear un vector con N elementos,
valores menores a N (0,1,1,2,3,5,8,13…) luego obtener el máximo y su posición
Inicio Inicio
A Leer N
Fin ? Fin
N>0?
Ejec? A
I=1…N
Leer N
Leer X
N>3 ? A
V(I) = X
F(1) = 0 F(2) = 1 I
Max = V(I) K = I
I = 3… N
I=1…N
F(I) = F(I-1) + F(I-2)
I Fin
Leer N
N>0?
I=1…N
Leer X
V(I) = X
I = 1 … N-1
J = 1 … N-I
Aux = V(J)
V(J) = V(J+1)
V(J+1) = Aux
I=1…N
Imprimir V(I)
Fin
Inicio Inicio
Leer N Leer N
N>0? N>0?
I=1…N I=1…N
Leer X Leer X
A(I) = X A(I) = X
I I
I=1…N I=1…N
Leer X Leer X
B(I) = X B(I) = X
I I
I=1…N I=1…N
I I
I=1…N I=1…N
I I
Fin Fin
30.- Crear una matriz de N filas por N 31.- Crear una matriz de N filas por M co-
columnas cuyos elementos sean ceros lumnas cuyas filas pares sean unos y las
impares sean ceros
Inicio
Inicio Leer N, M
Leer N
I=1…N
N>0 ? J=1…M
I=1…N
I Mod 2 = 0?
J=1…N A(I, J) = 0
A(I, J) = 0 A(I, J) = 1
J J
I I
Imprimir A Imprimir A
Fin Fin
32.- Crear una matriz N por N con la 33.- Crear una matriz N por M con
diagonal principal igual a 1 numeración correlativa ascendente
Inicio
Inicio
Def A(I, J), N, M, I, J, C
Def A(I, J), N, I, J
Leer N, M
Leer N
N>0 y M>0?
N>0?
C=0
I=1…N
I=1…N
J=1…N
J=1…M
I = J?
C=C+1
J J
I I
Imprimir A Imprimir A
Fin Fin
34.- Construir una matriz N por N con N 35.- Construir la matriz N por N 1 2 3 4
N impar y mayor a 2. Calcular la suma 2 4 2 2456
de la siguiente manera (suma = 17) 1 2 3 3567
279 4678
Inicio
Inicio
Def A(I, J), N, I, J
Def A(I, J),N,I,J,C,S,K
Leer N
Leer N
N>1 ?
N>2 y N Mod 2=1
I=1…N
I=1…N
A(1, I) = I
J=1…N A(I, 1) = I
Leer C I
A(I, J) = C I=2…N
J J=2…N
I A(I, J) = I + J
J
S=0 K = N\2 + 1
I=1…N I
I=1…N
S = S + A(I, K)
S = S + A(K, I) J=1…N
I
Imprimir A(I, J)
S = S – A(K, K)
J
Imprimir S
I
Fin
Fin
Inicio
Leer N
N>2?
F=1
C=N
R=0
B
J = F…C
R=R+1
A(F, J) = R
A
J=C-1…F+1, -1
J = F+1…C
R=R+1 R=R+1
A(J, C) = R A(J, F) = R
J J
J=C-1…F, -1 F=F+1
C=C-1
R=R+1
B R>NxN
A(C, J) = R
J Imprimir A(I, J)
A Fin
37.- Formar la matriz zigzag N por N 38.- Convertir un número decimal a binario
para N > 2
Inicio
Leer N
N>2?
C=0 Inicio
I = 1…N
Def A(I),N,M,I,J
J = 1…N
Leer M
C=C+1
A(I, J) = C M>0?
J
N=M
L=I+1 I=0
L>N? I=I+1
A(I) = N Mod 2
K= N…1, -1
N = N\2
C=C+1
N=0?
A(L, K) = C
J= I…1, -1
K
Imprimir A(J)
I
K
Imprimir A
Fin
Fin
39.- Sumar los elementos de cada fila y cada columna de una matriz N por M
Inicio
Leer N,M
A
N>1 M>1?
J = 1…M
I = 1…N
C(I) = 0
J= 1…M
I = 1…N
Leer R
C(J) = C(J)+A(I,J)
A(I, J) = R
I
J
J
I
I = 1…N
I = 1…N
Imprimir F(I)
F(I) = 0
I
J= 1…M
J = 1…M
F(I) = F(J)+A(I,J)
J Imprimir C(J)
I J
A Fin
Inicio
Inicio
Def A(I,J),B(I,J),C(I,J),N,M,I,J
Def A(I,J),T(I,J),N,I,J,M
Leer N,M
Leer N,M
N>1 M>1?
N>1 M>1?
I = 1…N
I = 1…N
J= 1…M
J= 1…M
Leer R
Leer R
A(I, J) = R
A(I, J) = R
J A
J
I
I I = 1…N
I= 1…N
J= 1…M J= 1…M
J= 1…M
I= 1…N C(I, J) = A(I,J)
+ B(I,J)
Leer R
T(J, I) = A(I, J)
Imprimir
Imprimir T(I,J) B(I, J) = R C(I, J)
J J J
I I I
Fin A Fin
Inicio
A
Leer M,N,O
I= 1…M
M>1 N>1 O>1?
J = 1…O
I= 1…M
C(I, J) = 0
J = 1…N
K= 1…N
Leer R
C(I, J) = C(I,J)+A(I,K)*B(K,J)
A(I, J) = R
K
J
J
I
I
I = 1…N
I= 1…M
J = 1…O
J = 1…O
Leer R
Imprimir
B(I, J) = R C(I, J)
J J
I I
A Fin
Ejemplos:
1) La ecuación: cos x – x = 0 se puede transformar en: cos x = x.
2) La ecuación: tan x – e-x = 0 se puede transformar en: x + tan x – e-x = x.
Dada la aproximación xi, la siguiente iteración se calcula con la fórmula:
xi 1 g ( xi )
Supongamos que la raíz verdadera es xr, es decir,
xr g ( xr )
Restando las últimas ecuaciones obtenemos:
x r xi 1 g ( x r ) g ( xi )
Por el Teorema del Valor Medio para derivadas, sabemos que si g(x) es continua en
g (b) g (a )
[a, b] y diferenciable en (a, b) entonces existe ξ Є (a, b) tal que g ' ( ) .
ba
En nuestro caso, existe ξ en el intervalo determinado por xi y xr y tal que:
g ( x r ) g ( xi )
g ' ( )
x r xi
De aquí tenemos que:
g ( x r ) g ( xi ) g ' ( )( x r xi )
O bien,
x r xi 1 g ' ( )( x r xi )
Tomando valor absoluto en ambos lados,
| x r xi 1 || g ' ( ) || x r xi |
Observe que el término |xr–xi+1| es precisamente el error absoluto en la
(i+1)ésima iteración, mientras que el término |xr-xi| corresponde al error absoluto en
la i-ésima iteración.
Por lo tanto, solamente si |g’(ξ)| < 1, entonces se disminuirá el error en la
siguiente iteración. En caso contrario, el error irá en aumento.
En resumen, el método de iteración del punto fijo converge a la raíz si |g’(x)| <
1 para x en un intervalo [a, b] que contiene a la raíz y donde g(x) es continua y
diferenciable, pero diverge si |g’(x)| > 1 en dicho intervalo.
Analicemos nuestros ejemplos anteriores:
En el ejemplo 1, g(x) = cos x y claramente se cumple la condición de que
|g’(x)| < 1. Por lo tanto el método sí converge a la raíz.
En el ejemplo 2, g(x) = x+tan x– e-x, en este caso |g’(x)| = |1 + sec2x + e-x| > 1.
Por lo tanto, el método no converge a la raíz.
Para aclarar el uso de la fórmula veamos dos ejemplos:
Ejemplo 1
Usar el método de iteración del punto fijo para aproximar la raíz de f(x) = cos x – x.
comenzando con x0 = 0 y hasta que |Єa| < 1%.
Solución
Como ya aclaramos anteriormente, el método sí converge a la raíz. Aplicando la
fórmula iterativa tenemos,
x1 = g(x0) = cos 0 = 1
Con un error aproximado de 100%
Aplicando nuevamente la fórmula iterativa tenemos,
x2 = g(x1) = cos 1 = 0.540302305
Y un error aproximado de 85.08%.
Intuimos que el error aproximado se irá reduciendo muy lentamente. En
efecto, se necesitan hasta 13 iteraciones para lograr reducir el error aproximado
menor al 1%. El resultado final que se obtiene es:
x13 = 0.7414250866
Con un error aproximado igual al 0.78%.
Ejemplo 2
Usar el método de iteración del punto fijo para aproximar la raíz de f(x) = x2 – 5x – ex.
comenzando con x0 = 0 y hasta que |Єa| < 1%.
Solución
Si despejamos la x del término lineal vemos que la ecuación equivale a
x2 ex
x
5
de donde,
x2 ex
g ( x)
5
2x e x
En este caso, tenemos que g ' ( x) . Un vistazo a la gráfica, nos convence
5
que |g’(x)| < 1, para x Є [-1, 1], lo que es suficiente para deducir que el método sí
converge a la raíz buscada.
Aplicando la fórmula iterativa, tenemos:
x1 = g(x0) = -0.2
Con un error aproximado del 100%.
i Xi % de Error f(x) = x2 – 5x – ex
0 0,0000000000 x2 ex 2x e x
1 -0,2000000000 100,000000 x g ( x) g ' ( x)
5 5
2 -0,1557461506 28,414089
3 -0,1663039075 6,348472 xi2 e xi
4 -0,1638263720 1,512293 xi 1
5 -0,1644100640 0,355022 5
De donde vemos que la aproximación buscada es: x5 = -0.164410064
Veremos a continuación un ejemplo del método de Punto Fijo con la siguiente
ecuación:
X3 + X + 16 = 0
Se ve que no converge
Este método, que es un método iterativo, es uno de los más usados y efectivos. El
método de Newton-Raphson no trabaja sobre un intervalo sino que basa su fórmula
en un proceso iterativo. En análisis numérico, el método de Newton-Raphson
Trazamos la recta tangente a la curva en el punto (xi, f(xi)); ésta cruza al eje x
en un punto xi+1 que será nuestra siguiente aproximación a la raíz xr.
Para calcular el punto xi+1, calculamos primero la ecuación de la recta
tangente. Sabemos que tiene pendiente
m = f’(xi)
Y por lo tanto la ecuación de la recta tangente es:
y – f(xi) = f’(xi)(x – xi)
Hacemos y = 0:
- f(xi) = f’(xi)(x - xi)
Y despejamos x:
f ( xi )
x xi
f ' ( xi )
Que es la fórmula iterativa de Newton-Raphson para calcular la siguiente
aproximación:
f ( xi )
xi 1 xi si f ' ( xi ) 0
f ' ( xi )
Ejemplo 1
Solución
En este caso, tenemos que
1
f ' ( x ) e x
x
De aquí tenemos que:
e xi ln( xi ) e xi ln( xi ) xi (e xi ln( xi )) xi ( xi e xi 1 e xi ln( xi ))
xi 1 xi xi xi xi
e xi 1
e xi 1 x i e 1 xi e xi 1
xi xi
i Xi % de Error
0 1,0000000000 xi ( xi e xi 1 e xi ln( xi ))
xi 1
1 1,2689414214 21,194156 xi e xi 1
2 1,3091084033 3,068270
3 1,3097993887 0,052755
De lo cual concluimos que la aproximación obtenida es: x3 = 1.309799389
Ejemplo 2
SOLUCIÓN
En este caso, tenemos que
1
f ' ( x) 1
1 x2
La cual sustituimos en la fórmula de Newton-Raphson para obtener:
arctan( xi ) xi 1
xi 1 xi
1
1
1 xi2
Comenzamos sustituyendo x0 = 0 para obtener:
arctan( x0 ) x0 1
x1 x0 0 .5
1
1
1 x02
0 .5 0
En este caso tenemos un error aproximado de a x100% 100%
0 .5
Continuamos con el proceso hasta lograr el objetivo.
Resumimos los resultados en la siguiente tabla:
i Xi % de Error arctan( xi ) xi 1
0 0,0000000000 xi 1 xi
1
1 0,5000000000 100,000000 1
1 xi2
2 0,5201957728 3,882341
3 0,5202689918 0,014073
Ejemplo 3
Solución
Sea R > 0.Queremos calcular x tal que x R ; elevando al cuadrado x2 = R, o bien:
x2 – R = 0
Esto nos sugiere definir la función f(x) = x2 – R de donde f’(x) = 2x. Al sustituir
estos datos en la fórmula de Newton-Raphson nos da:
x2 R
xi 1 xi i
2 xi
La cual simplificada nos da:
1 R
xi 1 xi
2 xi
Esta fórmula era conocida por los antiguos griegos (Herón).
Para fijar un ejemplo de su uso, pongamos R = 26 y apliquemos la fórmula
obtenida, comenzando con x0 = 5. Resumimos los resultados en la siguiente tabla:
i Xi % de Error
0 5,0000000000 1 R 1 26
1 5,1000000000 1,9607843 xi 1 xi xi 1 xi
2 xi 2 xi
2 5,0990196078 0,0192271
3 5,0990195136 0,0000018
De lo cual concluimos que 26 ≈ 5.0990195136, la cual es correcta en todos
sus dígitos.
La misma idea puede aplicarse para crear algoritmos que aproximen raíces n -
ésimas de números reales positivos.
Observe que cuando el método de Newton-Raphson converge a la raíz, lo
hace de una forma muy rápida y de hecho, observamos que el error aproximado
disminuye a pasos agigantados en cada paso del proceso. Aunque no es nuestro
objetivo establecer formalmente las cotas para los errores en cada uno de los
métodos que hemos estudiado, cabe mencionar que si existen estas cotas que
miden con mayor precisión la rapidez ó lentitud del método en estudio.
Veremos a continuación un ejemplo del método de Newton Raphson, con la
siguiente ecuación: X3 + X + 16 = 0.
Ejemplo 1
Usar el método de la secante para aproximar la raíz de f ( x) e x x , comenzando
2
Solución
Tenemos que f(x0) = 1 y f(x1) = -0.632120558, que sustituimos en la fórmula de la
secante para calcular la aproximación x2:
f ( x1 )( x0 x1 ) f ( xi )( xi 1 xi )
x 2 x1 0.612699837 xi 1 xi
f ( x0 ) f ( x1 ) f ( xi 1 ) f ( xi )
x 2 x1
Con un error aproximado de: a x100% 63.2%
x2
Como todavía no se logra el objetivo, continuamos con el proceso. Resumimos los
resultados en la siguiente tabla:
( e xi xi )( xi 1 xi )
2
i x(i) % Error Aprox
0 0,000000000 xi 1 xi
( e xi 1 xi 1 ) ( e xi xi )
2 2
1 1,000000000 100,00000
2 0,612699837 63,21206 Haciendo operaciones algebraicas se resume a:
xi e xi 1 xi 1e xi
2 2
3 0,653442133 6,23503
4 0,652917265 0,08039 xi 1
e xi 1 e xi xi xi 1
2 2
5 0,652918640 0,00021
Solución
Tenemos los valores f(x0) = 1 y f(x1) = -0.214601836, que sustituimos en la fórmula
de la secante para obtener la aproximación x2
f ( x1 )( x0 x1 ) f ( xi )( xi 1 xi )
x 2 x1 0.823315073 xi 1 xi
f ( x0 ) f ( x1 ) f ( xi 1 ) f ( xi )
x 2 x1
Con un error aproximado de: a x100% 21.46%
x2
Como todavía no se logra el objetivo, continuamos con el proceso. Resumimos los
resultados en la siguiente tabla:
i x(i) % Error Aprox (arctan( xi ) 2 xi 1)( xi 1 xi )
xi 1 xi
0 0,000000000 arctan( xi 1 ) 2 xi 1 1 (arctan( xi ) 2 xi 1)
1 1,000000000 100,00000
Haciendo operaciones algebraicas se llega a:
2 0,823315073 21,46018
3 0,852330280 3,40422
xi arctan( xi 1 ) xi 1 arctan( xi ) xi 1 xi
4 0,853169121 0,09832 xi 1
5 0,853164044 0,00060 arctan( xi 1 ) arctan( xi ) 2 xi 1 2 xi
i xi % de Error xi 1 xi f ( xi )( xi 1 xi )
0 -3,0000000000 f ( xi 1 ) f ( xi )
1 -2,0000000000 50,000000 Reemplazando las funciones y variables:
2 -2,3000000000 13,043478 x x ( xi3 xi 16)( xi 1 xi )
i 1
xi31 xi 1 16 ( xi3 xi 16)
i
3 -2,4029550034 4,284516
4 -2,3871468897 0,662218 Realizando operaciones algebraicas se tiene:
5 -2,3876833053 0,022466 x x 3 x x 3 16 xi 1 16 xi
xi 1 i i 1 3 i 1 3i
6 -2,3876865541 0,000136 xi 1 xi xi 1 xi
7 -2,3876865535 0,000000
El proceso se vuelve a repetir con el nuevo intervalo, hasta que: |Єa| < Єr, es decir,
x actual x previa
x100% r
x actual
Ejemplo 1
Aproximar la raíz de f(x) = e-x – ln x hasta que |Єa| < 1%
Solución
La única raíz de f(x) se localiza en el intervalo [1, 1.5]. Así que este intervalo es
nuestro punto de partida; sin embargo, para poder aplicar el método de bisección
debemos controlar que f(1) y f(1.5) tengan signos opuestos.
En efecto, tenemos que
f(1) = e-1 – ln 1 = e-1 > 0
(Sabemos que e = 2.71828182845905
Mientras que
f(1.5) = e-1.5 – ln (1.5) = -0.18233 < 0
Cabe mencionar que la función f(x) sí es continua en el intervalo [1, 1.5]. Así
pues, tenemos todos los requisitos satisfechos para poder aplicar el método de
bisección. Comenzamos:
Así, vemos que la raíz se encuentra en el intervalo [1.25, 1.375]. Calculamos el punto
medio,
1.25 1.375
xr 3 1.3125
2
Y calculamos el nuevo error aproximado:
x xr 2
a r 3 x100% 4.76%
xr 3
El proceso debe seguirse hasta cumplir el objetivo. Resumimos los resultados
que se obtienen en la siguiente tabla:
Ejemplo 2
Aproximar la raíz de f(x) =arctan x + x - 1 hasta que |Єa| < 1%.
Solución
Como vimos en el ejemplo 2 de la sección anterior, la única raíz de f(x) se localiza en
el intervalo [0, 1]. Para poder aplicar el método de bisección, es importante controlar
que se cumplen las hipótesis requeridas.
Sabemos que f(x) es continua en el intervalo [0, 1], y controlamos que f(0) y f(1)
tengan signos opuestos. En efecto,
f(0) = arctan 0 + 0 – 1 = -1 < 0
Mientras que,
f(1) = arctan 1 + 1 – 1 = 0.7853 > 0
Por tanto, sí podemos aplicar el método de bisección.
Calculamos el punto medio del intervalo [0, 1],
1 0
x r1 0 .5
2
Que es la primera aproximación a la raíz de f(x)
Puesto que f(0.5) y f(1) tienen signos opuestos, entonces la raíz se localiza en el
intervalo [0.5, 1]
En este punto, solo contamos con una aproximación, a saber xr1 = 0.5, que es el
primer punto medio calculado.
Repetimos el proceso, es decir, calculamos el punto medio ahora del intervalo [0.5, 1]
1 0 .5
xr 2 0.75
2
Que es la nueva aproximación a la raíz de f(x). Aquí podemos calcular el primer
error aproximado:
0.75 0.5
a x100% 33.33%
0.75
Puesto que no se cumple el objetivo, continuamos con el proceso.
Evaluamos f(0.75) = arctan(0.75) + 0.75 – 1 = 0.3935 > 0. y hacemos la tabla de
signos:
Puesto que f(0.5) y f(0.75) tienen signos opuestos, entonces la raíz se localiza en el
intervalo [0.5, 0.75]. Calculamos el punto medio,
0.5 0.75
xr 3 0.625
2
Y el nuevo error aproximado:
0.625 0.75
a x100% 20%
0.625
El proceso se debe continuar hasta que se logre el objetivo.
Resumimos los resultados que se obtienen en la siguiente tabla:
0 a32 a12 a31 / a11 a33 a13 a31 / a11 p3 q3 q1 a31 / a11
Simbólicamente los nuevos valores pueden representarse más simplificadamente así:
a11 a13 p1 q1
a12
0 ' ' '
a a 23
22 p 2 = q 2 3.3
0 a '
a33 p3 q3
32
' '
Donde las a’ y las q’ son los nuevos elementos que se obtienen de las
operaciones ya mencionadas, y donde p 1 se ha eliminado en la segunda y tercera
ecuaciones. Ahora, multiplicando la segunda ecuación de 3.3 por (-a32‘/a22’) y
sumando el resultado a la tercera ecuación de 3.3, se obtiene el sistema triangular:
a11 a12 a13 p1 q1
0 a' ' '
22 a 23 p2 = q2
0 0 a33 '
a 23
' '
a32 '
/ a 22 p3 q3' q 2' a32
' '
/ a 22
Que simbólicamente puede representarse más simplificadamente así:
a11 a12 a13 p1 q1
0 a' ' '
22 a 23 p 2 = q 2 3.4
0 0 a33 p3 q3
" "
Ejemplo 3.1
Resuelva por eliminación de Gauss el sistema:
4p1 - 9p2 + 2p3 = 500
2p1 - 4p2 + 6p3 = 300
p1 - p2 + 3p3 = 400
SOLUCIÓN
La matriz aumentada del sistema es
4 9 2 500
2 4 6 300
1 1 3 400
Triangularización.- Al sumar la primera ecuación multiplicada por (-2/4) a la
segunda, y la primera ecuación multiplicada por (-1/4) a la tercera, resulta:
4 9 2 500
0 0 .5 5 50
0 1.25 2.5 275
Obsérvese que en este paso la primera fila se conserva sin cambio.
Sumando la segunda fila multiplicada por (-1.25/0.5) a la tercera se obtiene la matriz
4 9 2 500
0 0 .5 5 50
0 0 10 150
Que en términos de sistemas de ecuaciones quedaría como
4p1 - 9p2 + 2p3 = 500
0.5p2 + 5p3 = 50 3.5
-10p3 = 150
Un proceso de sustitución regresiva produce el resultado buscado. La tercera
ecuación de 3.5 da el valor de p3 = -15. De la segunda ecuación se obtiene entonces
0.5P2 = 50 - 5p3 = 125. Y por tanto P2 = 250
Finalmente al sustituir p2 y p3 en la primera ecuación de la forma 3.5 resulta
4p1 = 500 + 9p2 – 2p3 = 2780. De modo que p1 = 695
Con la sustitución de estos valores en el sistema original se verifica la
exactitud de los resultados.
forma 3.3, la fila pivote utilizada fue la segunda. El coeficiente de la incógnita que se
va a eliminar en la fila pivote se llama pivote. En la eliminación que dio como
resultado el sistema de ecuaciones 3.4, los pivotes fueron a11 y a22’. Esta elección
natural de los pivotes a11, a22’, a33”, etc., es muy conveniente tanto para trabajar con
una calculadora como con una computadora, desafortunadamente falla cuando
alguno de esos elementos es cero, puesto que los multiplicadores quedarían
indeterminados (por ejemplo si a11 fuera cero, el multiplicador (a21/a11) no está
definido). Una manera de evitar esta posibilidad es seleccionar como pivote el
coeficiente de máximo valor absoluto en la columna relevante de la matriz reducida.
Como antes, se tomarán las columnas en orden natural de modo que se vayan
eliminando las incógnitas también en orden natural p1, p2, p3, etc. Esta técnica,
llamada pivoteo, se ilustra con la solución del siguiente sistema.
Ejemplo 3.2
Resuelva el sistema
10p1 + p2 - 5p3 = 100
-20p1 + 3p2 + 20p3 = 200 3.6
5p1 + 3p2 + 5p3 = 600
SOLUCIÓN
La matriz aumentada es
10 1 5 100
20 3 20 200
5 3 5 600
El primer pivote debe ser (-20), ya que es el elemento de máximo valor absoluto en la
primera columna. Vamos a la forma triangular en la eliminación. Para esto es
necesario, por ejemplo en la ecuación 3.6, intercambiar la segunda fila (donde se
encuentra el elemento de máximo valor absoluto) con la primera, con lo que se
obtiene
20 3 20 200
10 1 5 100
5 3 5 600
Se elimina entonces P1 de la segunda y tercera filas de la ecuación 3.6. Para
ello, se suma a la segunda fila la primera multiplicada por (-10/(-20)), y a la tercera
fila la primera multiplicada por (-5/(-20)). Con esto se reduce en la primera
eliminación a
20 3 20 200
0 2.5 5 200
0 3.75 10 650
El siguiente pivote debe seleccionarse entre la segunda y tercera filas (segunda
columna) y en este caso es (3.75). se intercambian la segunda y la tercera filas de la
matriz anterior para obtener
20 3 20 200
0 3.75 10 650
0 2.5 5 200
Ejemplo 3.3
Por eliminación de Jordan, resuelva el sistema
4p1 - 9p2 + 2p3 = 500
2p1 - 4p2 + 6p3 = 300
p1 - p2 + 3p3 = 400
SOLUCION
La matriz aumentada del sistema es
4 9 2 500
2 4 6 300
1 1 3 400
Como en la primera columna el elemento de máximo valor absoluto se
encuentra en la primera fila, ningún intercambio es necesario y el primer paso de
eliminación produce
4 9 2 500
0 0 .5 5 50
0 1.25 2.5 275
El elemento de máximo valor absoluto en la parte relevante de la segunda columna
(filas 2 y 3) es 1.25; por tanto, la fila 3 debe intercambiarse con la 2.
4 9 2 500
0 1.25 2.5 275
0 0.5 5 50
Sumando la segunda fila multiplicada por (-(-9)/1.25), a la primera fila y la segunda
multiplicada por (-0.5/1.25) a la tercera, se obtiene el nuevo arreglo
4 0 20 2480
0 1.25 2.5 275
0 0 4 60
Donde se han eliminado los elementos de arriba y abajo del pivote (nótese que en
este paso el primer pivote no se modifica porque sólo hay ceros debajo de él). Por
último, sumando la tercera fila multiplicada por (-20/4) a la primera fila y la tercera
multiplicada por (-2.5/4) a la segunda
4 0 0 2780
0 1.25 0 312.5
0 0 4 60
Que escrita de nuevo como sistema de ecuaciones da
4p1 = 2780
1.25p2 = 312.5
4p3 = -60
De donde el resultado final se obtiene fácilmente
p1 2780 / 4 695
p 2 312.5 / 1.25 250
p3 60 / 4 15
a3,1 a3,1
l 3,1
u1,1 a1,1
a3,1
a 3, 2 a1, 2
a3, 2 l3,1u1, 2 a1,1
l 3, 2 3.10
u 2, 2 a 2,1
a 2, 2 a1, 2
a1,1
a3,1
a 3, 2 a1, 2
u 3,3 a3,3 l3,1u1,3 l3, 2 u 2,3 a3,3
a3,1
a1,3
a1,1 a a 2,1 a
a1,1 a 2,1 2,3 a1,1 1,3
a 2, 2 a1, 2
a1,1
Las ecuaciones 3.8, 3.9 y 3.10, convenientemente generalizadas constituyen un
método directo para la obtención de L y U, con la ventaja sobre la triangularización
de que no se tiene que escribir repetidamente las ecuaciones o arreglos modificados
de Ap = q. A continuación se resuelve un ejemplo.
Ejemplo 3.4
Resuelva por el método de Doolitle el sistema
4p1 -9p2 + 2p3 = 5
2p1 -4p2 + 6p3 = 3
p1 -p2 + 3p3 = 4
SOLUCIÓN
Con l1,1 l 2, 2 l 3,3 1 , se procede al cálculo de la primera fila de U
u1,1 = 4 u1,2 = -9 u1,3 = 2
Cálculo de la primera columna de L
l1,1 = 1 (dato) l2,1 = 2/4 = 0.5 l3,1 = ¼ = 0.25
Cálculo de la segunda fila de U
u2,1 = 0 (recuérdese que U es triangular superior)
u2,2 = -4 –(2/4)(-9) = 0.5
u2,3 = 6 –(2/4)(2) = 5
Cálculo de la segunda columna de L
l1,2 = 0 (ya que L es triangular inferior)
l2,2 = 1 (dato)
l3,2 = (-1-(1/4)(-9))/(-4-(2/4)(-9)) = 2.5
Cálculo de la tercera fila de U, o más bien sus elementos faltantes, ya que por
ser triangular superior
u31 = u32 = 0
u3,3 = 3 -(1/4)(2) -[(-l-(1/4)(-9))/(-4-(2/4)(-9))](6-(2/4)(2)) = - 10
Con esto se finaliza la factorización.
Las matrices L y U quedan como sigue
1 0 0 4 9 2
L = 0 .5 1 0 U = 0 0.5 5
0.25 2.5 1 0 0 10
l i ,i 1 i = i+1, … ,n
0
Con la convención en las sumatorias que 0
k 1
Puede observarse al seguir las ecuaciones 3.8, 3.9, 3.10 o bien las ecuaciones
3.11, que una vez empleada ai,j el cálculo de ui,j o li,j según sea el caso, esta
componente de A no vuelve a emplearse como tal, por lo que las componentes de L
y U generadas pueden guardarse en A y ahorrar memoria de esa manera.
p1k
Si p (k ) p 2k 3.18
p3k
Ejemplo 3.5
Resuelva el siguiente sistema por los métodos de Jacobi y Gauss-Seidel
4 p1 p 2 5
p1 4 p 2 p3 9
3.23
p 2 4 p3 p 4 3
p3 4 p 4 15
SOLUCIÓN
Despejando p1 de la primera ecuación, p2 de la segunda, etc., se obtiene
p2 5
p1
4 4
p1 p3 9
p2
4 4 4 3.24
p2 p4 3
p3
4 4 4
p3 15
p4
4 4
Vector inicial.- Cuando no se tiene una aproximación al vector solución, se
emplea generalmente como vector inicial el vector cero, esto es p ( 0 ) 0 0 0 0
T
a) Método de Jacobi
El cálculo de p (1) en el método de Jacobi se obtiene remplazando p ( 0 ) en cada
una de las ecuaciones de 3.24
0 5 5
p1 1.25
4 4 4
0 0 9 9
p2 2.25
4 4 4 4
0 0 3 3
p3 0.75
4 4 4 4
0 15 15
p4 3.75
4 4 4
T
5 9 3 15
Y entonces p , , ,
(1)
4 4 4 4
Para calcular p se sustituye p (1) en cada una de las ecuaciones de 3.24. Para
( 2)
K p 2k p3k p 4k 1k 2k 3k 4k
k
0 p
0,000000
1 0,000000 0,000000 0,000000
1 1,250000 2,250000 -0,750000 3,750000 1,2500 2,2500 -0,7500 3,7500
2 1,812500 2,375000 0,750000 3,562500 0,5625 0,1250 1,5000 -0,1875
3 1,843750 2,890625 0,734375 3,937500 0,0313 0,5156 -0,0156 0,3750
4 1,972656 2,894531 0,957031 3,933594 0,1289 0,0039 0,2227 -0,0039
5 1,973633 2,982422 0,957031 3,989258 0,0010 0,0879 0,0000 0,0557
6 1,995605 2,982666 0,992920 3,989258 0,0220 0,0002 0,0359 0,0000
7 1,995667 2,997131 0,992981 3,998230 0,0001 0,0145 0,0001 0,0090
8 1,999283 2,997162 0,998840 3,998245 0,0036 0,0000 0,0059 0,0000
9 1,999290 2,999531 0,998852 3,999710 0,0000 0,0024 0,0000 0,0015
10 1,999883 2,999536 0,999810 3,999713 0,0006 0,0000 0,0010 0,0000
b) Método de Gauss-Seidel
Para el cálculo del primer elemento del vector p (1) , se sustituye p ( 0 ) en la primera
ecuación de 3.24, para simplificar la notación se han omitido los superíndices.
0 5 5
p1 1.25
4 4 4
Para el cálculo de p2 de p (1) , se emplea el valor de p1 ya obtenido (1/4) y los
valores de p2, p3 y p4 de p ( 0 ) . Así
1.25 0 9
p2 2.5625
4 4 4
Con los valores de p1 y p2 ya obtenidos y con p3 y p4 de p(0) se evalúa p3 de p (1) .
2.5625 0 3
p3 0.109375
4 4 4
Finalmente, con los valores de p1, p2 y p3 calculados previamente y con p4 de
p(0), se obtiene la última componente de p (1)
0.109375 15
p4 3.722656
4 4
Entonces p (1) [1.25 2.5625 0.109375 3.722656]T
Para la segunda iteración (cálculo de p 2) se procede de igual manera.
2.5625 5
p1 1.890625
4 4
1.890625 0.109375 9
p2 2.6953125
4 4 4
2.6953125 3.722656 3
p3 0.854492
4 4 4
0.854492 15
p4 3.963623
4 4
Con lo que p ( 2 ) [1.890625 2.695313 0.854492 3.963623]T
En la tabla 3.2 se presentan los resultados de las iteraciones subsecuentes.
p 2k p3k p 4k 1k 2k 3k 4k
K p1k
0 0,000000 0,000000 0,000000 0,000000
1 1,250000 2,562500 -0,109375 3,722656 1,2500 2,5625 -0,1094 3,7227
2 1,890625 2,695313 0,854492 3,963623 0,6406 0,1328 0,9639 0,2410
3 1,923828 2,944580 0,977051 3,994263 0,0332 0,2493 0,1226 0,0306
4 1,986145 2,990799 0,996265 3,999066 0,0623 0,0462 0,0192 0,0048
y 3.27
n
aii a ji 1 i n
j 1
j i
k
P2k P3k P4k
K P 1 DIFERENCIAS
0 0,00000 0,00000 0,00000 0,00000
1 -1,50000 1,83333 0,60000 0,16667 -1,50000 1,83333 0,60000 0,16667
2 -2,63333 1,35185 0,59556 0,64815 -1,13333 -0,48148 -0,00444 0,48148
3 -2,14963 1,08807 0,65798 0,91193 0,48370 -0,26379 0,06242 0,26379
4 -1,91705 0,88950 0,71811 1,11050 0,23258 -0,19856 0,06014 0,19856
5 -1,74856 0,72907 0,76865 1,27093 0,16849 -0,16043 0,05053 0,16043
6 -1,61339 0,59784 0,81025 1,40216 0,13516 -0,13124 0,04160 0,13124
7 -1,50296 0,49026 0,84439 1,50974 0,11044 -0,10758 0,03414 0,10758
8 -1,41245 0,40204 0,87239 1,59796 0,09051 -0,08822 0,02800 0,08822
9 -1,33824 0,32970 0,89535 1,67030 0,07422 -0,07234 0,02296 0,07234
10 -1,27737 0,27038 0,91418 1,72962 0,06086 -0,05932 0,01883 0,05932
11 -1,22747 0,22173 0,92962 1,77827 0,04991 -0,04865 0,01544 0,04865
12 -1,18654 0,18183 0,94229 1,81817 0,04093 -0,03990 0,01266 0,03990
13 -1,15297 0,14911 0,95267 1,85089 0,03356 -0,03272 0,01038 0,03272
14 -1,12545 0,12228 0,96119 1,87772 0,02753 -0,02683 0,00852 0,02683
15 -1,10287 0,10028 0,96817 1,89972 0,02257 -0,02200 0,00698 0,02200
Tabla 3.4. Aplicación del método de Gauss-Seidel al sistema 3.26, re-arreglando las ecuaciones para
obtener una aproximación a un sistema diagonal dominante.
ecuación resultante para esa incógnita en todas las demás ecuaciones; con esto el
sistema se reduce en una ecuación y una incógnita. Continúese de esta manera
hasta donde sea posible.
Por ejemplo, en el sistema
f1 ( x1 , x 2 ) 10( x 2 x12 ) 0
f 2 ( x1 , x 2 ) 1 x1 0
se despeja x1 en la segunda ecuación
x1 = 1
y se sustituye en la primera
10(x2 – 12) = 0
cuya solución, x2 = 1, conjuntamente con x1 = 1 proporciona una solución del sistema
dado, sin necesidad de resolver dos ecuaciones con dos incógnitas.
a la solución del Sistema dado. En caso contrario, habría que proponer un nuevo
valor de x3 y repetir el proceso.
Nótese la íntima relación que guarda este método con el método de punto fijo,
ya que un problema multidimensional se reduce a unidimensional en x3:
h(x3) = 0.
a) De consideraciones físicas
Sí el sistema de ecuaciones 3.28 tiene un significado físico, con frecuencia es
posible acotar los valores de las incógnitas a partir de consideraciones físicas. Por
ejemplo, si alguna de las variables xi, representa la velocidad de flujo de un fluido,
ésta no podrá ser negativa. Por tanto xi ≥ 0 En el caso de que xi represente una
concentración expresada como fracción peso o fracción molar de una corriente de
alimentación, se tiene que 0 ≤ xi ≤ 1.
b) De consideraciones geométricas
En caso de tener un sistema de dos ecuaciones con dos incógnitas
f1 ( x1 , x 2 ) 0
f 2 ( x1 , x 2 ) 0
Cada una define, en general, una curva en el plano x1 – x2 y el problema de
resolver el sistema puede verse como el problema de encontrar el punto o los puntos
de intersección de estas dos curvas. Graficando (puede usarse el software GC, el
Math-CAD, o un programa que grafique) pueden obtenerse buenos valores iniciales.
Sea por ejemplo el sistema
f1 ( x1 , x 2 ) 0
f 2 ( x1 , x 2 ) 0
Por último, resulta muy conveniente conocer bien las características de cada método
de solución del sistema 3.28 para efectuar la elección más adecuada del mismo.
repite el proceso, esperando que después de cada iteración los valores de xk, yk se
aproximen a la raíz buscada x, y , la cual cumple con
x g1 ( x , y )
y g 2 (x, y)
Por analogía con los casos discutidos, puede predecirse el comportamiento y
las características de este método de punto fijo multivariable.
Como se sabe, en el caso de una variable la manera particular de pasar de
f(x) = 0 a x = g(x), afecta la convergencia del proceso iterativo. Entonces debe
esperarse que la forma en que se resuelve para x = g1(x, y) y y = g2(x, y) afecte la
convergencia de las iteraciones 3.31.
Por otro lado, se sabe que el reordenamiento de las ecuaciones en el caso
lineal afecta la convergencia, por lo que puede esperarse que la convergencia del
método en estudio dependa de si se despeja x de f2 o de f1.
Finalmente, la convergencia - en caso de existir - es de primer orden, cabe
esperar que el método iterativo multivariable tenga esta propiedad.
Ejemplo 3.6
Encuentre una solución del sistema de ecuaciones no lineales
f1(x, y) = x2 – 10x + y2 + 8 = 0
f2(x, y) = xy2 + x – 10y + 8 = 0.
SOLUCIÓN
Con el despeje de x del término (-10x) en la primera ecuación y de y del término (-
10y) en la segunda ecuación, resulta
x2 y2 8
x
10
xy x 8
2
y
10
o con la notación de la ecuación 3.31
(xk )2 ( y k )2 8
x k 1
10
x (y ) xk 8
k k 2
y k 1
10
Con los valores iniciales x0 = 0, y0 = 0, se inicia el proceso iterativo
Primera iteración
02 02 8
x1 0 .8
10
0( 0) 2 0 8
y
1
0 .8
10
Segunda iteración
(0.8) 2 (0.8) 2 8
x
2
0.928
10
0.8(0.8) 2 0.8 8
y2 0.9312
10
Al continuar el proceso iterativo, se encuentra la siguiente sucesión de vectores
k xk yk
0 0.00000 0.00000
1 0.80000 0.80000
2 0.92800 0.93120
3 0.97283 0.97327
4 0.98937 0.98944
5 0.99578 0.99579
6 0.99832 0.99832
7 0.99933 0.99933
8 0.99973 0.99973
9 0.99989 0.99989
10 0.99996 0.99996
11 0.99998 0.99998
12 0.99999 0.99999
13 1.00000 1.00000
Para observar la convergencia del proceso iterativo, se pudieron usar criterios
como distancia entre dos vectores consecutivos o bien las distancias componente a
componente de dos vectores consecutivos. También existe un criterio de
convergencia equivalente que dice: Una condición suficiente aunque no necesaria,
para asegurar la convergencia es que
g1 g 2 g1 g 2
M 1 ; M 1 3.32
x x y y
para todos los puntos (x, y) de la región del plano que contiene todos los valores
(x , yk) y la raíz buscada ( x , y ) .
k
Ejemplo 3.7
Resuelva el sistema del ejemplo 3.6 utilizando el método de punto fijo multivariable
con desplazamientos sucesivos
f1(x, y) = x2 – 10x + y2 + 8 = 0
f2(x, y) = xy2 + x – 10y + 8 = 0.
Sugerencia: Se pueden seguir los cálculos con un pizarrón electrónico o se
programa una calculadora.
SOLUCIÓN
Al despejar x del término (-10x) y y del término (-10y) de la primera y segunda
ecuaciones, respectivamente, resulta
(xk )2 ( y k )2 8
x k 1 g1 ( x k , y k )
10
k 1 x ( y ) x k 1 8
k 1 k 2
y g2 (x , y ) k k
10
Al derivar parcialmente, se obtiene
g1 2 x k g1 2 y k
x 10 y 10
g 2 ( y ) 1
k 2
g 2 2 x k 1 y k
x 10 y 10
0
y evaluadas en x = 0 y en y = 0 0
g1 g 2 1
0 para x0 y y0
x x 10
g1 g 2
0 0 para x0 y y0
y y
con lo que se puede aplicar la condición 3.32
g1 g 2
+ = 0 + 1/10 = 1/10 < 1
x x
g1 g 2
+ =0+0=0<1
y y
Que se satisface; si los valores sucesivos de la iteración: x1, y1; x2, y2; x3, y3; …la
Primera iteración
02 02 8
x1 0 .8
10
0.8(0) 2 0.8 8
y
1
0.88
10
Cálculo de la distancia entre el vector inicial y el vector [x1, y1]T
x (1) x ( 0 ) (0.8 0.0) 2 (0.88 0.0) 2 1.18929
Segunda iteración
(0.8) 2 (0.88) 2 8
x2 0.94144
10
0.94144(0.88) 2 0.94144 8
y2 0.96704
10
Cálculo de la distancia entre [x2, y2]T y [x1, y1]T
x ( 2 ) x (1) (0.94144 0.8) 2 (0.96704 0.88) 2 0.16608
A continuación se muestran los resultados de las iteraciones
k xk yk |x(k+1)-xk|
0 0.00000 0.00000
1 0.80000 0.88000 1.18929
2 0.94144 0.96705 0.16608
3 0.98215 0.99006 0.04677
4 0.99448 0.99693 0.01411
5 0.99829 0.99905 0.00436
6 0.99947 0.99970 0,00135
7 0.99983 0.99991 0.00042
8 0.99995 0.99997 0.00013
9 0.99998 0.99999 0.00004
10 0.99999 1.00000 0.00001
11 1.00000 1.00000 0.00001
Nótese que se requirieron once iteraciones para llegar al vector solución (1, 1)
contra 13 del ejemplo 3.6, donde se usaron desplazamientos simultáneos.
A continuación se presenta un algoritmo para el método de punto fijo
multivariable en sus versiones de desplazamientos simultáneos y desplazamientos
sucesivos.
1 2 f 2 k 1 2 f 2 k 1 k 1 2 f 2 k 1
(x x ) 2
k 2
( x x )( y y )
k k
( y y k ) 2 ... 3.34
2! xx xy yy
Al igual que en la ecuación 3.33, todas las derivadas parciales de 3.34 están
evaluadas en (xk, yk)
Ahora supóngase que xk+1 y yk+1 están tan cerca de la raíz buscada ( x , y ) que
los lados izquierdos de las dos últimas ecuaciones son casi cero; además, asúmase
que xk y yk están tan próximos de xk+1 y yk+1 que pueden omitirse los términos a partir
de los que se encuentran agrupados en paréntesis rectangulares. Con esto las ecua-
ciones 3.33 y 3.34 se simplifican a
f 1 k 1 f
0 f1 ( x k , y k ) ( x x k ) 1 ( y k 1 y k )
x y
f f
0 f 2 ( x k , y k ) 2 ( x k 1 x k ) 2 ( y k 1 y k ) 3.35
x y
Para simplificar aún más, se cambia la notación con
xk+1 – xk = h
yk+1 – yk = j 3.36
y así queda la (k+1)-ésima iteración en términos de la k-ésimas, como se ve a con-
tinuación
xk+1 = xk + h
yk+1 = yk + j 3.37
La sustitución de la ecuación 3.36 en la 3.35 y el rearreglo dan como resultado
f 1 f
h 1 j f1 ( x k , y k )
x y
f 2 f
h 2 j f2 (xk , y k ) 3.38
x y
El cual es un sistema de ecuaciones lineales en las incógnitas h y j (recuérdese que
las derivadas parciales de la ecuación 3.38, así como f1 y f2 están evaluadas en
(xk,yk) y, por tanto, son números reales).
Este sistema de ecuaciones lineales resultante tiene solución única, siempre
que el determinante de la matriz de coeficientes o matriz jacobiana J no sea cero; es
decir si
f1 f1
x y
J 0
f 2 f 2
x y
Precisando: El método de Newton-Raphson consiste fundamentalmente en
formar y resolver el sistema 3.38, esto último por alguno de los métodos vistos. Con
la solución y la ecuación 3.37 se obtiene la siguiente aproximación.
Este procedimiento se repite hasta satisfacer algún criterio de convergencia
establecido.
Es interesante notar que como en el caso unidimensional, este método puede
obtenerse encontrando un plano tangente a cada f de la ecuación 3.29 en (xk,yk) y
luego encontrar el cero común de estos planos; es decir, hallar un plano tangente en
(xk, yk) tanto a la superficie f1 como a la superficie f2, y luego la intersección de cada
plano tangente con el plano x-y, con lo cual se obtienen dos líneas rectas en el plano
x-y y, por último, la intersección de estas dos líneas rectas, que da el cero común de
los planos tangentes.
Cuando converge este método, lo hace con orden dos, y requiere que el
vector inicial (x0, y0) esté muy cerca de la raíz buscada ( x , y )
Ejemplo 3.8
Use el método de Ncwton-Raphson para encontrar una solución aproximada del
sistema
f1(x, y) = x2 – 10x + y2 + 8 = 0
f2(x, y) = xy2 + x – 10y +8 = 0
con el vector inicial: [x0, y0]T = [0, 0]T
SOLUCIÓN
Primero se forma la matriz coeficiente del sistema 3.38, también conocida como
matriz de derivadas parciales
f1 f1
2 x 10 2y
x y
f 2 f 2
y2 1 2 xy 10
x y
que aumentada en el vector de funciones resulta en
2 x 10 2 y x 2 10 x y 2 8
y 2 1 2 xy 10 xy 2 x 10 y 8
Primera iteración
Al evaluar la matriz en [x0, y0]T se obtiene
10 0 8
1 10 8
que al resolverse por eliminación de Gauss da
h = 0.8, j = 0.88
al sustituir en la ecuación 3.37 se obtiene
x1 = x0 + h = 0 + 0.8 = 0.8
y1 = y0 + j = 0 + 0.88 = 0.88
Cálculo de la distancia entre x(0) y x(1)
x (1) x ( 0 ) (0.8 0) 2 (0.88 0) 2 1.18929
Segunda iteración
Al evaluar la matriz en [x1, y1]T resulta
8 .4 1.76 1.41440
1.7744 8.592 0.61952
que por eliminación gaussiana da como nuevos resultados de h y j
h = 0.36497, j = 0.1117
de donde
x2 = x1 + h = 0.8 + 0.36497 = 1.16497
y2 = y1 + j = 0.88 + 0.1117 = 0.9917
Cálculo de la distancia entre x(1) y x(2)
3.4.3.1 GENERALIZACIÓN
Ejemplo 3.9
Con el algoritmo 3.2, elabore un programa de propósito general para resolver
sistemas de ecuaciones no lineales. Luego úselo para resolver el Sistema
f1(x1, x2, x3) = 3x1 - cos(x2x3) - 0.5 = 0
f2(x1, x2, x3) = x12 - 625x22 = 0
f3(x1, x2, x3) = e-x1x2 + 20x3 + (10π -3)/3 = 0
SOLUCIÓN
La matriz jacobiana ampliada para el sistema es
3 x3 sen( x 2 x3 ) x 2 sen( x 2 x3 ) 3 x1 cos( x 2 x3 ) 0.5
2x 1250 x 2 0 x12 625 x 22
1
x 2 e x1x2 x1e x1x2 x1 x2 10 3
20 e 20 x3
3
Al ejecutar el programa con el vector inicial [1 1 1]T debe producir los siguientes
resultados
k x1 x2 x3 Distancia
0 1.00000 1.00000 1.00000
1 0.90837 0.50065 -0.502861.5863
2 0.49927 0.25046 -0.519040.47982
3 0.49996 0.12603 -0.520450.12444
4 0.49998 0.06460 -0.521990.61446E-01
5 0.49998 0.03540 -0.522720.29214E-01
6 0.49998 0.02335 -0.523020.12052E-01
7 0.49998 0.02024 -0.523090.31095E-02
8 0.49998 0.02000 -0.523100.23879E-03
9 0.49998 0.02000 -0.523100.14280E-05
f 2 y
donde ∂f2/∂y se evalúa en x1, y0. Se tiene ahora x1 y y1. Con estos valores se calcula
x2, después y2, y así sucesivamente.
Este método converge a menudo si x0, y0 está muy cerca de ( x , y ) , y requiere
la evaluación de sólo 2n funciones por paso (cuatro para el caso de dos ecuaciones
que se está manejando).
Nótese que se han empleado desplazamientos sucesivos, pero los
desplazamientos simultáneos también son aplicables.
Ejemplo 3.10
Resuelva el sistema
f1(x, y) = x2 – 10x + y2 + 8 = 0
f2(x, y) = xy2 + x – 10y +8 = 0
con el método Newton-Raphson modificado, usando los valores iniciales x0=0, y0=0.
SOLUCIÓN
Primero se obtiene
f 1 f 2
2 x 10 y 2 xy 10
x y
Primera iteración
Se evalúan f1 y ∂f1/∂x en [0,0]T
f1(0,0) = 8
y
f1 0
x 10
x 0
y
se sustituye
8
x1 0 0 .8
10
Para el cálculo de y1 se necesita evaluar f2 y ∂f2/∂y en x1, y0
F2(0.8, 0) = 0.8(0)2 + 0.8 - 10(0) + 8 = 8.8
f 2 1
x 2(0.8)(0) 10 10
y 0
y
se sustituye
8 .8
y1 0 0.88
10
Segunda iteración
f1 1
f1(0.8, 0.88) = 1.4144 y x 8.4
x 1
y
1.4144
x 2 0 .8 0.96838
8 .4
f 2 2
F2(0.96838, 0.88) = 0.91829 y x 8.29565
y 1
y
de donde
0.91829
y 2 0.88 0.99070
8.29565
Continuar las iteraciones y calcular las distancias entre cada dos vectores
consecutivos. Continuar hasta que xk ≈ 1 y yk ≈ 1. Comparar además la velocidad de
convergencia de este método con la velocidad de convergencia del método de
Newton-Raphson y el de punto fijo para este sistema particular.
En la aplicación de este método se pudo tomar f2 para evaluar x1 y f1 a fin de
evaluar y1, asÍ
f2 (x0 , y 0 )
x x
1 0
f 2 x
f1 ( x1 , y 0 )
y y
1 0
f1 y
Esto puede producir convergencia en alguno de los arreglos y divergencia en
el otro. Es posible saber de antemano si la primera o la segunda forma convergirán
para el caso de sistemas de dos ecuaciones, pero cuando n ≥ 3 las posibilidades son
varias (n!) y es imposible conocer cuál de estos arreglos tiene viabilidad de
convergencia, por lo cual la elección se convierte en un proceso aleatorio. Esta
aleatoriedad es la mayor desventaja de este método.
En general, para un sistema de n ecuaciones en n incógnitas: x1, x2, …, xn, el
algoritmo toma la forma:
k 1 f i ( x1k 1 , x 2k 1 ,..., xik11 , xik ,..., x nk )
x x k
i
f i k 1 k 1
i 1≤i≤n 3.43
( x1 , x 2 ,..., xik11 , xik ,..., x nk )
xi
4.1 INTRODUCCIÓN
a cos(ix) b sen(ix)
n n
4.3
a
0 i i
i 1 i 1
a e
i 0
i
ix
4.4
(d )
i 0
i
2
mínimo
encontrar un polinomio que pase por esos puntos sino más bien que pase entre ellos;
entonces, el método de mínimos cuadrados es aplicable.
Una vez que se obtiene el polinomio de aproximación, éste puede usarse para
obtener puntos adicionales a los existentes en la tabla, mediante su evaluación, lo
que se conoce como Interpolación. También puede derivarse o integrarse a fin de
obtener información adicional de la función tabular.
A continuación se describen distintas formas de aproximar con polinomios ob-
tenidos por ajuste exacto y su uso en la interpolación. Más adelante se describe la
aproximación polinomial por mínimos cuadrados y en los siguientes capítulos la
derivación y la integración.
Fig. 4.1 Aproximación polinomial con criterio de ajuste exacto (curva discontinua) y
con mínimos cuadrados (curva llena).
P(atm) 1 2 5 10 20 30 40
Tabla 4.1 Temperatura de ebullición de la acetona a diferentes presiones
Puntos 0 1 2 3
P(atm) 1 5 20 40
Tabla 4.2 Temperatura de ebullición de la acetona a diferentes presiones
Para obtener los (n + 1) coeficientes del polinomio de grado n (n > 0) que pasa
por (n +1) puntos, proporcionar los
DATOS: El grado del polinomio N y las N + 1 parejas de valores (X(I),
FX(I), I = 0, 1, …,N).
RESULTADOS: Los coeficientes A(0), A(1), …, A(N) polinomio de aproximación.
PASO1. Hacer I = 0
PASO2. Mientras I ≤ N, repetir los pasos 3 a 9.
PASO3. Hacer B(I, 0) = 1
PASO4. Hacer J = 1
PASO5. Mientras J ≤ N, repetir los pasos 6 y 7
PASO6. Hacer B(lJ, J) = B(I, J-1l) * X(I)
PASO7. Hacer J = J ± 1
PASO8. Hacer B(I, N + 1) = FX(I)
PASO9. Hacer I = I + 1
PASO10. Resolver el sistema de ecuaciones lineales Ba = fx de orden N + 1 con
alguno de los algoritmos conocidos
PASO11. IMPRIMIR A(0), A(1), …, A(N) y TERMINAR.
x x
0
0 1
y para hallar el valor de a1, se sustituye el valor de x con el de x1, con lo que
resulta
f ( x1 ) 4.14
a
x x
1
1 0
de tal modo que al sustituir las ecuaciones 4.13 y 4.14 en la 4.12 queda
Ing. Hermas Herrera Callejas Página: 5 de 21
Programación Aplicada Capítulo 4 – Aproximación funcional e interpolación
f (x ) f (x ) 4.15
f ( x) 0
(x x ) 1
(x x )
x x x x
1 0
0 1 1 0
n
4.22
f ( x)
n
L ( x) f ( x )
i i
i 0
donde
n (x x )
Li ( x)
j
4.23
j 0 ( xi x j )
j i
n
(sabiendo que ( x x ) ( x x )( x x
i 1
i 1 2 )...( x x n ) )
Ejemplo 4.1
Para la tabla que se presenta a continuación
a) Obtenga la aproximación polinomial de Lagrange con todos los puntos
b) Interpole el valor de la función f(x) para x = 1.8
I 0 1 2 3
f(xi) -3 0 5 7
Xi 0 1 3 6
SOLUCIÓN
a) Obsérvese que hay cuatro puntos en la tabla, por lo que el polinomio será de
tercer grado. Al sustituir los cuatro puntos en las ecuaciones generales 4.22 y 4.23 se
obtiene
( x 1)( x 3)( x 6) ( x 0)( x 3)( x 6)
f3 ( x) ( 3) ( 0)
(0 1)(0 3)(0 6) (1 0)(1 3)(1 6)
( x 0)( x 1)( x 6) ( x 0)( x 1)( x 3)
(5) (7)
( 3 0)( 3 1)( 3 6) ( 6 0)( 6 1)( 6 3)
al efectuar las operaciones queda
1 5 7
f ( x) ( x 10 x 27 x 18) ( x 7 x 6 x)
3
3 2
( x 4 x 3x)
3 2 3 2
6 18 90
y finalmente resulta
1 3 1 2 46
f 3 ( x) x x x3
30 30 30
b) El valor de x = 1.8 se sustituye en la aproximación pollnomial de Lagrarige de
tercer grado obtenida arriba y se tiene f(1.8) ≈ 2.2176
Obsérvese que si se reemplaza x con cualquiera de los valores dados en la
tabla, en la aproximación polinomial, se obtiene el valor de la función dado por la
misma tabla.
Ejemplo 4.2
Encuentre tanto la aproximación polinomial de Lagrange a la tabla 4.2 como el
valor de la temperatura para una presión de 2 atm utilizando esta aproximación.
SOLUCION
a) Aproximación polinomial de Lagrange mediante dos puntos (n = 1)
p p p p
T f ( p) f (p ) 1
f (p ) 0
p p p p
0 1
0 1 1 0
( p p )( p p ) ( p p )( p p ) ( p p )( p p )
2 0 1 2
0 1 0 2 1 0 1 2 2 0 2 1
( p 5 )( p 20 ) ( p 1 )( p 20 ) ( p 1 )( p 5 )
T f ( p) 56 . 5 113 181
(1 5 )( 1 20 ) ( 5 1 )( 5 20 ) ( 20 1 )( 20 5 )
2
( p p )( p p )( p p ) ( p p )( p p )( p p )
2 3
2 0 2 1 2 3 3 0 3 1 3 2
Ejemplo 4.3
Elabore un programa para aproximar la función f(x) = cos x en el intervalo [0,
8π], con polinomios de Lagrange de grado 1, 2, 3,..., 10. Use los puntos que se
requieran, distribuidos regularmente en el intervalo.
SOLUCIÓN
Para calcular el error máximo dividir el intervalo [0, 8π] en 20 subintervalos y
calcular el valor con el polinomio interpolante y el valor verdadero con la función cos
x, determinando el error absoluto. Se obtienen los siguientes resultados
Grado Error máximo
1 2.23627
2 2.23622
3 3.17025
4 2.23627
5 4.04277
6 4.1879
7 5.68560
8 33.74134
9 12.82475
10 35.95174
Se observa que al aumentar el grado del polinomio, el error absoluto máximo
va aumentando.
Puntos 0 1 2 … n
X x0 x1 x2 … xn
f(x) f(x0) f(x1) f(x2) … f(xn)
f ( x1 ) f ( x0 )
f ' ( ), ( x1 , x0 )
x1 x0
siempre y cuando f(x) satisfaga las condiciones de aplicabilidad de dicho teorema.
Para obtener aproximaciones de derivadas de orden más alto, se extiende el
concepto de diferencias divididas a órdenes más altos como se ve en la tabla 4.3,
donde para uniformar la notación se han escrito los valores funcionales en los
argumentos xi, 0 ≤ i ≤ n, como f[xi] y se les llama diferencias divididas de orden cero.
Información Diferencias Divididas
X f(x) Primeras Segundas Terceras
x0 f[x0]
f [x ] f [x ]
f [x , x ] 1 0
x x
0 1
1 0
f [x , x ] f [x , x ]
f [x , x , x ] 1 2 0 1
x x
0 1 2
x1 f[x1] 2 0
f [x ] f [x ] f [x , x , x ] f [x , x , x ]
f [x , x ] 2 1
f [x , x , x , x ] 1 2 3 0 1
...
2
x x x x
1 2 0 1 2 3
2 1 3 0
f [x , x ] f [x , x ]
f [x , x , x ] 2 3 1 2
x x
1 2 3
x2 f[x2] 3 1
f [x ] f [x ] f [x , x , x ] f [x , x , x ]
f [x , x ] 3 2
f [x , x , x , x ] 2 3 4 1 2
...
3
x x x x
2 3 1 2 3 4
3 2 4 1
f [x , x ] f [x , x ]
f [x , x , x ] 3 4 2 3
x x
2 3 4
x3 f[x3] 4 2
f [x ] f [x ] f [x , x , x ] f [x , x , x ]
f [x , x ] 4 3
f [x , x , x , x ] 3 4 5 2 3
...
4
x x x x
3 4 2 3 4 5
4 3 5 2
f [x , x ] f [x , x ]
f [x , x , x ] 4 5 3 4
x x
3 4 5
x4 f[x4] 5 3
f [x ] f [x ]
f [x , x ] 5 4
x x
4 5
5 4
x5 f[x5]
Tabla 4.3 Tabulación general de diferencias divididas
Ejemplo 4.4
Puntos 0 1 2 3 4 5
X -2 -1 0 2 3 6
f(x) -18 -5 -2 -2 7 142
A partir de ella, elabore una tabla de diferencias divididas.
SOLUCIÓN
Las primeras diferencias divididas mediante los puntos 0,1 y 1,2,
respectivamente, son
5 ( 18) 2 ( 5)
f [ x0 , x1 ] 13 ; f [ x1 , x2 ] 3
1 ( 2 ) 0 ( 1)
La segunda diferencia dividida mediante los puntos (0), (1) y (2) es
3 13
f [ x0 , x1 , x2 ] 5
0 ( 2 )
de igual manera se calculan las demás diferencias divididas, que se resumen
en la siguiente tabla
Puntos x f(x) 1er orden 2do orden 3er orden 4to orden
0 -2 -18
13
1 -1 -5 -5
3 1
2 0 -2 -1 0
0 1
3 2 -2 3 0
9 1
4 3 7 9
45
5 6 142
Debe notarse que todas las diferencias divididas de tercer orden tienen el
mismo valor, independientemente de los argumentos que se usen para su cálculo.
Obsérvese también que las diferencias divididas de cuarto orden son todas cero, lo
cual concuerda con que la tercera y cuarta derivada de un polinomio de tercer grado
son – respectivamente - una constante y cero, sea cual sea el valor del argumento x.
El razonamiento inverso también es válido: si al construir una tabla de diferencias
divididas en alguna de las columnas el valor es constante (y en la siguiente columna
es cero), la información proviene de un polinomio de grado igual al orden de las
diferencias que tengan valores constantes.
Supóngase que se tiene una función dada en forma tabular como se presenta
a continuación
Puntos 0 1 2 3 … n
X x0 x1 x2 x3 … xn
f(x) f[x0] f[x1] f[x2] f[x3] … f[xn]
y que se desea aproximarla preliminarmente con un polinomio de primer grado que
pasa por ejemplo por los puntos (0) y (1). Sea además dicho polinomio de la forma
f(x) = a0 + a1(x - x0), 4.27
donde x0 es la abscisa del punto (0) y a0, a1 son constantes por determinar. Para
encontrar el valor de a0 se hace x = x0 de donde a0 = f(xo) = f[x0] y a fin de encontrar el
valor de a1 se hace x = x1, de donde a1 = (f[x1] - f[x0])/(x1 – x0), o sea la primera
diferencia dividida f[x0, x1].
Al sustituir los valores de estas constantes en la ecuación 4.27 ésta queda
f(x) = f[x0] + (x – x0)f[x0, x1]
o sea un polinomio de primer grado en términos de diferencias divididas. Y si ahora
se desea aproximar la función tabular con un polinomio de segundo grado que pase
por los puntos (0), (1) y (2) y que tenga la forma
f2(x) = a0 + a1 (x - x0) + a2(x - x0)(x –x1) 4.28
donde x0 y x1 vuelven a ser las abscisas de los puntos (0) y (1) y a0, a1 y a2, son
constantes por determinar, se procede como en la forma anterior para encontrar
estas constantes; o sea
Ejemplo 4.5
Elabore una aproximación polinomial de Newton para la información tabular de
las presiones de vapor de la acetona (tabla 4.2) e interpole la temperatura para una
presión de 2 atm.
SOLUCIÓN
Para el cálculo de los coeficientes del polinomio de Newton, se construye la
tabla de diferencias divididas
Diferenciad divididas
Puntos P T=f(p) Primera Segunda Tercera
0 1 56,5
14,12500
1 5 113 -0,50482
4,53333 0,01085
2 20 181 -0,08167
1,67500
3 40 214,5
a) Para n = 1
T = f1(p) = a0 + a1(p-p0) = f[p0] + f[p0, p1](p-p0)
de la tabla se tiene f[p0] = 56.5 y f[p0, p1] = 14.125, de donde
T = f1(p) = 56.5 + 14.l25(p - 1) = 42.375 + 14.125p
ecuación que equivale a las obtenidas anteriormente (4.5 y 4.24).
Si p = 2, f(2) ≈ f1(2) = 56.5 + 14.125(2 - 1) = 70.6 0C
b) Para n = 2
T = f2(p) = a0 + a1(p-p0) + a2(p-p0)(p-p1) = f[p0] + f[p0, p1](p-p0) + f[p0, p1, p2](p-p0)(p-p1)
Ing. Hermas Herrera Callejas Página: 13 de 21
Programación Aplicada Capítulo 4 – Aproximación funcional e interpolación
Ejemplo 4.6
Aproxime la temperatura de ebullición de la acetona a una presión de 30 atm
usando aproximación polinomial de Newton de grado dos (véase Ej. 4.5).
SOLUCIÓN
Se hace pasar el polinomio de Newton por los puntos (1), (2) y (3), con lo que
toma la forma
T = f2(p) = a0 + a1(p – p1) + a2(p – p1)(p – p2)
con los coeficientes dados ahora de la siguiente manera
a0 = f[p1] = 113
a1 = f[p1, p2] = 4.53333
a2 = f[p1, p2, p3] = -0.08167
Al sustituir
T = f2(p) = f[p1] + f[p1, p2](p – p1) + f[p1, p2, p3](p – p1)(p – p2)
= 113 + 4.533(p - p1) – 0.08167(p – p1)(p – p2)
T = f2(p) = 82.16635 + 6.57508p – 0.08167p2
y al evaluar dicho polinomio en p = 30, se obtiene la aproximación buscada
T = f2(30) = 113 + 4.53333(30 - 5) - 0.08167(30 - 5)(30 - 20) = 205.92
El valor reportado en la tabla 4.1 es 205, por lo que la aproximación es buena.
Fig. 4.4 Aproximación polinomial que pasa por entre los puntos
i 1
f1 ( xi ) f ( xi ) d i mínimo
i 1
Para evitar problemas de derivabilidad más adelante, se acostumbra utilizar
las distancias di elevadas al cuadrado:
m m
f 1 ( xi )
i 1
f ( xi ) d i2 mínimo
2
i 1
En la figura 4.5 se observan los puntos tabulados, la aproximación polinomial
f1(x) y las distancias di, entre los puntos correspondientes, cuya suma hay que mini-
mizar.
Si se utiliza
f1(x) = a0 + a1x 4.31
para aproximar la función dada por la tabla, el problema queda como el de
minimizar
m
a a1 xi f ( xi )
2
0 4.32
i 1
Nótese que del número infinito de polinomios que pasan entre los puntos, se
selecciona aquel cuyos coeficientes a0 y a1 minimicen 4.32.
En el cálculo de funciones de una variable, se ha visto que para encontrar el
mínimo o máximo de una función, se deriva y se iguala con cero esa derivada.
Después se resuelve la ecuación resultante para obtener los valores de la variable
que pudieran minimizar o maximizar la función. En el caso en estudio, donde se tiene
una función por minimizar de dos variables (a0 y a1), el procedimiento es derivar
parcialmente con respecto a cada una de las variables e igualar a cero cada
derivada, con lo cual se obtiene un sistema de dos ecuaciones algebraicas en las
incógnitas a0 y a1 o sea,
m
a 0 i 1
(a 0 a1 xi f ( xi )) 2 0
4.33
m
a1 i 1
(a 0 a1 xi f ( xi )) 2 0
Se deriva dentro del signo de sumatoria
m
m
a a x f ( x ) 2
2a 0 a1 xi f ( xi )1 0
i 1 a 0
0 1 i i
i 1
m
m
a a x f ( x ) 2
2a 0 a1 xi f ( xi )xi 0
i 1 a1
0 1 i i
i 1
m m
ma 0 a1 xi f ( xi )
i 1 i 1
m m m
a 0 xi a1 xi2 f ( xi ) xi
i 1 i 1 i 1
El sistema se resuelve por la regla de Cramer y se tiene
m m 2 m m
f ( x i ) xi xi f ( xi ) xi
a 0 i 1 i 1 i 1 i 1
2
4.34
m
m
m xi2 xi
i 1 i 1
m
m
m
m f ( xi ) xi f ( xi ) xi
a1 i 1 i 1
2
i 1
m
m
m xi2 xi
i 1 i 1
que sustituidos en la ecuación 4.31 dan la aproximación polinornial de primer grado
que mejor ajusta la información tabulada. Este polinomio puede usarse a fin de
aproximar valores de la función para argumentos no conocidos en la tabla.
Ejemplo 4.7
En la tabla siguiente se presentan los alargamientos de un resorte corres-
pondientes a fuerzas de diferente magnitud que lo deforman
Puntos 1 2 3 4 5
Fuerza (kgf) x 0 2 3 6 7
Longitud del resorte (m) y 0,120 0,153 0,170 0,225 0,260
Determine por mínimos cuadrados el mejor polinomio de primer grado (recta)
que represente la función dada.
SOLUCION
Para facilitar los cálculos y evitar errores en los mismos, primero se construye
la siguiente tabla
Puntos Fuerza xi Longitud yi x 12 xiyi
1 0 0,120 0 0,000
2 2 0,153 4 0,306
3 3 0,170 9 0,510
4 6 0,225 36 1,350
5 7 0,260 49 1,820
Σ xi = 18 Σ yi = 0.928 Σ xi2 = 98 Σ xiyi = 3.986
Los valores de las sumatorias de la última fila se sustituyen en el sistema de
ecuaciones 4.34 y se obtiene
a0 = 0.11564 y a1 = 0.019434, de donde
f1(x) = 0.11564 + 0.019434x.
El grado del polinomio no tiene relación con el número de puntos usados y
debe seleccionarse de antemano con base en consideraciones teóricas que apoyan
el fenómeno estudiado, el diagrama de dispersión (puntos graficados en el plano x-y)
o ambos.
El hecho de tener la mejor recta que aproxima la información, no significa que
la información esté bien aproximada; quizá convenga aproximarla con una parábola o
una cúbica.
Para encontrar el polinomio de segundo grado f2(x) = a0 + a1x + a2x2 que mejor
aproxime la tabla, se minimiza
a
m
2
0 a1 xi a 2 xi2 f ( xi ) 4.35
i 1
donde los parámetros a0, a1 y a2 se obtienen al resolver el sistema de ecua-
ciones lineales que resulta de derivar parcialmente e igualar a cero la función por
minimizar con respecto a cada uno. Dicho sistema queda
m m m
ma 0 a1 xi a 2 xi2 f ( xi )
i 1 i 1 i 1
m m m m
a 0 xi a1 xi2 a 2 xi3 f ( xi ) xi 4.36
i 1 i 1 i 1 i 1
m m m m
a 0 xi2 a1 xi3 a 2 xi4 f ( xi ) xi2
i 1 i 1 i 1 i 1
cuya solución puede obtenerse por alguno de los métodos vistos anteriormente
Ejemplo 4.8
El calor específico Cp (cal/K gmol) del Mn304 varía con la temperatura de
acuerdo con la siguiente tabla
Puntos 1 2 3 4 5 6
T (oK) 280 650 1000 1200 1500 1700
Cp (cal/K gmol) 32,70 45,40 52,15 53,70 52,90 50,30
Aproxime esta información con un polinomio por el método de mínimos cuadrados.
SOLUCIÓN
El calor específico aumenta con la temperatura hasta el valor tabulado de
1200 oK, para disminuir posteriormente en valores más altos de temperatura. Esto
sugiere utilizar un polinomio con curvatura en vez de una recta, por ejemplo uno de
segundo grado, que es el más simple.
Para facilitar el cálculo de los coeficientes del sistema de ecuaciones 4.36, se
construye la siguiente tabla
Puntos i T xi Cp yi xi2 xi3 xi4 xi y i y i xi2
1 280 32,70 0,78x105 0,022x109 0,062x1011 9156 2,56x106
2 650 45,40 0,42x106 0,275x109 1,785x1011 29510 19,18x106
3 1000 52,15 1,00x106 1,000x109 1,000x1012 52150 52,15x106
4 1200 53,70 1,44x106 1,728x109 2,074x1012 64440 77,33x106
5 1500 52,90 2,25x106 3,375x109 5,063x1012 79350 119,03x106
6 1700 50,30 2,89x106 4,900x109 8,350x1012 85510 145,37x106
Σ Totales 6330 287,15 8,08x106 11,3x109 166,7x1011 320116 415,62x106
Los coeficientes se sustituyen en el sistema de ecuaciones 4.36 y se obtiene:
6a0 + 6330a1 + 8.08x106a2 = 287.15
Ing. Hermas Herrera Callejas Página: 18 de 21
Programación Aplicada Capítulo 4 – Aproximación funcional e interpolación
NOTA
Muchas de las calculadoras de mano cuentan con un programa interno para
obtener esta aproximación; por otro lado, puede usarse un pizarrón electrónico para
los cálculos (sumatorias, solución de ecuaciones, etcétera).
Ejemplo 4.9
Use la aproximación polinomial de segundo grado obtenida en el ejemplo
anterior para aproximar el calor especifico del Mn304 a una temperatura de 800 oK
SOLUCIÓN
Con la sustitución de T = 800 oK en el polinomio de aproximación se tiene
Cp (800) ≈ f2(800) = 22.4066 + 0.0458(800) - 1.694x10-5(800)2 = 48.2 cal/K gmol
Para obtener los N + 1 coeficientes del polinomio óptimo de grado N que pasa
entre M parejas de puntos, proporcionar los
DATOS: El grado del polinomio de aproximación N, el número de parejas de
valores (X(I), FX(l), I = 1, 2, …, M)
PASO 1 Hacer J = 0
PASO 2 Mientras J ≤ (2*N - 1), repetir los pasos 3 a 5.
PASO 3 Si J ≤ N Hacer SS(J) = 0. De otro modo continuar.
PASO 4 Hacer S(J) = 0
PASO 5 Hacer J = J + 1
PASO 6 Hacer I = 1
PASO 7 Mientras I ≤ M, repetir los pasos 8 a 15
PASO 8 Hacer XX = 1
PASO 9 Hacer J = 0
PASO 10 Mientras J ≤ (2*N - 1), repetir los pasos 11 a 14.
PASO 11 Si J ≤ N hacer SS(J) = SS(J) + XX*FX(l)
De otro modo continuar.
PASO 12 Hacer XX = XX*X(l)
PASO 13 Hacer S(J) = S(J) + XX
PASO 14 Hacer J = J + 1
PASO 15 Hacer I = I + 1
PASO 16 Hacer B(0, 0) = M
PASO 17 Hacer I = 0
Ing. Hermas Herrera Callejas Página: 19 de 21
Programación Aplicada Capítulo 4 – Aproximación funcional e interpolación
Con frecuencia se tienen funciones de más de una variable; esto es, f(u, v, z).
Si se sospecha una funcionalidad lineal en las distintas variables; es decir, si se
piensa que la función
y = a0 + a1u + a2v + a3z
puede ajustar los datos de la tabla siguiente
Puntos u v Z y
1 u1 v1 z1 F(u1, v1, z1)
2 u2 v2 z2 F(u2, v2, z2)
3 u3 v3 z3 F(u3, v3, z3)
. . . . .
. . . . .
. . . . .
M um vm zm f(um, vm, zm)
se puede aplicar el método de los mínimos cuadrados para determinar los
coeficientes a0, a1, a2 y a3 que mejor aproximen la función de varias variables
tabuladas. El procedimiento es análogo al descrito anteriormente y consiste en
minimizar la función
m
(a a1u i a 2 vi a3 z i ) y i
2
0
i 1
que derivada parcialmente con respecto de cada coeficiente por determinar:
a0, a1, a2 y a3 e igualada a cero cada una, queda
m m
m m
a1 i 1 i 1
m m
0 1i 2i 3i i
( a a u a v a z ) y 2
2 (a 0 a1u i a 2 vi a3 z i y i )vi 0
a 2 i 1 i 1
m m
a3 i 1 i 1
Ejemplo 4.10
A partir de un estudio experimental acerca de la estabilización de arcilla muy
plástica, se observó que el contenido de agua para moldeo con densidad óptima
dependía linealmente de los porcentajes de cal y puzolana mezclados con la arcilla.
Se tuvieron así los resultados que se dan abajo. Ajuste una ecuación de la forma:
y = a0 + a1u + a2v a los datos de dicha tabla.
Agua(%) Cal(%) Puzolana(%)
y U v
27.5 2.0 18.0
28.0 3.5 16.5
28.8 4.5 10.5
29.1 2.5 2.5
30.0 8.5 9.0
31.0 10.5 4.5
32.0 13.5 1.5
SOLUCIÓN
El sistema lineal por resolver es una modificación del sistema de ecuaciones
4.37 para una función y de dos variables u y v
n a0 + a1Σu + a2Σv = Σy
a0Σu + a1Σu2 + a2Σuv = Σuy
a0Σv + a1Σvu + a2Σv2 = Σvy
Con objeto de facilitar el cálculo del sistema anterior se construye la siguiente tabla:
I ui vi yi ui2 uivi vi2 uiyi viyi
1 2,0 18 27,5 4 36 324 55 495
2 3,5 16,5 28,0 12,25 57,75 272,25 98,00 462,00
3 4,5 10,5 28,8 20,25 47,25 110,25 129,60 302,40
4 2,5 2,5 29,1 6,25 6,25 6,25 72,75 72,75
5 8,5 9,0 30,0 72,25 76,50 81,00 255,00 270,00
6 10,5 4,5 31,0 110,25 47,25 20,25 325,50 139,50
7 13,5 1,5 32,0 182,25 20,25 2,25 432,00 48,00
Σ Totales 45,0 62,5 206,4 407,50 291,25 816,25 1367,85 1789,65
Los coeficientes se sustituyen en el sistema de ecuaciones y al aplicar alguno
de los métodos de solución para sistemas de lineales, se obtiene
a0 = 28.69, a1 = 0.2569, a2 = 0.09607
al sustituir estos valores se tiene
y = 28.69 + 0.2569u + 0.09607v
5.1 INTRODUCCION
En este capitulo se abordan los temas clásicos de integración con procesos finitos de
aproximación.
Una vez que se ha determinado un polinomio pn(x)* de manera que aproxime
satisfactoriamente una función dada f(x) sobre un intervalo de interés, puede espe-
rarse que al integrar pn(x) en forma definida, también aproxime satisfactoriamente la
integral definida correspondiente a f(x).
b
Para estimar I a
f ( x ) dx , los métodos de Newton-Cotes funcionan en
general en dos pasos
f (x ) s f ( x 0 ) dx
b x1
a
f ( x ) dx
x0
0 5.2
0
0
Al integrar se tiene
s2 1 f (x0 )
f (x0 )
s f ( x 0 ) ds h sf ( x 0 )
1
h f ( x0 ) 0 h f ( x0 )
0
2 2
como Δf(x0) = f(x0 + h) – f(x0), Se llega finalmente a:
a f ( x ) dx 2 f ( x 0 ) f ( x 1 )
b h
5.3
El algoritmo del método trapezoidal
Nótese que el lado derecho de la ecuación 5.3 es el área de un trapezoide de
altura h y lados paralelos de longitud f(x0) y f(x1) (véase Fig. 5.2).
Antes de empezar a resolver ejercicios, es conveniente observar que los
métodos vistos y los siguientes sirven también cuando la función f(x) está dada
analíticamente y las técnicas estudiadas en el cálculo integral no dan resultado, o
bien cuando esta función es imposible de integrar analíticamente. En esos casos, la
tabla de puntos se elabora evaluando la función del integrando en abscisas
seleccionadas adecuadamente.
Ejemplo 5.1
5
b) Aproxime A 2 0
( 2 3 x ) dx
4
(1 2 x 3 x 2
c) Aproxime A 3 ) dx
2
/ 2
d) Aproxime A 4 0
sen ( x ) dx
SOLUCION
x2
para realizar la integración x0
p 2 ( x ) dx , se usa la formula de Newton en
diferencias finitas hacia adelante para expresar p2(x)
s ( s 1) 2
p 2 ( x) p 2 ( x 0 sh) f ( x 0 ) sf ( x 0 ) f (x 0 )
2!
al sustituir p2(x) y expresar toda la integral en términos de la nueva variable s, queda
b x2 2
a
f ( x ) dx x0
p 2 ( x ) dx h 0
p 2 ( x 0 sh ) ds
2 2 s ( s 1) 2
h 0
p 2 ( x 0 sh ) ds h 0
( f ( x0 ) s f ( x0 )
2!
f ( x 0 )) ds
2
s2 s3 s2 1
h sf ( x0 ) f ( x0 ) 2 f ( x0 ) 2 f ( x0 ) h 2 f ( x0 ) 2f ( x0 ) 2 f ( x0 )
2 3! 4 0 3
De la definición de la primera y segunda diferencia hacia adelante se tiene
Δf(x0) = f(x0 + h) – f(x0) = f(x1) – f(x0) y
Δ2f(x0) = f(x0 + 2h) – 2f(x0 + h) + f(x0) = f(x2) – 2f(x1) + f(x0)
que sustituidas en la última ecuación dan lugar a
h
f ( x 0 ) 4 f ( x1 ) f ( x 2 )
b
a
f ( x ) dx
3
5.4
que es el algoritmo de Simpson.
Ejemplo 5.2
SOLUCION
0
d) h 2 x0 = 0 x1 = 0+π/4 = π/4 x2 = π/2
2 4
/4
A4 ( sen(0) 4 sen( / 4) sen( / 2)) 1.0023
3
Se recomienda la comparación y discusión de los resultados obtenidos en forma
analítica [casos de los incisos (b), (c) y (d)] y con los obtenidos en el ejemplo 5.1.
n 5 n 4 11n 3 n 2 4
( ) f ( x0 ) ter min os _ fal tan tes ) 5.5
120 16 72 8
A continuación se dan las formulas de Newton-Cotes para integrar cuando n =
1, 2, 3, 4, 5 y 6. El lector puede verificarlas sustituyendo el valor seleccionado de n y
las diferencias correspondientes en términos de sus valores funcionales en la
ecuación 5.5.
n=1
x1 h
x 0 f ( x ) dx 2 ( f ( x 0 ) f ( x 1 )) trapezoidal
n=2
x2 h
x0 f ( x ) dx 3 ( f ( x 0 ) 4 f ( x1 ) f ( x 2 )) Simpson 1/3
n=3
x3 3h
x0
f ( x ) dx
8
( f ( x 0 ) 3 f ( x 1 ) 3 f ( x 2 ) f ( x 3 )) Simpson 3/8
n=4
x4 2h
x0
f ( x ) dx
45
( 7 f ( x 0 ) 32 f ( x1 ) 12 f ( x 2 ) 32 f ( x 3 ) 7 f ( x 4 )) 5.6
n=5
x5 5h
x0
f ( x ) dx
288
(19 f ( x 0 ) 75 f ( x 1 ) 50 f ( x 2 ) 50 f ( x 3 ) 75 f ( x 4 ) 19 f ( x 5 ))
n=6
x6 h
x0
f ( x)dx
140
(41 f ( x0 ) 216 f ( x1 ) 27 f ( x 2 ) 272 f ( x3 ) 27 f ( x 4 ) 216 f ( x5 ) 41 f ( x6 ))
Por ejemplo, en vez de aproximar la integral de f(x) en [a, b] por una recta
(véase Fig. 5.4 a), conviene dividir [a, b] en n sub-intervalos y aproximar cada uno
por un polinomio de primer grado (véase Fig. 5.4 b). Una vez hecho esto, se aplica la
formula trapezoidal a cada sub-intervalo y se obtiene el área de cada trapezoide, de
tal modo que la suma de todas ellas da la aproximación al área bajo la curva f(x).
Esto es
Ejemplo 5.3
Mediante el algoritmo trapezoidal compuesto, aproxime el área bajo la curva
de la siguiente función dada en forma tabular, entre x = - 1 y x = 4.
Puntos 0 1 2 3 4 5
X -1 0 1 2 3 4
F(x) 8 10 10 20 76 238
SOLUCION
Si se toman todos los puntos de la tabla, se puede aplicar cinco veces el
método trapezoidal. Como todos los intervalos son del mismo tamaño (h = 1), se usa
la ecuación 5.8 directamente
1
A (8 2(10 10 20 76) 238) 239
2
Compárese este resultado con la solución analítica (los datos de la tabla
corresponden a la función f(x) = x4 – 2x2 + x + 10).
Ejemplo 5.4
Mediante el algoritmo de Simpson de integración, aproxime el área bajo la
curva del ejemplo 5.3.
SOLUCION
Con los puntos dados de la tabla, se puede aplicar la regla de Simpson en dos
ocasiones; por ejemplo, una vez con los puntos (0), (1) y (2) y otra con los puntos (2),
(3) y (4). Como la integración debe hacerse de x = - 1 a x = 4, se integra entre los
puntos (4) y (5) con el método trapezoidal y la suma será la aproximación buscada:
a) Método de Simpson aplicado dos veces: h1 = h2 = h3 = h4 = 1, entonces puede
usarse la ecuación 5.10
1
A1 (8 4(10 20) 2(10) 76) 74.666
3
b) Método trapezoidal aplicado a los puntos (4) y (5)
1
A2 (76 238) 157
2
por lo tanto, la aproximación al área es
A = 74.666 + 157 = 231.666
Compare este resultado con el obtenido en el ejemplo 5.3 y el resultado de la
solución analítica (la función tabulada es f(x) = x4 – 2x2 + x + 10).
En la figura 5.6 se tiene la curva de la función f(x) que se desea integrar entre
los límites a y b. La parte (a) de la figura muestra cómo se integraría usando un
trapezoide: uniendo el punto A de coordenadas (a, f(a)) con el punto B (b, f(b))
mediante un segmento de recta p1(x). Esto forma un trapezoide de base h = (b-a),
cuya área es:
T = h/2[f(a) + f(b)]
y que podría escribirse como:
T = w1f(a) +w2f(b) 5.11
Donde w1 = w2 = h/2
(De hecho, cualquiera de las fórmulas de integración desarrolladas anteriormente
n
f ( x)dx wi f ( xi ) , donde, por ejemplo, la regla de
b
puede ponerse en la forma a
i 1
Simpson aplicada una vez tendría w1 = w3 = h/3 y w2 = 4h/3. Véase la Ec. 5.4)
El área del trapezoide calculada T, aproxima el área bajo la curva f(x).
Por otro lado, en la aplicación de la cuadratura de Gauss, en lugar de tomar
los dos puntos A y B en los extremos del intervalo, se escogen dos puntos interiores
C y D (véase la parte b de la Fig. 5.6).
Se traza una línea recta por estos dos puntos, se extiende hasta los extremos
del intervalo y se forma el trapezoide sombreado. Parte del trapezoide queda por
encima de la curva y parte por abajo. Si se escogen adecuadamente los puntos C y
D, cabe igualar las dos zonas de modo que el área del trapezoide sea igual al área
bajo la curva y el cálculo del área del trapezoide resultante dé la integral exacta. El
método de Gauss consiste esencialmente en seleccionar los puntos C y D
adecuados. La técnica se deduce a continuación:
Considérese primero, sin que esto implique perder generalidad, que se desea
integrar la función mostrada en la figura 5.7 entre los límites -1 y +1 (si los límites son
distintos, se hace un cambio de variable para pasarlos a -1 y +1). Los puntos C y D
se escogen sobre la curva y se forma el trapezoide con vértices E, F, G, y H.
I 4 w1 z13 w2 z 23 0
De la primera ecuación se tiene que w1 + w2 = 2; nótese también que si
w1 = w2
y
z1 = -z2
Se satisfacen la segunda y la cuarta ecuaciones. Entonces se elige
w1 = w2 = 1
y
z1 = -z2
y al sustituir en la tercera ecuación se obtiene
2
z12 ( z1 ) 2
3
o bien
1
z12
3
de donde
1
z1 0.57735...
3
y queda entonces
z1 = - 0.57735...
z2 = 0.57735...
con lo que se tiene la formula:
1
F ( z )dz w F ( z ) w F ( z
1
1 1 2 2 ) F (0.57735...) F (0.57735...) 5.13
que, salvo porque se tiene que calcular el valor de la función en un valor
irracional de z, es tan simple como la regla trapezoidal; además, trabaja
perfectamente para funciones cúbicas, mientras que la regla trapezoidal la hace solo
para líneas rectas.
Anteriormente se comentó que para integrar en un intervalo distinto de [-1, 1],
se requiere un cambio de variable a fin de pasar del intervalo de integración general
[a, b] a [-1, 1] y así aplicar la ecuación 5.13; por ejemplo, si se desea obtener
5
e
x
dx
0
2
se puede cambiar a z x 1 , de modo que si x = 0, z = -1 y si x = 5, z = 1. El resto
5
de la integral se pone en términos de la nueva variable z y se encuentra que
e x e 5( z 1) / 2 y
5 5
dx d ( ( z 1)) dz
2 2
entonces la integral queda
5 5 1 5( z 1) / 2
0 e dx 2 1e
x
dz
de modo que las condiciones de aplicación del método de Gauss quedan
satisfechas. Al resolver se tiene:
5 1 5( z 1) / 2 5
2 1
e dz ( w1 F (0.57735...) w2 F (0.57735...))
2
5
((1)e 5( 0.577351) / 2 (1)e 5( 0.577351) / 2 ) 0.91752
2
Esto es
5
e
x
dx 0.91752
0
El valor exacto de esta integral es 0.99326.
b
En general, si se desea calcular a
f ( x)dx aplicando la ecuación 5.13, se
cambia el intervalo de integración con la siguiente fórmula
2 x ( a b) ba ab
z de donde x z 5.14
ba 2 2
(Solo es aplicable cuando los límites de integración a y b son finitos)
ya que si x = a, z = -1; y si x = b, z = 1
El integrando f(x)dx en términos de la nueva variable queda
ba ab
f ( x) F ( z )
2 2
y
ba ab ba
dx d ( z ) dz
2 2 2
Por lo que la integral queda finalmente como
b ba 1 ba ab
a f ( x)dx 2 1F ( 2 z 2 )dz
ba ba ab ba ab 5.15
(F ( (0.57735) ) F( (0.57735) ))
2 2 2 2 2
Una cuestión importante es que el método de Gauss puede extenderse a tres
o más puntos; por ejemplo, si se escogen tres puntos no equidistantes en el
segmento de la curva F(z) comprendida entre -1 y 1, se podría pasar una parábola
por los tres como en la regla de Simpson, excepto en que dichos puntos se
escogerían de modo que minimicen o anulen el error. Similarmente es factible elegir
cuatro puntos y una curva cúbica, cinco puntos y una curva cuártica, etc. En general,
el algoritmo tiene la forma:
1
F ( z )dz A w F ( z ) w F ( z
1
1 1 2 2 ) w3 F ( z 3 ) ... wn F ( z n ) 5.16
donde se han calculado los valores de wi y zi por usar y la tabla 5.2 da valores
hasta para seis puntos.
Con dos puntos, el método de Gauss está diseñado para obtener exactitud en
polinomios cúbicos; con tres, se tendrá exactitud en polinomios de cuarto grado y así
sucesivamente.
Los coeficientes y abscisas dadas en la tabla 5.2 sirven para integrar sobre
todo el intervalo de interés, o bien puede dividirse el intervalo en varios sub-intervalos
(como en los métodos compuestos de integración) y aplicar el método de Gauss a
cada uno de ellos.
Número
Coeficientes wi Abscisas zi
de puntos
2 w1 = w2 = 1.0 -z1 = z2 = 0.5773502692
Ejemplo 5.5
x2
1 2
Integre la función e en el intervalo (-0.8, 1.5) por cuadratura de Gauss.
2
SOLUCION
a) Con dos puntos
Cambio de límites de la integral con la ecuación
2 x ( a b ) 2 x 0 .7
z
ba 2 .3
Si x = -08, z = -1; si x = 1.5, z=1
Con el cambio de la función en términos de la nueva variable z queda:
x2
1 1 .5
I
2
0 .8
e 2
dx
1 .5 ( 0 .8 )
1 1
(
1 .5 ( 0 .8 )
z
0 .8 1 .5 2
2 .3 1
dz
) / 2
( 2 .3 z 0 .7 ) 2 / 8
( )e 2 2
e dz
2 1 2 2 2 1
Ejemplo 5.6
2
Halle 0
sen ( x ) dx , por el método de la cuadratura de Gauss utilizando tres puntos.
SOLUCION
Se cambian variable y límites de integración con la expresión
2 x ( a b)
z
ba
como a = 0, b = 2π, entonces
2 x 2 x x
z 1
2
se despeja x: x = π z + π de donde dx = π dz
Se sustituye en la integral
2 1 1
0
sen ( x ) dx 1
sen ( z ) dz
1
sen ( z ) dz
Con el empleo de la ecuación 5.16 con n = 3 y los valores de la tabla 5.2 queda
A ≈ π{0.55555…[sen(π(-0.7745966692) + π)] + 0.88888…[sen(π(0) + π)]
+ 0.55555…[sen(π(0.7745966692) + π)]}
Comparar este resultado con la solución analítica.
La expresión 5.15 puede ponerse en forma más general y adecuada para
programarla así:
b a n (b a ) z i b a
b
a f ( x ) dx
2 i 1
wiF
2
5.17
Para aproximar el área bajo la curva de una función analítica f(x) en el intervalo [a, b],
proporcionar la función a integrar F(X) y los
6.1 INTRODUCCION
dx dx
n
d f ( x) d p ( x)
n
6.1
n
n
n
dx dx
Al diferenciar la fórmula fundamental de Newton dada anteriormente se tiene
n n
d f ( x) d p ( x) d R ( x)
n
n
n
n
n
n
6.2
dx dx dx
d n Rn ( x) n
d n p n ( x)
donde es el error que se comete al aproximar d f (n x) por
dx n n
dx dx
Si las abscisas dadas x0, x1, …, xn están espaciadas regularmente por
intervalos de longitud h, entonces pn(x) puede escribirse en términos de diferencias
divididas. Al sustituir f[x0], f[x0, x1] etcétera en la ecuación de diferencias divididas en
términos de diferencias divididas, se obtiene
f [ x 0 ] 2 f [ x 0 ] 3 f x 0
p n ( x) f [ x 0 ] ( x x 0 ) ( x x 0 )( x x1 ) ( x x 0 )( x x1 )( x x 2 ) ...
h 2! h 2 3! h 3
n f [ x 0 ]
( x x 0 )( x x1 )( x x 2 )...( x x n 1 ) y se tendrá
n! h n
f [ x0 ] 2 f [ x0 ]
d (( x x0 ) ) d (( x x0 )( x x1 ) )
df ( x) dp n ( x) d ( f [ x0 ]) h 2!h 2
dx dx dx dx dx
3 f x0 n f [ x0 ]
d ((x x0 )(x x1 )(x x2 ) ) d ((x x0 )(x x1 )(x x2 )...(x xn1 ) )
3!h 3 ... n!h n 6.3
dx dx
Se desarrollan algunos de los primeros términos y se tiene
f [x0 ] 2 f [ x0 ]
d (( x x 0 ) ) d (( x 2 x 0 x x 1 x x 0 x 1 ) )
dp n ( x ) df [ x 0 ] h 2! h 2
dx dx dx dx
3 f x 0
d (( x ( x 0 x1 x 2 ) x ( x 0 x1 x 0 x 2 x1 x 2 ) x x 0 x1 x 2 )
3 2
)
3! h 3 ...
dx
La primera derivada de f(x) en todo el intervalo [x0, x1] queda aproximada por el
valor constante (f(x1) – f(x0)) / h, el cual es muy diferente del valor verdadero df(x)/dx
en general.
Es importante recordar que hay una estrecha relación entre las diferencias
divididas y las derivadas.
Ejemplo 6.1
Calcule ∂P/∂v cuando v = 2300 cm3 y compárelo con el valor de la derivada analítica.
SOLUCION
Al usar la ecuación 6.6 con los puntos (0), (1) y (2) se obtiene
P 2v v0 v1 2h 2v 4v 2v1 2h 2v v0 v1
P0 0 P1 P2 ; con h = 200
v 2h 2
2h 2
2h 2
2 * 2300 2000 2200 2 * 200 2 * 2000 4 * 2300 2 * 2200 2 * 200
2
13.782 12.577
2 * 200 2 * 200 2
2 * 2300 2000 2200
11.565 0.00506
2 * 200 2
La derivada analítica es
P RT 2a 82.1 * 350 2(3.6 x10 6 )
0.005048
v (v b) 2 v 3 (2300 42.8) 2 2300 3
Nótese que la aproximación es muy buena (error relativo = - 0.24%) a pesar
de haber aplicado un polinomio de segundo grado para aproximar la ecuación de
Van der Waals que, como se sabe, es un polinomio de tercer grado en v.
Ejemplo 6.2
SOLUCION
n n x xj
De la ecuación p n ( x) f ( xi ) se deriva con respecto a x
i 0 j 0 xi x j
j i
p n ( x) n
d n x xj
f ( xi )
x dx j 0 xi x j
i 0
j i
Se hace
n x xj
y y se toman logaritmos en ambos lados, con lo que se tiene
j 0 xi x j
j i
n x xj n x xj
ln y = ln x ln
j 0 i xj j 0 xi x j
j i j i
Se despeja dy / dx
n
dy 1
y
dx j 0 x x j
j i
y finalmente
dp n ( x) n
n x xj n 1
f ( xi )
dx j 0 xi x j j 0 x x j
i 0
j i
j i
Obsérvese que esta ecuación no sirve para evaluar la derivada en una de las
abscisas de la tabla, ya que significaría dividir entre cero en la sumatoria dentro del
paréntesis. Sin embargo, manipulado algebraicamente el lado derecho puede
escribirse en la forma
n n n
dp n ( x) f ( xi )
n ( x x j ) 6.16
dx i 0
k 0 j 0
( x i x j ) k i j k ,i
j 0
j i
La cual ya no tiene la limitante mencionada.
Ejemplo 6.3
T(K)
P (kg/cm2) (T0) 273 (T1) 300 (T2) 325 (T3) 360
1 0.99 0.97 0.96 0.93
2 0.88 0.82 0.79 0.77
8 0.62 0.51 0.48 0.45
15 0.56 0.49 0.46 0.42
20 0.52 0.44 0.41 0.37
SOLUCION
C A
Lo que se busca es en si , que se puede evaluar con la ecuación 6.16.
T T 300 , P 8
2 2 2
dp2 ( x ) f ( xi )
Al desarrollarla para n = 2 se tiene 2 ( x x j
)
i 0
dx k 0 j 0
( xi x j ) k i j k ,i
j 0
j i
dp2 ( x ) (( x x2 ) ( x x1 )) f ( x0 ) (( x x2 ) ( x x0 )) f ( x1 ) (( x x1 ) ( x x0 )) f ( x2 )
dx ( x0 x1 )( x0 x2 ) ( x1 x0 )( x1 x2 ) ( x2 x0 )( x2 x1 )
dp 2 ( x) (2 x x1 x 2 ) f ( x0 ) (2 x x0 x 2 ) f ( x1 ) (2 x x0 x1 ) f ( x 2 )
dx ( x0 x1 )( x0 x 2 ) ( x1 x0 )( x1 x 2 ) ( x 2 x0 )( x 2 x1 )
donde f(x) representa a CA y x a T; de tal modo que sustituyendo los tres puntos
enmarcados de la tabla queda
dp 2 ( x) C A (2 * 300 300 325)0.62 (2 * 300 273 325)0.51
dx T T 300 (273 300)(273 325) (300 273)(300 325)
P 8
Ejemplo 6.4
Obtenga la primera y segunda derivadas evaluadas en x = 1 para la siguiente
función tabulada
Puntos 0 1 2 3 4
x -1 0 2 5 10
f(x) 11 3 23 143 583
SOLUCION
Al construir la labia de diferencias divididas se tiene
Tabla 6.3 Diferencias divididas de la función
Diferencias divididas
Puntos X f(x) Primeras Segundas
0 -1 11
-8
1 0 3 6
10
2 2 23 6
40
3 5 143 6
88
4 10 583
p2(x) = f[x0] + (x – x0)f[x0, x1] + (x2 – x0x - x1x + x0x1)f[x0, x1, x2]
que al derivarse da
f x0 , x1 (2 x x0 x1 ) f x0 , x1 , x 2
dp 2 ( x)
dx
y al derivarlo nuevamente se obtiene
d 2 p 2 ( x)
2 f x0 , x1 , x 2
dx 2
con la sustitución de valores finalmente resulta
dp 2 (1)
8 (2 * 1 (1) 0) * 6 10
dx
y
d 2 p 2 (1)
12
dx 2
ALGORITMO 6.1 Derivación con polinomios de Lagrange
Para obtener una aproximación a la primera derivada de una función tabular
f(x) en un punto x, proporcionar los
PASO 1. Hacer DP = 0
PASO 2. Hacer I = 0
PASO 3. Mientras I ≤ N, repetir los pasos 4 a 21.
PASO 4. Hacer P = 1
PASO 6. Hacer I = 0
PASO 6. Mientras J ≤ N, repetir los pasos 7 a 8.
PASO 7. Si I <> J Hacer P = P*(X(I) - X(J))
PASO 8. Hacer J = J + 1
PASO 9. Hacer S = 0
PASO 10. Hacer K = 0
PASO 11. Mientras K ≤ N, repetir los pasos 12 a 19.
PASO 12. SI I < > K, realizar los pasos 13 a 18.
PASO 13. Hacer P1 = 1
PASO 14. Hacer J = 0
PASO 16. Mientras J ≤ N, repetir los pasos 16 a 17.
PASO 16. SI J<>I y J<>K Hacer P1 =
P1*(XD - X(J))
PASO 17. Hacer J = J + 1
PASO 18. Hacer S = S + P1
PASO 19. Hacer K = K + 1
PASO 20. Hacer DP = DP + FX(I) / P * S
PASO 21. Hacer l = I + 1
PASO 22. IMPRIMIR DP y TERMINAR.
Lee Gonzalez y Eakin (1966) presentaron una relación semi empírica para
calcular la viscosidad del gas natural. Los autores expresaron la viscosidad del
gas en función a la temperatura del reservorio, densidad del gas, y el peso
molecular del gas. La expresión general de esta correlación esta dada por:
4 Y
g 10 K exp( X g ) ..................................................................................... (1)
1 .5
(9.379 0.01607 M w )T
K ............................................................................. (2)
209.2 19.26 M w T
X 3.448
986.4 0.01009M .........................................................................(3)
T w
donde
donde
, son constantes
es la temperatura absoluta.
En forma logarítmica, esta ecuación es
Barr
donde
es la viscosidad cinemática en
Ecuación de la viscosidad ASTM
fue deducida por Walther:
Los petróleos crudos que se extraen de los diferentes campos petrolíferos son
de naturalezamuy variada, incluso en su apariencia externa. Así por ejemplo
existen petróleos de color amarillento, de gran volatilidad y fluidez, otros de
color negro de menor fluidez; otros de color negro-castaño oscuro, viscosos y
de extrema dificultad para fluir, algunos otros que incluso solidifican a
temperatura ambiente, dando lugar a una masa de consistencia semi-sólida
etc.
A pesar de estas diferencias externas, en algunos casos muy pronunciadas,
considerados químicamente se asemejan grandemente unos a otros ya que
son fundamentalmente mezclas de hidrocarburos, es decir combinaciones de
los elementos químicos Carbono (C) e Hidrógeno (H). De estas combinaciones
surge una enorme variedad de posibilidades y de formación de compuestos
análogos, denominados “familias” de hidrocarburos, que se van formando
según la cantidad de átomos de carbonos que formen la molécula.
Las propiedades físicas, tales como punto de ebullición, peso específico, punto
de escurrimiento, viscosidad, etc. varían siempre proporcionalmente a medida
que aumenta el número de átomos de carbono, por lo que en la medida que las
“cadenas” de carbono se hacen más numerosas y complejas, esos valores
aumentan.
La viscosidad es una propiedad física que se puede definir como la resistencia
interna que un fluido opone al desplazamiento de sus partículas, resistencia
que se debe al frotamiento de las moléculas cuando se deslizan unas contra
otras. En forma práctica, la viscosidad es la resistencia que ofrece un líquido
para moverse.
La viscosidad es una propiedad característica de un determinado fluido, es un
criterio particularmente importante para definir el diseño de los conductos. No
sólo se trata de resistencias internas del fluido sino que este efecto físico
también interactúa con las paredes del caño, de tal manera que un fluido muy
viscoso necesitará un mayor esfuerzo para moverse que uno de baja
viscosidad y por consiguiente el efecto de fricción con las paredes del tubo será
mucho mayor en el de mayor viscosidad. La consecuencia final de esta
situación es que frente a una mayor viscosidad es mayor la fricción y mayor el
gradiente de presión necesario para movilizar dicho fluido dentro de esa
cañería.
La viscosidad de un fluido depende fundamentalmente de la temperatura y de
la cantidad de gas disuelto que contenga. Al contener mayor cantidad de gas
disuelto o al suministrarle calor a un fluido, el efecto que se obtiene es
aumentar la energía interna de sus moléculas, por lo que aumenta la actividad
molecular y las mismas se alejan unas de otras dando por resultado una menor
fricción y por lo tanto una menor resistencia interna al movimiento, o sea una
menor viscosidad.
Por esta razón, en general y sobre todo en los líquidos, la viscosidad disminuye
cuando aumenta la temperatura del fluido, y aumenta cuando su temperatura
disminuye.
En el fondo del pozo el petróleo está a mayor temperatura y a mayor presión
que en superficie, por lo que es mayor también la cantidad de gas disuelto que
contiene y como consecuencia su viscosidad será menor que en superficie.
Cuando sube a superficie, la temperatura del petróleo disminuye y la presión a
la que está sometido también, por lo que el gas disuelto contenido es cada vez
menor, ya que éste se va liberando. Esta condición (menor temperatura y
menor gas disuelto), hace que aumente gradualmente su viscosidad a medida
que llega a superficie.
Una vez en la cañería de conducción y cuando es necesario desplazar un
petróleo viscoso a baja temperatura, es cuando se requiere una mayor
diferencia de presión.
En lugares de baja temperatura ambiental será conveniente (y a veces
imprescindible) disminuir la viscosidad de alguna manera, lo que generalmente
se logra calentando el fluido.
Para obtener los valores de viscosidad se hace necesario realizar ensayos en
laboratorio utilizando viscosímetros de alta calidad. Para expresar estos
resultados se utilizan diferentes unidades, según sea el sistema de que se
trate. En el CGS, la unidad para expresar la viscosidad absoluta (que es la del
propio interior del líquido) es el POISE (o el centipoise), y para la viscosidad
cinemática se utiliza el STOKES (o el centistokes) que se obtiene dividiendo la
viscosidad dinámica por la densidad.
Junto a esta medida absoluta de viscosidad existen en los diferentes países
varias escalas empíricas de la misma, como por ejemplo los sistemas Engler,
Redwood y Saybolt. Los viscosímetros correspondientes determinan el tiempo
que tarda en fluir una cantidad conocida de líquido a través de un orificio
calibrado. Así por ejemplo en la escala Saybolt existen dos tipos distintos de
orificios calibrados. Uno estándar, en el que el tiempo se expresa en Segundos
Saybolt Universales (0.176 cm) y otro de mayor diámetro, en el que el tiempo
se expresa en Segundos Saybolt Furol (0.315 cm), para los líquidos muy
viscosos.
Como la viscosidad es dependiente de la temperatura, es necesario siempre
relacionar ambos valores, la viscosidad y la temperatura a la que se tomó. Es
también conveniente contar con dos o tres lecturas a distintas temperaturas
para observar el comportamiento y la variación que se produce.
Un error medio de 0.64% con una desviación normal de 13.53% era informado
para la correlación cuando probó contra los datos usados para su el desarrollo.
Sutton y Farshad (1980) informó un error de 114.3% cuando la correlación se
probó contra 93 casos de la literatura.
7.2.3.3 La Correlación de Glasso
Glaso (1980) propuso una relación matemática generalizada para computando
la viscosidad de petróleo muerto. La relación se desarrolló de las dimensiones
experimentales en 26 muestras de petróleo crudo. La correlación tiene la
forma siguiente:
5. DESARROLLO
El desarrollo del presente trabajo se basa en el cálculo del factor de
desviación de los gases reales del factor Z por el método de Beggs y Brill.
El cálculo del factor Z se los realiza con la siguiente ecuación:
1 A
z A B C * Ppr
D
(24)
e
El cálculo de los parámetros de la anterior ecuación se lo realiza con las
siguientes ecuaciones:
A 1,39 * (Tpr 0,92) 0,5 0,36 * Tpr 0,101 (25)
0,066 0,32
B (0,62 0,23 * T pr ) * Ppr ( 0,037) * Ppr2 9(Tpr 1) * Ppr6 (26)
T pr 0,86 10
C 0,132 0,32 * log T pr (27)
( 0 , 3016 0 , 49T pr 0 ,1824 T pr
2
D 10
)
(28)
Los datos de Tpr y Ppr pueden ser introducidos directamente o ser
calculados en el mismo programa por medio de una base de datos:
Los gráficos de Obert dan errores menores del 6% salvo en el punto crítico.
La temperatura de Boyle es aquella para la cual:
(13)
donde,
(14)
(15)
(16)
(17)
(18)
F = A11.(0,27.sP r/sT r) ² (19)
Para encontrar el valor de Z que sea solución de la ecuación (13) se aplica el
método de Newton - Raphson que involucra los siguientes pasos:
Se calcula sT r y sP r
Se calculan las constantes A - F.
Se escribe la ecuación (13) como:
(20)
Se supone un valor de Z(Z0), se recomienda mayor que 1, y se chequea si
hace F(Z0) = 0 dentro de la tolerancia requerida.
Si F(Z0) = 0, el valor supuesto es el correcto y es el valor que se está
buscando.
si F(Z0) ≠ 0 se busca un nuevo valor de Z(Z1) de la siguiente manera:
Z1 = Z0 - F(Z0)/F'(Z0) (21)
Donde F'(Z0) es la derivada de F(Z), dada por:
(22)
y calculada en Z0.
Con Z1 se chequea si F(Z1) = 0 y si
no lo es se calcula un valor Z2
usando la ecuación (22) cambiando
Z1 por Z2 y Z0 por Z1.
El procedimiento continua hasta
encontrar un valor Zn que haga F(Z)
= cero; después del primer valor
supuesto para Z(Z0) los demás
valores usados se obtienen a partir
de la ecuación (21) usando Zn en
lugar de Z1 y Zn-1 en lugar de Z0.
COMPORTAMIENTO REAL
Se presenta a altas P donde la ec de estado antes presentada ya no es aceptada por los errores
que presenta. Esta desviación del comportamiento ideal aumenta con la T y la complejidad de
la mezcla de gases, esta desviación se da por las consideraciones planteadas en la teoría
cinética de los gases.
Para describir este comportamiento real, distintas ecuaciones de estado fueron desarrolladas,
sin embargo el caso que se aplica en mayor proporción en la industria petrolera es el de la
ecuación de estado corregida con la inserción del factor de desviación z, esto es:
PV = znRT
FACTOR “Z”
Este factor es adimensional y generalmente varía entre 0,7 y 1,2, donde para el caso en que z=
1, representa el comportamiento ideal, matemáticamente se define por:
Para el caso del gas natural considerando distintas composiciones, se observo que este factor
puede ser generalizado con suficiente exactitud para propósitos de ingeniería cuando es
expresado en función de términos pseudo reducidos de P y T, los cuales a su vez se calculan con
ayuda de las propiedades críticas de los compuestos y las pseudo-crítas de la mezcla.
PROPIEDADES CRÍTICAS
Es el conjunto de condiciones físicas de presión, temperatura y volumen, a las cuales la
densidad y otras propiedades del líquido y gas se vuelven idénticas, es decir, es un punto a una
presión y temperatura dada donde físicamente no puede diferenciarse si se trata de gas o
líquido. Estas propiedades críticas son únicas (una sola presión, una sola temperatura) para
una sustancia dada o mezcla, en cuyo caso son propiedades pseudo-criticas. Estas se requieren
para la determinación de otras propiedades del sistema en análisis ya sea sustancia pura o
mezcla gaseosa. Como en la determinación de los términos reducidos
La presión crítica, Pcr, y la temperatura crítica, Tcr, son medidas en el laboratorio y usualmente
son desconocidas por lo que se requiere para su determinación el empleo de correlaciones,
como:
TPR Y PPR
La presión y temperaturas reducidas son términos adimensionales que en el caso de mezclas
gaseosas, están definidos por:
Donde:
Estas propiedades pseudo críticas no son exactamente las propiedades de la mezcla, sin
embargo para efectos de estimación de otras propiedades se las emplea como propiedades
generadoras aceptables.
Con estos datos es posible determinar el valor de z, ya sea mediante
La grafica presentada por Standing-katz
Métodos directos de cálculo de Z, como:
o Hall-Yarborough
o Dranchuk-Abu-Kassem
o Dranchuk-Purvis-Robinson
Hall-Yarborough
Se caracteriza por:
Método basado en la ecuación de estado de Starling-Carnahan.
Coeficientes obtenidos del ajuste realizado a los datos de la grafica de Standing-Katz.
No recomendado para rangos de Tpr < 1
La expresión matemática presentada es:
Donde:
Ppr Presión pseudo-reducida
t Tpc/T
Y densidad reducida, obtenida de:
Donde:
Dranchuk-Abu-Kassem
Se caracteriza por:
Método basado en la densidad critica, definida por:
Coeficientes A1 a A11 fueron obtenidos del ajuste realizado a los datos de la grafica de
Standing-Katz:
Dranchuk- Purvis-Robinson
Se caracteriza por:
Método basado en la ecuación de estado de Benedict-Webb-Rubin:
Donde la densidad reducida se define por:
Donde:
(Tc)C7+ : Temperatura critica del C7+
(Pc)C7+: Presión critica del C7+
yC7+ : Fracción molar del C7+
(3) Ajustar los parámetros J y K con EJ y EK.
COMP SIMBOLO M PC TC
METANO C1 16,042 673,1 343,2
ETANO C2 30,068 708,3 549,9
PROPANO C3 44,094 617,0 666,0
BUTANO C4 58,120 529,0 734,5
ISO-BUTANO i-C4 58,120 550,0 765,7
PENTANO C5 72,146 483,5 829,6
ISO-PENTANO i-C5 72,146 489,8 846,2
HEXANO C6 86,720 440,1 914,2
HEPTANO+ C7+
Luego se calcula z de forma directa con la ayuda de los métodos analíticos presentados en el
marco teórico, considerando sus rangos de aplicación para Ppr y Tpr, lo cual determina el
método a utilizar en el programa.
El cálculo de z se lo realiza por sub rutinas incrustadas en un modulo que contiene el código
correspondiente para cada método.
Los resultados finales se presentan tanto para z corregida como para z sin corrección por C7+,
en cuadros que resaltan esta diferencia.
Donde:
Pci = Presiones críticas de cada componente [psia]
Yi = Composición de cada componente
DESVENTAJAS DE LA GRÁFICA DE Z
PASOS A SEGUIR:
T PC
i 1
y i * Tc i
TPC B (1 B )
PC
Donde:
Tpc = temperatura pseudocrítica [R]
Ppc = Presión pseudos crítica, [psia]
T’pc= temperatura pseudocrítica corregida [R]
P’pc= presión pseudocrítica corregida [psia]
P-4: Cálculo de las propiedades pseudo-reducidas:
presión del sistema psia
PPr
presión pseudocrítica psia
temperatura del sistema R
TPr
temperatura pseudocríticaR
P-5: Leer el valor del factor z en gráfica
Standing-Katz.
…………(2)
Donde:
Los autores propusieron una ecuación de estado de once constantes para calcular la densidad reducida del
gas
DRANCHUK-PURVIS-ROBINSON:
Desarrollaron un correlación basada en un tipo de
ecuación de estado de Benedict-Webb-Rubin.
Haciendo coincidir la ecuación con 1500 puntos de
la gráfica de z se optimizaron 8 coeficientes de
la ecuación propuesta:
……………………(5)
Factores geológicos.-
Factores físicos.-
Comportamiento de flujo.-
Relación que puede ser resuelta par los diferentes valores de Qi en [bbl/Dìa].
Zona intermedia.-
+
De manera analítica la ecuación anterior puede expresarse como
Parcial
Limitado
Condición de borde externo
Infinito El efecto de la declinación de presión no se siente en el borde externo.
La presión en el borde externo es igual a pi
Regímenes de flujo
Existen tres regimenes de flujo que influencian la tasa de influjo de agua hacia
el yacimiento:
Geometrías de flujo
Empuje lateral
Empuje lineal
Empuje de fondo
Para un Gas Ideal, el factor de compresibilidad es unitario, mientras que para Gases Reales es
mayor o menor que 1.
Ejemplos para el H2O, CO2 y O2 gaseosos:
Ecuación de Van der Waals; Es la ecuación de estado “por excelencia” de los Gases Reales.
Van der Waals atribuyó las desviaciones de los gases de la idealidad debido a:
- El volumen de las moléculas sí importa, no es despreciable
- Las fuerzas de interacción entre moléculas de los gases influye
Efecto del Volumen de las Partículas
b = covolumen (volumen efectivo ocupado por 1 mol de gas)
V = volumen total (ocupado por el gas)
Vdisponible = (Vreal – nb) nb = volumen ocupado por “n” moles de gas
Reemplazando en la Ley Ideal:
P = nRT/(V – nb)
Efecto de las Fuerzas de Interacción
Cuando las moléculas interactúan, la presión total disminuye en un factor proporcional a la
densidad de moléculas.
a = parámetro de interacción, que indica cuan fuertes son las atracciones
P = nRT/(V – nb) – an2 /V2
Con lo que se llega a la Ecuación de Van der Waals, para Gases Reales con desviaciones
moderadas de la Idealidad:
Donde a y b son las “constantes de Van der Waals”, conocidas para los distintos gases.
Unidades de los parámetros de Van der Waals: a [atm l2 /mol2]; b [l/mol]
Ejemplo
Nota: los valores grandes de “a” indican gran interacción entre las moléculas
A parte de la ec. de Van der Waals, existen una serie de ecuaciones de estado que definen el
comportamiento de los Gases Reales para determinadas condiciones:
ECUACIÓN DE REDLICH-KWONG
Difiere de la ec. de Van de Waals al expresar el potencial de atracción (o de interacción) como
una función más complicada de la temperatura y el volumen molar:
[P + an2 / (T1/2V(V+b))] (V – nb) = nRT
Si llenamos una olla a presión con agua pura y calentamos lentamente, pasaremos de agua
líquida (1 fase) a la “campana de ebullición”, donde coexistirán agua líquida y vapor de agua en
equilibrio (2 fases). Si la temperatura y presión siguen aumentando, llegaremos a un punto
superior de la campana denominado punto crítico, sobre el cual el líquido y el gas pierden sus
límites de fase y se transforman en una sola fase, denominada fluido supercrítico. Este estado
de la materia es muy importante en Geología. En espacio PVT, diagrama de fases tiene el
siguiente aspecto:
El diagrama PVT en 3 dimensiones nos permite observar el estado crítico en los espacios PV y
PT:
Por fortuna, las ecuaciones termodinámicas de estado para gases reales también pueden
utilizarse para los fluidos supercríticos. Cada especie de gas puro (agua, dióxido de carbono,
oxígeno, etc.) tiene un único punto crítico, sobre el cual desaparecen los límites entre el líquido
y el vapor. Para el punto crítico se tiene que:
(∂p/∂V)Tc = 0 y (∂2p/∂V2)Tc = 0
Si aplicamos estas condiciones a la Ecuación de Van der Waals, podremos expresar las
constantes a y b en función de la presión y temperatura críticas (variables conocidas y tabuladas
para cada especie):
aVdW = 0.4219 (R2T2c)/Pc bVdW = 0.1250 (RTc)/Pc
Y para el caso de la ECUACIÓN DE REDLICH-KWONG, también es válido expresar los
parámetros a y b en función de las propiedades críticas:
aR-K = 0.4275 (R2T2.5c)/Pc bR-K = 0.0866 (RTc)/Pc
ECUACIÓN DE REDLICH-KWONG
Por tanto el factor acéntrico se define como la diferencia logarítmica evaluada a la temperatura,
Tr = 0,7.
Por definición, el valor del factor acéntrico w es cero para el argón, criptón y xenón; los datos
experimentales que permiten calcular el factor de compresibilidad Z para estos tres fluidos se
correlacionan por la misma curva cuando Z se representa como función de Tr y pr. Por ésto, el
teorema de estados correspondientes con tres parámetros indica que todos los fluidos que tienen
el mismo valor de w, poseen
el mismo valor de Z cuando se comparan a las mismas Tr y pr.
DESARROLLO Y PLANTEAMIENTO
Que tiene tres raíces para el volumen, de las que dos pueden ser complejas; físicamente, los
valores significativos de v siempre son reales, positivos y mayores que la constante b.
En la siguiente figura se observa que:
- Cuando, T > Tc, cualquier valor positivo de p conduce a una solución con una sola raíz real
positiva.
- Cuando, T = Tc, lo antes dicho sigue siendo verdad, excepto a la presión crítica donde existe
una raíz triple,vc.
- Para, T < Tc, a presiones elevadas solo existe una raíz real positiva, pero en el intervalo de
presiones bajas se tienen tres raíces reales positivas (puntos A, D, B; en este caso, la raíz
intermedia no tiene significado físico, la raíz menor es un volumen líquido o casi-líquido y la
raíz mayor es un volumen de vapor o casi-vapor.
Los volúmenes de líquidos y vapores saturados están dados por la raíz menor y mayor,
respectivamente, cuando p es la presión de saturación o presión de vapor. Aunque las raíces de
una ecuación cúbica de estado se pueden encontrar explícitamente, es más frecuente que se
empleen técnicas iterativas, que resultan prácticas solamente si convergen en la raíz deseada.
No se puede dar una seguridad absoluta a este respecto, pero con frecuencia, las consideraciones
que presentamos a continuación resultan efectivas en la ecuación de Redlich/Kwong.
y así sucesivamente se continúa hasta que se llega a un valor Z con un error menor que un cierto
valor preestablecido.
ρ = P.M/Z.RT. = m/V
donde,
P pr: presión pseudoreducida de la mezcla.
T pr: temperatura pseudoreducida de la mezcla.
En el caso de que no se tenga a disposición los datos de composición del gas natural, los
parámetros psudoreducidos pueden ser predecidos a través del valor de la gravedad específica
y las ecuaciones de Brown (1948), presentados en forma gráfica con una aproximación
conveniente de las propiedades pseudocríticas de los gases cuando solo tenemos este dato
disponible. Las expresiones matemáticas que pueden emplearse para efectuar el cálculo son
las siguientes:
Donde:
Tpc =Temperatura pseudo crítica, °R
Ppc = Presión Pseudo crítica, Psia
g = Gravedad específica de mezcla de gas
Limitaciones:
Máximo.
5% N2
2% CO2
2% H2S
Donde:
Tpc = temperatura pseudocritica, °R
ppc = presión pseudocritica, psia
T’pc = Temperatura pseudocritica corregida, °R
P’pc = Presión pseudocritica corregida, psia
B = fracción molar del H2S en la mezcla de gas
ε = Factor de temperature pseudocrítica, definido como:
A= yH2S + yCO2
Donde
T’pc = temperatura pseudocritica, °R
Tpc = temperatura pseudocritica no ajustada, °R
yCO2 = fracción molar de CO2
yH2S = fracción molar de H2S
yN2 = fracción molar de N2
p’pc = presión pseudocritica ajustada, psia
ppc = presión pseudocritica no ajustada, psia
Paso 3. Usando los valores pseudocríticos ajustados calcular las propiedades
pseudoreducidos.
Paso 4. Calcular el factor Z de la gráfica de Standing Katz o por correlaciones.
El Método de Hall-Yarborough
Hall y Yarborough (1973) presentaron una ecuación de estado que representa en forma
cercana al gráfico de Z de Standing y Katz. La expression propuesta se basa en la ecuación de
estado de Starling-Carnahan. Los coeficientes de la correlación han sido determinados de un
ajuste de los mismos con los datos del gráfico de Standing Katz.
La ecuación de estado de Starling presenta la siguiente forma
(11)
donde, las constantes Ai tienen los siguientes valores:
A1 = 0,3265
A2 = -1,0700
A3 = -0,5339
A4 = 0,01569
A5 = -0,05165
A6 = 0,5475
A7 = -0,7361
A8 = 0,1844
A9 = 0,1056
A10 = 0,6134
A11 = 0,7210
Reemplazando s ρ r por su expresión en la ecuación (11) se tiene:
(13)
donde,
(14)
(15)
(16)
(17)
(18)
F = A11.(0,27.sP r/sT r) ² (19)
La ecuación reajustada de Yarborough y Hall es la siguiente:
Donde:
ppr = presión pseudoreducida
t = recíproco de la temperature pseudoreducida, Tpc/T
Y = densidad reducida que puede obtenerse de la siguiente expresión:
Donde:
Esta ecuación no lineal puede ser convenientemente resuelta gracias al método iterativo de
Newton-Raphson
El procedimiento computarizado para resolver la ecuación a condiciones específicas de ppr y
Tpr de acuerdo a los siguientes pasos:
Paso 1. Asumir un valor inicial del parámetro Yk, donde k es el contador de la iteración, de
acuerdo a la siguiente relación:
Paso 2. Sustituya el valor inicial en la ecuación y evaluar la función no lineal, hasta que el valor
correcto de Y haya sido seleccionado, la ecuación correspondiente tendrá un valor determinado
de F(Y)
Paso 3. Un Nuevo valor estimado de Y, como Yk1 es cálculado con la siguiente expresión:
Paso 4. Los pasos 2 y 3 son repetidos n veces hasta tener un valor de error menor al
predeterminado o también denominado tolerancia.
Paso 5. El valor corregido de y es usado para evaluar la ecuación de z.
Hall y Yarborough puntualizaron que este método no es susceptible de ser aplicado cuando la
temperatura pseudoreducida es menor a uno.
El Método de Dranchuk-Abu-Kassem
Dranchuk y Abu-Kassem (1975) derivaron una expresión analítica para calcular la densidad
pseudoreducida que puede ser usada para estimar el factor de compresibilidad del gas. La
densidad reducida del gas es definida como la razón de la densidad de gas a una presión y
temperatura específica:
A A 5 2
7 8 2 r 2
A9 r A10 (1 A11 r ) 3 exp( A11 r )
T T 2 Tr
r r
Donde los valores numéricos de sus constantes están dados por:
A1 = 0.3265 A5 = -0.05165 A9 = 0.1056
A2 = -1.0700 A6 = 0.5475 A10 = 0.6134
A3 = -0.5339 A7 = -0.7361 A11 = 0.7210
A4 = 0.01569 A8 = 0.1844
El factor de compresibilidad crítico zc es aproximadamente 0.27 que conduce a la siguiente
ecuación simplificada:
Paso 2. Sustituya este valor inicial en la ecuación correspondiente para evaluar la función no
lineal, hasta que se determine el valor correcto de densidad pseudoreducida.
Paso 3. Un Nuevo valor estimado de densidad pseudoreducida, como ρrk 1, es calculado de
la siguiente expresión:
Paso 4. Los pasos 2 y 3 se repiten n veces hasta que el error sea menor a la tolerancia
Paso 5. El valor corregido de la densidad es empleado en la ecuación del factor de
compresibilidad:
El Método de Dranchuk-Purvis-Robinson
Dranchuk, Purvis, y Robinson (1974) desarrollaron una correlación basada en la ecuación de
estado de Benedict-Webb-Rubin. Con una aproximación efectuada en base a 1500 datos del
gráfico de Standing and Katz, optimizando ocho de los coeficientes de la ecuación propuesta:
Donde se tiene:
JUSTIFICACION
la producción de los pozos, entonces el problema se formula como la minimización del costo
de la producción sujeto a las restricciones técnicas debidas a la topología y a la demanda en
sí. Para evitar problemas de almacenamiento, el sistema no produce mas que lo
demandado.
En nuestro caso nos limitamos solamente a determinar el algoritmo y el programa para el
calculo de estos parámetros pseudo-critico, factor z, peso molecular aparente, gravedad
especifica y densidad con el objetivo de determinar la magnitud de desviación de los gases
reales en función a su temperatura y presión.
OBJETIVOS
MARCO TEORICO
Estudios del factor de compresibilidad de los gases para gases naturales de varias
composiciones han demostrado que el factor de compresibilidad puede ser generalizado con
suficiente precisión para estudios de ingeniería cuando estos son expresados en términos de
las siguientes propiedades adimensionales:
presión- pseudo reducida
Temperatura-pseudo reducida
Estos términos adimensionales son definidos por las siguientes expresiones:
Donde
P = presión del sistema, psia
ppr = presión-pseudo reducida, adimensional
T = temperatur del sistema, °R
Tpr = temperatura-pseudo reducida, adimensional
ppc, Tpc = presión y temperatura-pseudo critica, respectivamente, y
definido por la siguiente relación:
Sin embargo estas propiedades pseudo-criticas, i.e., ppc y Tpc, no representan las
propiedades críticas actuales de la mezcla de gas.
Estas propiedades son usadas como parámetros de correlación dentro de la generación de
propiedades de gas.
En casos donde la composición del gas natural no es disponible, las propiedades pseudo-
criticas i.e., ppc y Tpc, pueden ser predecidas únicamente de la gravedad especifica de los
gases. Brown et al. (1948) presento un método grafico para una aproximación conveniente
de la presión y temperatura pseudo-critica, respectivamente, de gases cuando se cuenta
solamente con la gravedad especifica de los gases. La correlación se dada en la figura 2-
2.Standing (1977) expreso esta correlación grafica en las siguientes formas matemáticas:
Donde
Tpc = temperatura pseudo-critica, °R
ppc = presión pseudo critica, psia
gg = gravedad especifica de la mezcla de gas
1. Objetivo de La Aplicación:
Realizar el diseño para dimensionar separadores de tres fases (horizontales y verticales), en función a
parámetros conocidos de análisis previos de laboratorio, datos de producción estimaciones.
2. Conceptos Fundamentales:
a) Fuerza centrífuga
b) Fuerza de gravedad
c) Fuerza de choque
3. Diseño:
e) Factor de compresibilidad Z:
(ver subrutinas, 3.5)
h) Constante de separación:
b) Capacidad de líquido:
R = (Lcc * 12) / D
La capacidad de gas está basada en la velocidad de gas, v, que debe ser menor a la
velocidad final de la partícula, Vt.
a) Capacidad de gas:
El diámetro mínimo que satisface la capacidad de gas está expresado por la ecuación:
El diámetro mínimo que satisface el asentamiento de las gotas de líquido esta dada
por la ecuación:
La altura del líquido (agua + petróleo) para los tiempos de retención se puede
determinar mediante la ecuación: para un diámetro interno dado D:
R = (Lcc / D)*12
3.5. Subrutinas:
a) Factor de Compresibilidad Z:
Esta es una ecuación no lineal y puede ser resuelta por el método de Newton
Raphson (técnica interactiva):
f ( x)
Yx 1 Yx
f ' ( x)
1 T1 r T2 r2 T3 r5 T4 r2 (1 A8 r2 )e ( A8 r )
2 T5
r
0
con
A A
T1 A1 2 33
T pr T pr
A
T2 A4 5
T pr
A A
T3 5 6
T pr
A
T4 37
T pr
0.27 p pr
T5
T pr
donde ρr se define por la Ecuación 5.42 y los coeficientes A1 al A8 tienen los
siguientes valores:
A1 = 0.31506237 A3 = -0.57832720 A5 = -0.61232032 A7 = 0.68157001
A2 = -1.0467099 A4 = 0.53530771 A6 = -0.10488813 A8 = 0.68446549
El procedimiento de solución de la Ecuación es similar a la de Dranchuk y
Abu-Kassem.
El método es válido dentro de los siguientes rangos de temperatura y presión
seudo reducida:
1.05 ≤ Tpr < 3.0
0.2 ≤ ppr ≤ 3.0
Para resolver este método se debe emplear dos tipos de iteración usando el método
de Newton - Raphson.
f ( rk )
rk 1 rk
f ' ( rk )
donde
f ' ( r ) ( R1 ) 2( R2 ) r 5( R3 ) r4 2( R4 ) r e A11 r 1 2 A11 r2 (1 A11 r2
2
Paso 4. Los pasos 2 y 3 se repiten n veces, hasta que el error, o sea, |ρrk-
ρrk+1|, llegue a ser más pequeño que una tolerancia pre-establecida, por
ejemplo de 10-12.
Paso 5. El valor correcto de ρr es entonces usado para evaluar la Ecuación
para el factor de compresibilidad, esto es:
0.27 p pr
z
r T pr
La correlación propuesta se reportó que duplicó los factores de
compresibilidad de la gráfica de Standing y Katz con un error promedio
absoluto de 0.585 % y es aplicable a los rangos:
0.2 ≤ ppr < 30
1.0 < Tpr ≤ 3.0
donde
(9.4 0.02 M a )T 1.5
K
209 19 M a T
986
X 3 .5 0.01M a
T
Y 2 .4 0 .2 X
ρg = densidad del gas a presión y temperatura del reservorio, lb/ft3
T = temperatura del reservorio, °R
Ma = peso molecular aparente de la mezcla del gas
c) Constante de separación:
ρ = ρl - ρg
- Número de Reynols:
Nre = (0.0049 * ρg * v) / μg
- Factor de fricción:
case contrario:
- Constante de separación:
4. Aplicación Computarizada:
5. Análisis de Resultados:
De acuerdo a los resultados adjuntos obtenidos para datos de prueba, se debe tomar en
cuenta que para separadores verticales, según Lockhart (1986) establece dentro de otros
factores que la relación Longitud / Diámetro debe estar en el rango de 1.5 a 3. En tai sentido
se puede escoger razonablemente un separador vertical entre 84” y 96”, dependiendo de las
disponibilidades en la industria.
En forma similar, para separadores horizontales, se debe tomar en cuenta la relación Longitud
/ Diámetro en el rango de 3 a 5, por lo que se puede utilizar separadores entre 72” y 84” de
diámetro, dependiendo de las disponibilidades de la industria, sin embargo el diámetro del
separador debe ser menor que el diámetro máximo calculado.
6. Bibliografía:
- Bizanti and Yu Hancheng, “Program Acelerates Three-phase Oil-Gas Separator Desing”,
Petroleum Engineer,May/1988
- Sanjay Kumar, “Gas Production Engineering”, Gulf Publishing Co., 1987.
- Tarek Ahmed, “Hydrocarbon Phase Behavior”, Gulf Publishing Co., 1989.
Simbología
RESULTADOS
CAPACIDAD DE LIQUIDO
El mercado del gas y sus derivados, en forma directa como gas al usuario o
en forma de líquido embotellado que sale como gas, tiene sus características
propias, modalidades y normas para su utilización.
En resumen, las operaciones de exploración, perforación, producción,
transporte y procesamiento del gas se han convertido en una importantísima
industria dentro de la industria petrolera global.
Las conducciones de transporte llegan a tener diámetros entre 42 y 48
pulgadas (unidad aceptada internacionalmente para esta industria), equivalentes a
1 y 1,20 m, mientras que las de distribución oscilan entre 18 y 22 pulgadas (40 y 70
cm) para los cuales se hace importante el calculo de los caudales de gas que se
van a transportar por el gasoducto ya sea que tenga ramales o no, lo importante es
saber cuando flujo de gas va poder ser atravesado hasta llevarlo al final de entrega
del gasoducto y para esto se ha hecho necesario encontrar ecuaciones de flujo
horizontal a través de tuberías que se mencionan en el presente trabajo, las cuales
son las ecuaciones de Panhandle A, Panhandle B y de Weymouth, mismas que han
sido y son necesitadas para el calculo del flujo de gas en gasoductos que
actualmente están llevando el gas.
MARCO TEORICO
¿A qué se conoce como gas natural?
A la mezcla de hidrocarburos que comúnmente se emplea para propósitos
energéticos y que, por lo general, se utiliza para fines domésticos e industriales.
Es conveniente aclarar que - en su composición - no aparecen únicamente
los hidrocarburos, sino también las impurezas (ver Figura1), como el agua, el dióxido
de carbono y el sulfuro de hidrógeno, citados a título de ejemplo. Adicionalmente,
el personal que trabaja en este tipo de operaciones debe vigilar la presencia de
arena, que produce la erosión. Las parafinas y los asfaltenos se depositan y crean
problemas. Cuando el agua está en forma líquida y en presencia de sulfuro de
hidrógeno (H2S), forma ácidos que corroen las instalaciones. El ingeniero encuentra
su razón de ser en la solución de estos problemas.
Figura 5. Esquema de la balanza de Edward, utilizada para medir la gravedad específica de los gases
Factor de Compresibilidad “z”
Una de las características de los gases es que al aplicarles presión pueden ser
comprimidos y, por ende, pueden ser almacenados o confinados en recipientes de
determinados volúmenes.
Las condiciones o postulados en que se basa la teoría cinética de los gases
no se pueden cumplir y la situación en que más se aproximan a ellas es cuando la
presión y la temperatura son bajas; cuando éstas son altas el comportamiento del
gas se aleja de tales postulados, especialmente en lo relacionado a que no hay
interacción entre las moléculas de tipo gravitacional, eléctrica o electromagnética
y a que el volumen ocupado por las moléculas es despreciable comparado con el
volumen total ocupado por el gas; en este caso no se habla de gases ideales sino
de gases reales.
Como el gas real no se ajusta a la teoría cinética de los gases tampoco se
ajusta a la ecuación de estado y se hace necesario establecer una ecuación de
estado para gases reales.
Las relaciones de composición, presión, volumen y temperatura detalladas
antes e incluidas en la fórmula que define la ley sobre gases perfectos, todavía no
está completa porque falta tomar en cuenta el factor de compresibilidad (Z).
El físico Juan Van Der Waals (1837- 1923), estudió la atracción molecular y el
tamaño de las moléculas de los gases e introdujo en la fórmula el factor de
corrección, para que en su forma final la ecuación más sencilla y la más conocida
para analizar el comportamiento de los gases reales presenta la siguiente forma:
PV = znRT (1)
P: presión absoluta.
V: volumen.
R: constante universal de los gases.
T: temperatura absoluta.
De manera que para un determinado gas y n = 1:
z = PV
RT
Z se puede considerar como un factor de corrección para que la ecuación
de estado se pueda seguir aplicando a los gases reales. En realidad Z corrige los
valores de presión y volumen leídos para llevarlos a los verdaderos valores de presión
y volumen que se tendrían si el mol de gas se comportara a la temperatura T como
ideal. Z se conoce como factor de supercompresibilidad, y depende del tipo de gas
y las condiciones de presión y temperatura a que se encuentra; cuando éstas son
bajas, próximas a las condiciones normales, Z se considera igual a uno.
Cuando se trata de gases reales, la presión indicada por el registrador de
presión es menor que la presión a la que se encontraría el gas si fuera ideal pues hay
que descontar las interacciones entre las moléculas y por otra parte el volumen
disponible para el movimiento de las moléculas es menor que el volumen del
recipiente pues no se puede despreciar el volumen ocupado por las moléculas.
Z es adimensional y depende de las presiones y temperaturas a las que sea
sometido el gas. Por tanto, valores de Z pueden determinarse por experimentación.
De allí que en la industria existen catálogos, tablas y manuales de consultas sobre
infinidad de muestras y análisis del gas natural.
Sin embargo, a través del conocimiento de la temperatura y presiones
críticas, determinadas por experimentos, correspondientes a cada uno de los
componentes que forman el gas natural se pueden calcular presiones y
temperaturas “reducidas” que facilitan la obtención de supuestas “seudo presión
crítica” y “seudo temperatura crítica” para tomar en consideración la contribución
porcentual de cada componente, de acuerdo a la composición del gas.
El siguiente ejemplo hipotético servirá para calcular el factor de
compresibilidad.
Mezclas de Gases Reales
(3)
Para calcular la densidad (ρ) se aplica la ecuación:
ρ = P*M/z*R*T = m/V (4)
temperatura se le llama presión crítica, que a la vez representa la presión más alta
que los valores del líquido pueden ejercer.
Una vez conocidas las condiciones críticas se obtienen las condiciones
reducidas, que se definen como:
Pr = P/Pc (5)
Tr = T/Tc (6)
Donde:
Pr: presión reducida.
Tr: temperatura reducida.
Como se ve son adimensionales.
Los cálculos para el ejemplo dado muestran que la pseudo temperatura
crítica dio 198 °K (columna E) y la pseudo presión crítica resultó ser 45,78 atms. abs.
Si se desea obtener el factor de compresibilidad del gas en cuestión, a
determinada presión y temperatura, entonces se procede a calcular los valores de
presión y temperatura reducidas, Pr y Tr. Sea el caso que se desee conocer el valor
de Z a temperatura de 44 °C y a presión de 50 atms. abs.
50 317
Pr 1,90 Tr 1,60
45.78 198
Con estos dos valores se recurre a un gráfico de pseudo temperatura
reducida y pseudo presión reducida para determinar el valor de z = 0,90. Para esto
ver el Anexo 1 que muestra el grafico para obtener el factor de corrección z
utilizando valores de seudo presión y seudo temperatura reducidos, calculados
previamente.
b) Obtención de Z para mezclas
También se utiliza la correlación de Standing pero en este caso las
condiciones reducidas no se pueden obtener de tablas porque las mezclas no son
compuestos puros, además cuando se trata de mezclas no se habla de
condiciones críticas o reducidas sino de condiciones pseudocríticas y
pseudoreducidas.
Para obtener las condiciones pseudocríticas se debe conocer la composición
de la mezcla o la gravedad específica.
Cuando se tiene la composición se puede aplicar el procedimiento de Kay
para obtener las condiciones pseudocríticas.
1. Procedimiento de Kay
Se aplica la siguiente ecuación:
sP c = Σxi.P ci (7)
sT c = Σxi.T ci (8)
Donde:
sP c: presión pseudocríticas de la mezcla.
sT c: temperatura pseudocríticas de la mezcla.
xi: fracción molar de cada componente en la mezcla.
P ci: presión crítica de cada componente en la mezcla.
T ci: temperatura crítica de cada componente en la mezcla.
Una vez calculados los valores de sT c y sP c, se calculan las condiciones
pseudoreducidas:
sP r = P/sP c (9)
sT r = T/sT c (10)
Donde:
sP r: presión pseudoreducida de la mezcla.
sT r: temperatura pseudoreducida de la mezcla.
Ecuación de Starling
La ecuación de estado de Starling presenta la siguiente forma:
(11)
Donde:
s ρ r: se conoce como densidad pseudoreducida y está dada por:
sρr = 0,27.sP r/Z.sT r (12)
Las constantes Ai tienen los siguientes valores:
A1 = 0,3265 A2 = -1,0700 A3 = -0,5339 A4 = 0,01569
A5 = -0,05165 A6 = 0,5475 A7 = -0,7361 A8 = 0,1844
A9 = 0,1056 A10 = 0,6134 A11 = 0,7210
Reemplazando s ρ r por su expresión en la ecuación (11) se tiene:
(13)
Donde:
(14)
(15)
(16)
(17)
(18)
(20)
- Se supone un valor de Z(Z0), se recomienda mayor que 1, y se chequea si
hace F(Z0) = 0 dentro de la tolerancia requerida.
- Si F(Z0) = 0, el valor supuesto es el correcto y es el valor que se está
buscando; si F(Z0) ≠ 0 se busca un nuevo valor de Z(Z1) de la siguiente manera:
Z1 = Z0 - F(Z0)/F'(Z0) (21)
Donde F'(Z0) es la derivada de F(Z), dada por:
(22)
y calculada en Z0.
- Con Z1 se chequea si F(Z1) = 0 y si no lo es se calcula un valor Z 2 usando la
ecuación (22) cambiando Z1 por Z2 y Z0 por Z1.
- El procedimiento continua hasta encontrar un valor Z n que haga F(Z) = cero;
después del primer valor supuesto para Z(Z0) los demás valores usados se obtienen a
partir de la ecuación (21) usando Zn en lugar de Z1 y Zn-1 en lugar de Z0.
Poder calorífico del gas natural
Una de las características del gas natural es su poder calorífico, el cual se
determina por análisis de laboratorio, utilizando uno de los varios tipos de
calorímetros disponibles. Además, el poder calorífico del gas se considera para
determinar su calidad como combustible y, por ende, su precio.
En el sistema angloamericano a la unida de caloría se le llama Unidad
Térmica Británica (BTU) y se define como la cantidad de calor requerida para
aumentar la temperatura de 1 libra (453,592 gramos) de agua a un grado Fahrenheit
hasta la temperatura de su máxima densidad que es 39,2 °F. Una BTU es,
aproximadamente, igual a 0,252 kilocalorías.
El gas natural puede tener de 8.000 a 11.115 kilocalorías/metro cúbico, lo que
equivale a 900 y 1.250 BTU/pie cúbico, respectivamente.
El petróleo crudo tiene poder calorífico que va de 8.500 a 11.350 calorías por
gramos o 15.350 a 22.000 BTU por libra.
Así que, por medio del poder calorífico del gas natural en general o de
sus componentes en particular, y el poder calorífico de los crudos, es posible hacer
cálculos que permiten determinar que tantos metros cúbicos o pies cúbicos de gas
equivalen a un metro cúbico o barriles de petróleo. Este tipo de equivalencia es de
referencia común en la industria.
Específicamente, el precio que se le asigna a determinado gas se basa en
una unidad de volumen: metro cúbico o pie cúbico. Sin embargo, como los
volúmenes de entrega por lo general son muy grandes se opta por el millar de
metros o pies cúbicos. También se emplea el poder calorífico, expresado en millones
de calorías o de BTU. En el caso de gases licuados, en vez del volumen o del poder
calorífico, se hace referencia al peso en kilos o libras.
Viscosidad del gas natural
Así como la viscosidad es una característica física importante de los líquidos,
también lo es para los gases. La unidad de medida en ambos casos es el poise, en
honor al médico y físico francés J.L.M. Poiseuille († 1869).
La definición de poise se deriva de la determinación de la fuerza requerida
por centímetro cuadrado para mover a velocidad de un centímetro por segundo
un plano móvil y paralelo a otro plano fijo distantes un centímetro entre sí y cuyo
espacio está lleno del líquido o fluido objeto de la medición de viscosidad.
La viscosidad del gas natural es expresión de su resistencia al flujo y tiene
aplicaciones importantes en la producción, procesos de acondicionamiento y
mercadeo.
- No se requiere el almacenamiento.
- Se distribuye por ductos subterráneos de polietileno y acero, materiales probados
en zonas sísmicas: México, Tokio, San Francisco, etc.
- Todas las instalaciones cumplen las Normas Oficiales Mexicanas.
- Las distribuidoras de gas natural supervisan y monitorean constantemente las redes
de distribución y cuentan con equipos técnicos disponibles las 24 horas, los 365 días
del año.
Ecológico
- Comparado con otros hidrocarburos, poseé la menor relación de hidrógeno-
carbón en su composición química, por ello su combustión es más limpia y la que
menos emisiones contaminantes libera al ambiente.
- La combustión de gas natural no produce residuos contaminantes.
- No genera partículas sólidas ni emite residuos tóxicos.
- El protocolo de Kyoto reconoce al gas natural como el combustible fósil más
amigable con el medio ambiente.
- Son innecesarias instalaciones de almacenamiento masivo, ni vehículos
transportadores y/o repartidores.
Económico
- Es el combustible más económico del mercado, no importa el destino del
consumo, con gas natural solo paga lo que consume.
Comodidad
- El gas natural llega por tubería directamente a los hogares, al comercio y a la
industria.
Tranquilidad
- No tiene que cargarse ni recarga incómodos ni pesados cilindros.
USOS
El gas natural tiene diversas aplicaciones en la industria, el comercio, la
generación eléctrica, el sector residencial y el transporte de pasajeros. Ofrece
grandes ventajas en procesos industriales donde se requiere de ambientes limpios,
procesos controlados y combustibles de alta confiabilidad y eficiencia.
En el siguiente cuadro se presentan algunas de las aplicaciones más comunes
de gas natural:
Sector Aplicaciones/Procesos
Industrial Generación de vapor
Industria de alimentos
Secado
Cocción de productos
cerámicos
Fundición de metales
Tratamientos térmicos
Temple y recocido de
metales
Generación eléctrica
Producción de
petroquímicos
Sistema de calefacción
Hornos de fusión
Comercio y Servicios Calefacción central
Aire acondicionado
Cocción/preparación de
alimentos
Agua caliente
Energía Cogeneración eléctrica
Centrales térmicas
Residencial Cocina
Calefacción
Agua caliente
Aire acondicionado
Transporte de pasajeros Taxis
Buses
TABLA Combustibles que el Gas Natural puede sustituir
Adicionalmente, el gas natural es utilizado como materia prima en diversos
procesos químicos e industriales. De manera relativamente fácil y económica puede
ser convertido a hidrógeno, etileno, o metanol; los materiales básicos para diversos
tipos de plásticos y fertilizantes.
Uno de los usos del gas natural que ha tomado mayor importancia en el
último tiempo en nuestro país, dada la disminución en el recurso hídrico, es la
generación eléctrica.
Debido a que el gas natural puede ser utilizado con grandes beneficios en un
amplio número de aplicaciones, puede sustituir a los energéticos alternativos que se
señalan a continuación:
TABLA Combustibles que el Gas Natural puede sustituir
Sector Energía y/o combustible que puede sustituir
Industrial Carbón
Electricidad
Diesel
Fuel Oil
Gas licuado
Gasolina
Kerosene
Leña
Generación eléctrica Carbón
Fuel Oil
Comercio Carbón
Electricidad
Fuel Oil
Gas natural
Gas licuado
Kerosene
Residencial Electricidad
Gas natural
Gas licuado
Kerosene
Leña
Transporte de pasajeros Gasolina
Petróleo Diesel
RESERVAS MUNDIALES DE GAS NATURAL
Evolución de las reservas mundiales probadas de gas natural por zonas (Billones de m 3)
Fuente: SEDIGAS
Con los datos disponibles hoy en día las reservas evaluadas de gas natural son
suficientes para abastecer al mundo, con un consumo como el de 2005, durante
más de 65 años.
Reservas mundiales probadas de gas natural (Billones de m 3)
Fuente: SEDIGAS
2004 2005 2006
América del Norte 7 7 7
América Central y Sur 7,4 7,3 7,3
Europa-OCDE 6,3 6,2 6,5
Europa Oriental 57,8 57,8 57,8
África 13,9 14,1 14,4
Oriente Medio 71,5 73,3 74,6
Asia-Oceanía 13,7 14,3 14,5
TOTAL MUNDIAL 177,6 180 182,1
CONSUMO DE GAS NATURAL
La parte final del manejo del gas la constituye el transporte desde las
instalaciones de los campos y las entregas de volúmenes determinados a los
mercados en ruta.
Estas dos fases representan en la práctica el mercadeo y la comercialización
del gas. De acuerdo con las modalidades mundiales para este tipo de operaciones
cabe mencionar aspectos interesantes:
• Se da el caso de que existen empresas integradas cuyas operaciones
(exploración, perforación, producción, transporte y mercadeo) están dedicadas
exclusivamente al gas y no producen petróleo. Son empresas especializadas en el
negocio del gas.
• Existen otras empresas integradas que se dedican mayoritariamente al
petróleo y que pueden disponer de grandes volúmenes de gas asociado y de gas
libre que las pueden inducir a comercializar el gas parcialmente o totalmente. Esto
es que venden su gas a otras empresas y no se ocupan del mercadeo o podría
optar por transportar, distribuir y vender gas directamente.
• Hay casos en que el gas lo manejan varias empresas. Primero, la que lo
produce y acondiciona. Segundo, la que lo transporta y es dueña del sistema de
gasoductos, y tercero, la que se encarga de la distribución y venta del gas en
determinados mercados de su competencia.
GASODUCTO
Definicion
Un gasoducto es un conducto que transporta o transmite gas natural, en
general a largas distancias y grandes volúmenes y cuya presión de diseño es igual o
mayor a 40 bares
Un gasoducto es una conducción que sirve para transportar gases
combustibles a gran escala. Es muy importante su función en la actividad
económica actual.
Impropiamente, y puede que por analogía con el oleoducto, se le llama con
frecuencia gaseoducto.
Construcción
Consiste en una conducción de tuberías de acero, por las que el gas circula a
alta presión, desde el lugar de origen. Se construyen enterrados en zanjas y se
entierran a una profundidad típica de 1 metro. Excepcionalmente, se construyen
sobre la superficie. Si la distancia es larga, puede haber estaciones de compresión a
intervalos.
Por razones de seguridad, las regulaciones de todos los países establecen que
a intervalos determinados se sitúen válvulas en los gasoductos mediante las que se
pueda cortar el flujo en caso de incidente. Además, si la longitud del gasoducto es
importante, pueden ser necesarias estaciones de compresión a intervalos.
El inicio de un gasoducto puede ser un yacimiento o una planta de
regasificación, generalmente situada en las proximidades de un puerto de mar al
que llegan buques (para el gas natural, se llaman metaneros) que transportan gas
natural licuado en condiciones criogénicas a muy baja temperatura (-161 ºC).
Para cruzar un río en el trazado de un gasoducto se utilizan principalmente
dos técnicas, la perforación horizontal y la perforación dirigida. Con ellas se consigue
que tanto la flora como la fauna del río y de la ribera no se vean afectadas. Estas
técnicas también se utilizan para cruzar otras infraestructuras importantes como
carreteras, autopistas o ferrocarriles.
El tendido por mar se hace desde barcos especialmente diseñados, los cuales
van depositando sobre el lecho marino la tubería una vez soldada en el barco.
Regulaciones estatales en muchos países requieren que los gasoductos
enterrados estén protegidos de la corrosión. A menudo, el método más económico
es revestir el gasoducto con algún tipo de polímero de modo que la tubería queda
eléctricamente aislada del terreno que la rodea. Generalmente se reviste con
pintura y polietileno hasta un espesor de 2-3 mm. Para prevenir el efecto de posibles
fallos en este revestimiento, los gasoductos suelen estar dotados de un sistema de
protección catódica, utilizando ánodos de sacrificio que establecen la tensión
galvánica suficiente para que no se produzca corrosión.
El impacto ambiental que producen los gasoductos, se centra en la fase de
construcción. Una vez terminada dicha fase, pueden minimizarse todos los impactos
asociados a la modificación del terreno, al movimiento de maquinaria, etc. Queda,
únicamente, comprobar la efectividad de las medidas correctivas que haya habido
que tomar en función de los impactos, repoblaciones, reforestaciones, protección
de márgenes, etc.
Unidades Constructivas que componen la construccion de un gasoducto
- Replanteo
- Apertura de pista de trabajo (desbroce y movimiento de tierras).
- Excavación de la zanja.
- Distribución y manipulación de la tubería.
- Soldadura y radiografiado.
- Revestimiento de juntas de soldadura.
- Puesta en zanja.
- Pretapado, tapado y restitución de los terrenos.
- Cruces especiales.
- Pruebas hidráulicas.
- Restitución de los terrenos
Circulación del gas
La presión a la que circula en gas por el gasoducto es normalmente de 72 bar
para los de la redes básicas de transporte y 16 bar en las redes de distribución.
Para llevar el gas hasta los hogares y comercios, es preciso bajar la presión de
transporte hasta límites razonablemente seguros. Esto se consigue instalando
estaciones de regulación a lo largo del gasoducto en las que se baja la presión
hasta la presión habitual de distribución.
El cambio de presiones se hace de forma análoga a las redes eléctricas (alta
tensión/baja tensión), en este caso se utilizan estaciones de regulación y medida,
por medio de reguladores de presión de membrana se regula la presión de salida
que se necesite.
Obstáculos en el gasoducto
Los grandes exportadores no pueden acaparar el mercado global de gas
porque, para su mala suerte, éste no existe.
"A largo plazo vamos hacia una OPEP del gas", dijo Chakib Khelil, ministro de
energía de Argelia en una reciente reunión de grandes productores en Qatar. Es una
idea espeluznante, pues el suministro mundial de gas está concentrado en aún
menos países que el de petróleo. Qatar, Rusia e Irán controlan casi 60% y, con
excepción de Qatar, no están entre los aliados más confiables de Occidente (ver
gráfica).
Sin embargo, una crisis del gas similar a la del petróleo ocurrida en los años 70
es un escenario improbable por una razón: crear el citado cártel llevaría mucho
tiempo. Hasta ahora sus eventuales miembros se han conformado con elaborar un
vago aunque ligeramente inquietante "estudio de precios". Aun si se creara una
organización de países exportadores de gas (OPEG), existen varias formas de
generar electricidad y calefacción para los hogares; así, si el gas se vuelve
demasiado caro, los consumidores pueden recurrir a alternativas. En contraste, la
gasolina es el único combustible con el cual funciona la mayoría de los autos, por lo
cual la demanda de petróleo seguirá siendo relativamente inelástica a largo plazo.
Afección Medioambiental de un Gasoducto
Cualquier canalización de gas, al ser una infraestructura enterrada, no tiene
incidencia apreciable sobre el medio ambiente, limitándose su afección casi
exclusivamente al período de construcción, y siendo extremadamente reducido si se
compara con el de otras obras lineales, tales como carreteras, vías férreas o tendido
de líneas eléctricas.
Las principales afecciones que puede producir el trazado del gasoducto
sobre el medio natural durante la fase de obras son:
- Impactos sobre la vegetación
- Impactos sobre la fauna
- Impactos sobre los cursos de agua
- Impacto sobre zonas húmedas
- Impactos sobre el Patrimonio Cultural
- Impactos sobre el paisaje
a) Impactos sobre la vegetación
La vegetación se verá afectada principalmente por el despeje y desbroce en
la apertura de la pista de trabajo. La afección sobre la flora, depende de la
facultad, por la cobertera vegetal, de reconstruirse después de finalizados los
trabajos. Esto, salvo en el caso de árboles de tallo o tronco alto, es siempre posible.
Cuando esto no fuera posible, se minimiza la afección replantando especies
adecuadas preservando únicamente el pasillo inmediato a la conducción (4 metros
aproximadamente).
b) Impactos sobre la fauna
La afección sobre la fauna se limita a las molestias producidas durante el
período de construcción, es decir a una temporada o estación, por la presencia de
personal y funcionamiento de maquinarias en lugares donde la fauna no esta
acostumbrada a encontrarlos.
Mediante la aplicación de medidas preventivas como la reducción al mínimo
de presencia de personal, elección de fechas de intervención en función de las
costumbres de vida de la especie, etc, se reduce al máximo y en medida de lo
posible la duración de la obra, por lo que la afección a la fauna se verá disminuida.
c) Impactos sobre los cursos de agua
- Afección sobre la fauna piscícola durante los trabajos: Generalmente suelen
tomarse medidas adecuadas de acuerdo con la Administración, Organismos
Oficiales o Asociaciones de pesca para reducir o incluso suprimir la posible afección.
Cada travesía de río o curso de agua presenta normalmente circunstancias
particulares que hay que analizar específicamente.
- Incremento de los sólidos en suspensión en las aguas superficiales como
consecuencia de las obras: Durante la fase de construcción podría producirse un
aumento de los sólidos en suspensión en los cursos de agua superficial cercanos
Longitud
Para ductos de grandes diámetros y longitudes y elevadas presiones una de
las ecuaciones que mejor se aproxima al comportamiento del gas es la ecuación
Panhandle- B, en unidades inglesas:
0.51
2 0.0375GE(h1 h2 ) Pav2
To
1.02 P1 P2
2
Q 737Eff Di2.53
zavTav
0.961
Po
GE zavTav L
La presión promedio es:
2 P *P
Pav P1 P2 1 2
3 P1 P2
Ecuacion de Panhandle A
0.5394
2 0.0375* GE* (h2 h1 )Pav2
Tb
1.0778 2 2
P1 P2
Di2.6182
T z
Q 435.87 av av
Eff
Pb GE0.8538 * zav *Tav * L
Donde:
Q= Flujo transportado en Pies cubicos estandar por dia [STBD]
Di= Diametro interno del ducto en pulgadas [inch]
L= Longitud en millas
zav= Factor de compresibilidad del gas natural que es adimennsional calculado a las
condiciones de Presion promedio y Temperatura promedio
TO y Po= Temperatura y Presion a las condiciones base de medicion en [ oR] grados
Rankine y [psia]
P1 y P2= Presiones al inicio y al final del tramo del ducto en [psia]
Pav = Presion promedio en el tramo, es la presion utilizada en el calculo del factor de
compresibilidad
GE= gravedad especifica del gas , es decir, la densidad relativa a la del aire
h1 y h2=alturas sobre el nivel del mar a los extremos del ducto: inicio y final
Tav = Temperatura promedio en el tramo, es la temperatura utilizada en el calculo del
factor de compresibilidad que normalmente es constante en los ductos subterraneos
Eff= Eficiencia de flujo adimensional depende principalmente de la rugosidad y
edad del ducto, generalmente en el sistema de ductos PGPB se tienen una
eficiencia del 85% pues resulta coherente con los datos de campo
Ecuacion de Weymouth
La ecuación de Weymouth se ha utilizado desde hace muchos años como
ecuación de análisis de ductos, y por lo general proporciona resultados muy
conservadores. Los volúmenes mostrados en los diagramas de flujo se expresan en
millones de pies cúbicos medidos a 14.73 psia. y a 60° F.
0.5
Tb 2.667
2
P12 P22 Eff
Q 433 Di
Pb GE * zav * Tav * L
Donde:
Q= Flujo transportado en Pies cubicos estandar por dia [STBD]
Tb y Pb= Temperatura y Presion a las condiciones base de medicion en [ oR] grados
Rankine y [psia]
GE= gravedad especifica del gas , es decir, la densidad relativa a la del aire
zav= Factor de compresibilidad promedio del gas natural
L= Longitud del gasoducto en millas
Di= Diametro interno del ducto en pulgadas [inch]
P1 y P2= Presiones de entrada y salida del tramo del ducto en [psia]
Eff= Eficiencia de flujo adimensional depende principalmente de la rugosidad y
edad del ducto