69

2-. FUNDAMENTOS MATEMTICOS

=
(
(
(
(
(
(

n n nn n n n
n n
n n
n n
b
b
b
b
x a x a x a x a
x a x a x a x a
x a x a x a x a
x a x a x a x a
.
....
. . . . .
...
...
...
3
2
1
3 3 2 2 1 1
3 3 33 2 32 1 31
2 3 23 2 22 1 21
1 3 12 2 12 1 11
........................................................ (2.1)

En la literatura este sistema de ecuaciones tambin se representa en forma
matricial por

b X A = o por b X A = ................................................................................. (2.2)

Donde A ( A ) es la matriz de los coeficientes, X es el vector de las variables o
incgnitas y b es el vector de los trminos independientes.

El sistema de ecuaciones tiene solucin si A es una matriz singular, o sea si su
determinante es diferente de cero.

Los mtodos de solucin para un sistema de ecuaciones lineales pueden ser
mtodos directos o de ensayo y error. Entre los mtodos directos, tambin
conocidos como de eliminacin, se tienen: Mtodo de Cramer, Mtodo de
eliminacin de Gauss y de Jordan - Gauss.

70
2.2-Mtodos para Resolver Sistemas de Ecuaciones Utilizando Matrices
(Mtodos Directos).

1-. Regla de Cramer.

La regla de Cramer dice que cada uno de los valores de las incgnitas de un
sistema de ecuaciones se puede obtener de la siguiente manera:

A
A
A Det
DetA
x
1
1
1
= = .................................................................................
...........................................................................................................
A
A
x
2
2
= .............................................................................................

A
A
x
i
i
= ............................................................................................. (2.3)

donde A es la matriz de los coeficientes del sistema de ecuaciones y A
i
es la
matriz de los coeficientes en la cual se ha cambiado la columna i por el vector de
los trminos independientes del mismo sistema de ecuaciones.
A=
nn n n
n
n
a a a
a a a
a a a
... ...
... ... ... ... ...
... ... ... .... ...
... ...
... ...
2 1
2 22 21
1 12 11
....................................................................
.....................................................................................................................
nn n n
n
n
a a b
a a b
a a b
A
... ...
... ... ... ... ...
... ... ... ... ...
... ...
... ...
2
2 22 2
1 12 1
1
=
71
nn n n
n
n
i
a b a
a b a
a b a
A
... ...
... ... ... ... ...
... ... ... ... ...
... ...
... ...
1
2 2 21
1 1 11
= ................................................................... (2.4)

El mtodo de Cramer se aplica cuando el sistema de ecuaciones es pequeo,
pocas ecuaciones e incgnitas, y como se puede ver el sistema de ecuaciones no
tendr solucin si el determinante de la matriz de los coeficientes es nulo.

2-. Mtodo de Eliminacin de Gauss.

Este mtodo se basa en operaciones que se realizan sobre el sistema original de
ecuaciones, basados en las siguientes propiedades de un sistema de ecuaciones
lineales:
- Las ecuaciones se pueden ordenar arbitrariamente
- El sistema de ecuaciones no cambia si una o varias de las ecuaciones se
multiplica por un escalar.
- El sistema de ecuaciones no cambia si una o varias de las ecuaciones se
reemplazan por la suma de varias de las ecuaciones originales

Se toma la matriz aumentada (la matriz formada por los coeficientes y el vector de
los trminos independientes) que se representa por | | b A. y por medio de
operaciones elementales en esta matriz se llega a una matriz | | ' b U donde U es
una matriz triangular superior o sea que los elementos de por debajo de la
diagonal principal son nulos y b es el vector de los trminos independientes
transformado por las operaciones fundamentales realizadas sobre la matriz | | b A. .

72
(
(
(
(
(
(

=
1 ... 0 0 0
. . . .
... 1 0 0
... 1 0
... 1
x
x
x x
x x x
U

El procedimiento con este mtodo es el siguiente:

Sea un sistema de ecuaciones cuya matriz aumentada es

n nn n n
n
n
b a a a
b a a a
b a a a
Ab
... ... ...
... ... ... ... ... .... ...
... ... ... ... ... ... ...
... ... ... ... ... ... ...
.... .... ....
2 1
2 2 22 21
1 1 12 11
=

i-. Se eliminan los elementos
1 , i
a desde i=2 hasta n de la siguiente manera:

Se reordenan las filas de la matriz de tal manera que en la fila 1 se tenga una
ecuacin cuyo primer coeficiente sea 1. Si esto no es posible se divide la primera
fila por a
11
quedando la matriz de la siguiente forma

n nn n n
n
n
b a a a
b a a a
a b a a a a
Ab
... ... ...
... ... ... ... ... .... ...
... ... ... ... ... ... ...
... ... ... ... ... ... ...
/ / .... .... .... / 1
2 1
2 2 22 21
11 1 11 1 11 12
'
= .......................................... (2.5)

73
A la fila 2 se le resta la nueva fila 1 multiplicada por el coeficiente a
21
; a la fila 3 se
le resta la nueva fila 1 multiplicada por el coeficiente a
13
, y as sucesivamente
hasta la fila n a la cual se le restar la nueva fila 1 multiplicada por a
1n
.

La nueva matriz ser entonces

' ' '
2
'
3
'
3
'
32
'
2
'
2
'
22
'
1
'
1
'
12
) 1 (
... ... ... 0
... ... ... ... ... .... ...
... ... ... ... ... ... ...
... ... ... 0
... ... ... 0
.... .... .... 1
n nn n
n
n
n
b a a
b a a
b a a
b a
Ab = .................................................... (2.6)
donde
11 1
'
1
/ a a a
j j
= ,
11 1
'
1
/ a b b = ,
1
'
1
'
i j ij ij
a a a a = y 1 1 ,
1
'
1
'
> . > = j i a b b b
i i i


ii-. Se eliminan los elementos
'
2 i
a desde i=3 hasta n de la matriz
) 1 (
Ab de la
siguiente manera:

Se divide la fila 2 por el elemento
1
22
a , o se reorganizan las filas de la matriz
) 1 (
Ab
de la fila 2 en adelante, de tal forma que en la fila 2 quede una fila donde el
elemento
1
2 i
a sea 1. En este caso la matriz aumentada queda as:


' ' '
2
'
3
'
3
'
32
'
22
'
2
'
22
'
2
'
22
'
23
'
1
'
1
'
12
) 1 (
... ... ... 0
... ... ... ... ... .... ...
... ... ... ... ... ... ...
... ... ... 0
/ / ... ... ... / 1 0
.... .... .... 1
'
n nn n
n
n
n
b a a
b a a
a b a a a a
b a
Ab =


74
Se eliminan los elementos
'
2 i
a desde i=3 hasta n de la matrz
) 1 (
Ab , restando a
cada fila de la matriz anterior desde i=3 hasta n la fila 2 de la misma matriz anterior
multiplicada por
'
2 i
a . Al recorrer todas las filas desde I=3 hasta n la nueva matriz
quedar as:

2 2 2
3
2
3
2
3
2
33
2
2
2
2
2
23
'
1
1
1
1
12
) 2 (
... ... .. . 0 0
... ... ... ... ... .... ...
... ... ... ... ... ... ...
... ... .. . 0 0
... ... ... 1 0
.... .... .... 1
n nn n
n
n
n
b a a
b a a
b a a
b a
Ab = ....................................................... (2.7)

donde
1
22
1
2
2
2
/ a a a
j j
= ,
1
22
1
2
2
2
/ a b b = y
2
2
1
2
1 2
*
j i ij ij
a a a a = y
2
2
1
12
1 2
*b a b b
i i
= para i=3
hasta n y j=3 hasta n.

iii-. Se eliminan los elementos
2
3 i
a desde i=4 hasta n en la matriz
) 2 (
Ab siguiendo
dos pasos similares a los seguidos para eliminar los elementos
1 i
a de la fila 2 hasta
la fila n en la matriz Ab y los elementos
1
2 i
a de la matriz
) 1 (
Ab de la fila 3 hasta la
fila n. En este caso la matriz que resulta es

3 3 3
4
3
3
3
3
3
34
2
2
2
2
2
23
'
1
1
1
1
12
) 3 (
... ... .. 0 . 0 0
... ... ... ... ... .... ...
... ... ... ... ... ... ...
... ... .. 1 . 0 0
... ... ... 1 0
.... .... .... 1
n nn n
n
n
n
b a a
b a a
b a a
b a
Ab = (2.8)

Donde
2
33
2
3
3
3
/ a a a
j j
= ,
2
33
2
3
3
3
/ a b b = y
3
3
2
3
2 3
*
j i ij ij
a a a a = y
3
3
2
3
2 3
*b a b b
i i i
= para i=4
hasta n y j=4 hasta n

75
iv-. El proceso continua con las filas 4 hasta la n, repitiendo en cada caso los dos
pasos que se indicaron en los numerales i iii. Al terminar el proceso con la fila n
se tendr finalmente la siguiente matriz

n
n
n
n
n
n n
n
n
n
n
b o
b a
b a a
b a a
b a
Ab
1 ... ... .. 0 . 0 0
.. . .. 1 0 0 0 0
... ... ... ... ... ... ...
... ... .. 1 . 0 0
... ... ... 1 0
.... .... .... 1
1
1
1
, 1
3
3
3
3
3
34
2
2
2
2
2
23
'
1
1
1
1
12
) (

= ................................. (2.9)

v-. El sistema de ecuaciones original se ha transformado entonces en
1 ... ... .. 0 . 0 0
.. 1 0 0 0 0
... ... ... ... ... ...
... ... .. 1 . 0 0
... ... ... 1 0
.... .... .... 1
1
, 1
3
3
3
34
2
2
2
23
1
1
1
12
o
a
a a
a a
a
n
n n
n
n
n

*
n
n
x
x
x
x
x
1
3
2
1
...

n
n
n
n
b
b
b
b
b
1
1
3
3
2
2
1
1
...

= ............................ (2.10)

O sea que de acuerdo con las operaciones entre matrices se puede establecer:

n
n n
b x = .......................................................................................................... (2.11)
n
n
n n
n
n n
x a b x *
1
, 1
1
1 1

=
n
n
n n n
n
n n
n
n n
x a x a b x * *
2
, 2 1
2
1 , 2
2
2 2

+ =
= =
n
i j
j
i
ij
i
i i
n i x a b x
1
1 , 1 ................................. (2.12)

La siguiente es una propuesta de cdigo para programar la solucin de un sistema
de ecuaciones por el mtodo de eliminacin de Gauss:

DIMENSION S(30,31), X(30)
76
C S ES LA MATRIZ AUMENTADA Y X ES EL VECTOR DE LAS INCGNITAS
10 READ (8,20) N
20 FORMAT (xx,xx)
C
C
C
NP1=N+1
READ (8,50) ((S(I,J), J=1, NP1)i=1,N)
50 FORMAT (xx,xx)
DO 70 I=1,N
IF S(I,I) .NE.0 THEN 55
ELSE
DO L=I+1,N
IF S(L,I).NE.0 THEN
DO J=1,N
SS(I,J)=S(I,J)
S(I,J)=S(L,J)
S(L,J)=SS(I,J)
END DO
ELSE
END DO
END IF
END IF
55 SII=S(I,I)
DO J=1, NP1
S(I,J)=S(I,J)/SII
END DO
DO 70 K=I+1, N
SKI=S(K,I)
IF (SKI.EQ.0) THEN GO TO 70
DO 70 L=I+1, NP1
70 S(K,L)=S(K,L)-SKI*S(I,L)
77
X(N)=S(N,NP1)
DO I=1,N-1
X(N-I)=S(N-I,NP1)
DO L=N-I+1,N
X(N-I)=X(N-I)-S(N-I,L)*X(L)
END DO
END DO
WRITE (xx,xx) (X(I), I=1,N)
Xx FORMAT ( )
END.

Un caso especial de eliminacin de Gauss se presenta cuando la matriz de los
coeficientes es una matriz tridiagonal en la cual solo los elementos que estn en la
diagonal principal y en las diagonales a la derecha e izquierda de esta son
diferentes de cero, todos los dems son cero. En este tipo de sistema de
ecuaciones el proceso de eliminacin de Gauss es bastante sencillo en cuanto a
los clculos que hay que realizar pues el proceso de eliminacin que hay que
hacer con cada fila restndola de las dems filas siguientes, solo se realiza con la
fila inmediatamente siguiente ya que en las dems filas el elemento a eliminar en
el paso correspondiente es cero Una matriz de este tipo es comn en algunos
mtodos que se aplican para resolver un simulador. Este sistema de ecuaciones
se puede representar matricialmente de la siguiente manera:

n n n n
d
d
d
d
x
x
x
x
b a
c b a
c b a
c b
... ...
... 0 0
... ... ... ... ...
0 0
0 0 0
3
2
1
3
2
1
3 3 3
2 2 2
1 1
= ............................................................. (2.13)

O sea que la matriz aumentada es

78
1 1 1
2 2 2 2
3 3 3 3
0
0 ...
... ... ... ... ... ...
... ... ...
n n n
b c d
a b c d
a b c d
a b d


Al aplicar el primer paso la matriz se transforma en

n n n
d b a
d c b a
a d c a b
... ... ...
... ... ... ... ... ...
... 0
0
1
3 3 3 3
1 2 2 2 1 2 2
1 1
|
|

................................................. (2.14)

donde
1 1 1
/ b c = | y
1 1 1
/ b d = .

Al aplicar el segundo paso la nueva matriz es

n n n
d b a
a d c a b
... ... ...
... ... ... ... ... ...
... 0 0
1 0
1
2 3 3 3 2 3 3
2 2
1 1
|
|
|
................................................ (2.15)

donde
1 2 2
2
2
|
|
a b
c

= y
1 2 2
1 2 2
2
|

a b
a d

=

Cuando se aplique el paso i la matriz resultante ser

79
n n n
i i
d b a ... ... ...
1 0 0 0
... ... ... .... ... ...
1 0
1
2 2
1 1
|
|
|
................................................................... (2.16)

donde
1

=
i i i
i
i
a b
c
|
| y
1
1

=
i i i
i i i
i
a b
a d
|

y al terminar el paso n la matriz final es
n
i i

|
|
|
1 0 ... ... ...
1 0 0 0
... ... ... .... ... ...
1 0
1
2 2
1 1
..................................................................... (2.17)

Donde
1
1

=
n n n
n n n
n
a b
a d
|

Entonces el sistema de ecuaciones original se ha transformado en:

n n
x
x
x
x

|
|
|
... ...
*
1 0 ... 0 0
... ... ... ... ...
... 1 0 0
... ... 1 0
... ... ... 1
3
2
1
3
2
1
3
2
1
= ........................................... (2.18)

Y de acuerdo con la definicin de la multiplicacin de matrices se puede
establecer:

n n
x = ............................................................................................... (2.19)
80
n n n n
x x *
1 1 1
= |

y en general

1 , 1 , *
1
= =
+
n i x x
i i i i
| ................................................................ (2.20)

Este caso particular del mtodo de gauss para este tipo de matriz tridiagonal, se
conoce como el algoritmo de Thomas y se puede resumir en los siguientes pasos:

i-. Para cada ecuacin i llame b
i
el coeficiente de la incgnita i, a
i
el coeficiente de
la incgnita i-1 y c
i
el coeficiente de la incgnita i+1.
ii-. Para i=1 haga
1
1
1
b
c
= | y
1
1
1
b
d
= .............................................................................. (2.21)
Y para i=2 hasta n, haga secuencialmente:

1

=
i i i
i
i
a b
c
|
| y
1
1

=
i i i
i i i
i
a b
a d
|

......................................................... (2.22)

iii-. Para i=n haga

n n
x = ............................................................................................... (2.23)

Y para i=n-1 hasta 1, haga secuencialmente

1 , 1 , *
1
= =
+
n i x x
i i i i
| ................................................................ (2.24)

Obsrvese que en el paso ii, se calcula primero
1
y
1
, luego
2
y
2
, luego
3
y
3

y as sucesivamente hasta llegar a
n
y
n
, y en el paso iii-. se calcula primero x
n
,
luego x
n-1
y as sucesivamente hasta llegar a x
1
.
81

3-. Mtodo de Jordan Gauss.

El mtodo de Jordan Gauss es un proceso de eliminacin en dos etapas con el
fin de convertir la matriz de los coeficientes en la matriz I; la primera etapa es la
eliminacin de Gauss o sea es una eliminacin hacia delante donde se van
eliminando gradualmente en cada paso i, desde la ecuacin i+1 hasta la ecuacin
n la variable i; con este proceso se tiene al final la matriz de los coeficientes
transformada en la matriz triangular superior con todos los elementos de la
diagonal principal iguales a 1, dada por la ecuacin (2.9); la segunda etapa es un
proceso de eliminacin hacia atrs donde partiendo de la matriz triangular del final
de la primera etapa en cada paso i se elimina la incgnita x
n-i+1
de la ecuacin (n-i)
hasta la ecuacin (1). Los pasos de la segunda etapa seran entonces los
siguientes:

Paso 1
Tomando la matriz final de la primera etapa de eliminacin a cada fila i desde i=(n-
1) hasta 1, se le resta la fila n multiplicada por
i
n i
a
,
; o sea que al final de este paso
la matriz se convierte en la matriz
1 , n
Ab representada por la ecuacin (2.25) en la
cual los
1 , i
i
b estn dados por
n
n
i
n i
i
i
i
i
b a b b *
,
1 ,
= y 1 , 1 1 , 1
,
1 ,
,
= . = = n j n i a a
i
j i
i
j i



n
n
n
n
n
b o
b
b a
b a
b a
Ab
1 ... ... .. 0 . 0 0
.. . .. 0 1 0 0 0 0
... ... ... ... ... ... ...
0 ... ... .. 1 . 0 0
0 ... ... ... 1 0
0 .... .... .... 1
1 , 1
1
1 , 3
3
1 , 3
34
1 , 2
2
1 , 2
23
1 , '
1
1 , 1
12
) 1 , (

= ............................................. (2.25)


82
Paso 2.
Partiendo de la matriz
1 , n
Ab , dada por la ecuacin (2.25) se procede a eliminar
desde la fila n-2 hasta la fila 1, la incgnita x
n-1
restando a cada fila

i la fila (n-1)
multiplicada por
i
n i
a
1 ,
; o sea que al final la matriz se convierte en la ecuacin (2.26)

n
n
n
n
n
b o
b
b a
b a
b a
Ab
1 ... ... .. 0 . 0 0
.. . .. 0 1 0 0 0 0
... ... 0 ... ... ... ...
0 0 ... .. 1 . 0 0
0 0 ... ... 1 0
0 0 .... .... 1
1 , 1
1
2 , 3
3
3
34
2 , 2
2
2
23
2 , 1 '
1
1
12
) 2 , (

= ............................................... (2.26)

donde
1 , 1
1 1 ,
1 , 2 ,
*


=
n
n
i
n i
i
i
i
i
b a b b para i= (n-2) hasta 1.

Paso j

Partiendo de la matriz
1 , i n
Ab se procede a eliminar desde la fila (n-i) hasta la fila 1,
la incgnita x
n-i+1
restando a cada fila

la fila (n-i+1) multiplicada por
i
j n i
a
1 , +
; o sea
que al final la matriz se convierte en

n
n
n
n
j j n
j n
j n
j j n
i n
b
b
b
b a
b a
Ab
1 ... ... 0 .. 0 . 0 0
.. . .. 0 1 0 0 0 0
... ... 0 ... 1 ... ... ...
.... . 0 0 ... 0 .. 1 . 0 0
0 0 ... ... ....... 1 0
0 0 .... .... 1
1 , 1
1
, 1
1
2 , 2
2
2
1 , 2
2 , 1 '
1
1
12
) , (

= ................................ (2.27)

donde
1 , 1
1 1 ,
1 , ,
*

+ +

=
n
j n
i
i i
j i
i
j i
i
b a b b para i= (n-2) hasta 1.

4-. Mtodo de la Matriz Inversa.
83
Requiere obtener la inversa de la matriz de los coeficientes.

La expresin matricial para un sistema de ecuaciones es la ecuacin (2.2)

b X A = (2.2)

la cual por la definicin de inversa de la matriz A se puede transformar en

X b A X A A = =
1 1
(2.3)

De acuerdo con la ecuacin (2.3) y la definicin de multiplicacin de matrices se
puede establecer:

=
n
s
i s s i
x b a
1
1
,
* (2.4)

Con la ecuacin (2.4) y una vez obtenida la inversa de la matriz de los coeficientes
se pude obtener el valor de cada una de las incgnitas x
i
desde I=1 hasta n.

5-. Mtodo de Factorizacin de la Matriz de Coeficientes.

Este mtodo consiste en factorizar la matriz de los coeficientes A en dos matrices
L y U, tal que L sea una matriz triangular inferior y la matriz U es una matriz
triangular superior con todos los elementos de la diagonal principal iguales a uno.

El sistema de ecuaciones a resolver, como ya vimos, se puede representar
matricialmente por

b X A = o por b X A = ................................................................................. (2.2)

84
y por tanto al factorizar la matriz de los coeficientes se tendr

( ) b X LU =

y aplicando la ley asociativa del producto de matrices se tiene

( ) b X U L = * *

Definamos ahora el siguiente vector

G X U = * ..................................................................................................... (2.27)

o sea que

b G L = * ....................................................................................................... (2.28)

De acuerdo con la definicin del producto de matrices

=
=
n
s
j s s i j i
u l a
1
, , ,
* ............................................................................................ (2.29)

y la sumatoria de la ecuacin (2.29) se puede descomponer en


+ =

=
+ + =
n
j s
j s s i j j j i
j
s
j s s i j i
u l u l u l a
1
, , , ,
1
1
, , ,
*

Recordando la definicin de la matriz U se tiene que U
j,j
=1 y que U
s,j
=0 para todo
s>j , de la ecuacin anterior se tiene

85

=
=
1
1
, , , ,
j
s
j s s i j i j i
u l a l ....................................................................................... (2.30)

Por otra parte, la sumatoria de la ecuacin (2.29) tambin se puede expresar
como


+ =

=
+ + =
n
i s
j s s i j i i i
i
s
j s s i j i
u l u l u l a
1
, , , ,
1
1
, , ,
*

Recordando la definicin de la matriz L se tiene que los l
i,s
=0 para todo s>i y por
tanto

|
.
|

\
|
=

=
1
1
, , ,
,
,
*
1
i
s
j s s i j i
i i
j i
u l a
l
u ............................................................................ (2.31)

Las ecuaciones (2.30) y (2.31) permiten calcular los elementos de las matrices U y
L de la siguiente manera:

- Se calculan los elementos de la primera columna (j=1) de la matriz L
usando la ecuacin (2.30), en este caso la sumatoria vale cero. Luego se
pasa a calcular los elementos de la primera fila (i=1) de la matriz U usando
la ecuacin (2.31), en este caso la sumatoria vale cero.

- Se pasa luego a calcular los elementos de la segunda columna de la matriz
L, usando la ecuacin (2.30), en este caso la sumatoria se puede calcular
porque posee un elemento que est dado por elementos de la columna 1
de la matriz L y de la fila 1 de la matriz U, ya calculados. Se pasa luego a
calcular los elementos de la fila 2 de la matriz U usando la ecuacin (2.31),
en este caso la sumatoria se puede calcular y posee un elemento dado por
elementos de las columnas 1 de la matriz L y elementos de la fila 1 de la
matriz U, ya calculados.

- Se procede o a calcular los elementos de la columna 3 de la matriz L y
luego los de la fila 3 de la matriz U usando las ecuaciones (2.30) y (2.31)
respectivamente, y as se continua en cada paso i calculando los elementos
de la columna i de la matriz L y luego los elementos de la fila i de la matriz
86
U hasta llegar a calcular los elementos de la columna n de la matriz L y los
elementos de la fila n de la matriz U.

Una vez obtenidas las matrices L y U se procede a obtener el vector G de la
siguiente manera:

De acuerdo con la ecuacin (2.28) y recordando la definicin de multiplicacin
de matrices, ecuacin (2.29)

b G L = * .................................................................................................. (2.28)

=
=
n
j
j j i i
g l b
1
,
* ..........................................................................................
................................................................................................................

= + =
+ + =
1
1 1
, , ,
* * *
i
j
n
i j
j j i i i i j j i
g l g l g l

Y como l
i,j
=0 para j>i, la segunda sumatoria de la expresin anterior es igual a
cero y por tanto se puede escribir

(

=
1
1
,
,
* *
1
i
j
j j i i
i i
i
g l b
l
g ......................................................................... (2.32)

Con la ecuacin (2.32) se pueden ir encontrando los elementos del vector G,
los g
i
, en orden ascendente desde i=1 hasta n.

Conocidos los elementos del vector G se procede a encontrar los valores de
las x
i
, los elementos del vector X, usando la ecuacin (2.27) de acuerdo con la
cual y usando la definicin de multiplicacin de matrices se puede escribir

87

=
=
n
j
j j i i
x u g
1
,
*

= + =
+ + =
1
1 1
, , ,
* * *
i
j
n
i j
j j i i i i j j i
x u x u x u

Recordando la definicin de la matriz U, los u
i,i
=1 y los u
i,j
=0 para todo j<i, o sea
que la primera sumatoria de la expresin anterior es cero y por tanto

+ =
=
n
i j
j j i i i
x u g x
1
,
* ........................................................................ (2.33)

La ecuacin (2.33) permite obtener los valores de x
i
de manera secuencial
decreciente empezando por x
n
hasta llegar a x
1
.

Cuando la matriz de los coeficientes es una matriz tridiagonal este mtodo de
solucin del sistema de ecuaciones tambin se simplifica apreciablemente, al igual
que en el caso del procedimiento de Gauss y el mtodo resultante tambin se
conoce como algoritmo de Thomas. La matriz de los coeficientes tridiagonal se
representa por


n n
b a
c b a
c b a
c b
A
... 0 0
... ... ... ... ...
0 0
0 0 0
3 3 3
2 2 2
1 1
=

y el sistema de ecuaciones por la ecuacin (2.13)

88

n n n n
d
d
d
d
x
x
x
x
b a
c b a
c b a
c b
... ...
... 0 0
... ... ... ... ...
0 0
0 0 0
3
2
1
3
2
1
3 3 3
2 2 2
1 1
= (2.13)
Al descomponer la matriz A de los coeficientes en sus factores L y U de acuerdo
con las ecuaciones

LU A=

=
=
1
1
, , , ,
j
s
j s s i j i j i
u l a l (2.30)

|
.
|

\
|
=

=
1
1
, , ,
,
,
*
1
i
s
j s s i j i
i i
j i
u l a
l
u (2.31)

Y teniendo en cuenta la nomenclatura que se ha usado para la matriz de los
coeficientes tridiagonal, se encuentra que las matrices L y U tiene la siguiente
forma:


n n
l a
a
l a
l
L
... ...
...
3
2 2
1
= (2.34)


1
...
... ...
1
1
1
2
1

=
n
u
u
u
U (2.35)

89
Donde


1 1
b l = (2.36)


1 2 2
1
1
2 2 2
* u a b
b
c
a b l = = (2.37)

2 3 3
1 1 2 2
2
3 3 3
/ *
* u a b
b c a b
c
a b l =

= (2.38)


1 1 1
/ b c u = (2.39)


( )
1 2 2
2
1 1 2 2
2
2
/ * u a b
c
b c a b
c
u

= (2.40)


( )
2 3 3
3
1 1 2 2 2 3 3
3
3
/ * / * u a b
c
b c a b c a b
c
u

=

= (2.41)

De acuerdo con las ecuaciones (2.37) y (2.38) se puede establecer

n i u a b l
i i i i
, 2 ;
1
= =

(2.42)

De acuerdo con las ecuaciones (2.40) y (2.41)

n i
l
c
u a b
c
u
i
i
i i i
i
i
, 2 ,
1
1
1
1
= =

(2.43)

Las ecuaciones (2.36) y (2.4) permiten obtener la matriz L y la ecuacin (2.43)
permite obtener la matriz U; obsrvese adems que la diagonal inferior de la
matriz L es la misma de diagonal inferior de la matriz A.

90
Para obtener el vector G definido por la ecuacin (2.27) se parte de la ecuacin
(2.28), la cual usando la nomenclatura usada para la matriz tridiagonal se expresa
como

d G L = * (2.44)
O sea que la ecuacin (2.32) se convierte en


(

=
1
1
,
,
* *
1
i
j
j j i i
i i
i
g l d
l
g (2.45)

Y recordando la configuracin de la matriz L , ecuacin (2.34) se pueden obtener
las siguientes expresiones para calcular los componentes del vector G:


1
1
1
l
d
g = (2.47)

n i
l
g c d
g
i
i i i
i
, 2 ,
1
=

=

(2.48)


2.3-. Mtodos Indirectos (Ensayo y Error) para Resolver Sistemas de
Ecuaciones.

1-. Mtodo de Jacobi

La expresin matricial para el sistema de ecuaciones es

b X A = ......................................................................................................... (2.2)

La matriz A se puede descomponer en dos sumandos as
91

C D A = ..................................................................................................... (2.34)

donde D es una matriz diagonal cuyos elementos son todos cero con excepcin de
los elementos de la diagonal principal que son igual esa los elementos de la
diagonal principal de A; o sea

(
(
(
(
(
(

n n
n n
a
a
a
a
D
,
, 1
22
11
.... 0 0 0
0 ... 0 0
... ... ... ... ...
0 .... 0 0
0 .... 0 0
...................................................................... (2.35)

y la matriz C es una matriz cuyos elementos de la diagonal principal son cero y c
i,j

= -a
i,j
para todos los i,j; o sea

(
(
(
(
(
(

=
0 ....
.... .... .... .... ....
....
.... 0
.... 0
3 , 2 , 1 ,
, 3 2 , 3 1 , 3
, 2 3 , 2 21
, 1 3 , 1 2 , 1
n n n
n
n
n
a a a
a o a a
a a a
a a a
C
.................................................................. (2.36)

Teniendo en cuenta la ecuacin (2.34) y la propiedad distributiva del producto de
matrices se tiene:

b X C X D = * *

de donde se puede despejar

b X C X D + = * * .......................................................................................... (2.37)

92
Cuando X es la solucin del sistema de ecuaciones se cumple la ecuacin (2.37).
Pero la ecuacin (2.37) tambin nos permite calcular un posible vector solucin del
sistema de ecuaciones,
1 + m
X , a partir de un vector solucin supuesto
m
X , o sea

b X C X D
m m
+ =
+
* *
1
.................................................................................... (2.38)

Cuando se cumpla que la diferencia entre
1 + m
X y
m
X es mnima, se ha
encontrado la solucin del sistema de ecuaciones.

El mtodo de Jacobi consiste entonces en lo siguiente:

- Se supone un vector solucin
m
X y a partir de este vector solucin se
calcula un nuevo vector solucin
1 + m
X , usando la ecuacin (2.38). De
acuerdo con la definicin de las matrices D y C y del producto entre
matrices, el elemento
1 + m
i
x del vector
1 + m
X se obtiene de la siguiente
expresin al aplicar la ecuacin (2.38)

(

= + =
+
1
1 1
, ,
1
,
i
j
n
i j
m
j j i
m
j j i i
m
i i i
x a x a b x a

................................................................................................................
(

= + =
+
1
1 1
, ,
,
1
*
1
i
j
n
i j
m
j j i
m
j j i i
i i
m
i
x a x a b
a
x , n i , 1 = ....................................... (2.39)

- Con la expresin (2.39) se calculan todos los elementos del vector
1 + m
X y
una vez calculados todos los
1 + m
i
x se verifica si se cumple la siguiente
igualdad

( ) c *
1
2
1
n x x
n
i
m
i
m
i
<

=
+
............................................................................ (2.40)

93
donde c es la diferencia que se considera aceptable entre el valor supuesto
para una variable y el valor calculado. Si no se cumple la ecuacin (2.40) el
nuevo vector solucin calculado
1 + m
X se toma como un nuevo vector solucin
supuesto
1 + m
X y se procede a calcular un nuevo vector solucin
2 + m
X .

- El proceso continua hasta que se cumpla la ecuacin (2.40).
Convergencia del Mtodo Iterativo de Jacobi: Como en cualquier procedimiento
de ensayo y error, se debe demostrar que el mtodo de Jacobi es convergente, o
sea que es posible que teniendo un sistema de ecuaciones si se supone un vector
solucin para el mismo y se calcula un nuevo vector solucin aplicando la
ecuacin (2.39), despus de hacer esto durante un nmero m de iteraciones se
llegar a una situacin donde se cumpla la ecuacin (2.40). Para demostrar esto
retomemos la ecuacin (2.38)

b X C X D
m m
+ =
+
* *
1
.................................................................................... (2.38)

de la cual se puede concluir que

b D X C D X
m m
* *
1 1 1 +
+ = .................................................................... (2.41)

Cuando se tiene la solucin exacta X , la ecuacin (2.41) es

b D X C D X * *
1 1
+ = .................................................................................. (2.42)

94
Restando la ecuacin (2.42) de la (2.41) y definiendo M como D
-1
C, se tiene
.....................................................................................................................
( ) ( ) X X M X X
m m
=
+
*
1
............................................................................ (2.43)

En la ecuacin (2.43), ( ) X X
m

+1
es el error que se comete en la iteracin (m+1) y
se representa por
m+1
, y ( ) X X
m
es el error que se comete en la iteracin m y se
representa por
m
; o sea que la ecuacin tambin se puede escribir como
.....................................................................................................................
2 1 1
* * * * * *
+
= = =
m m m m
M M M M M M c c c c ............................................

Finalmente, de acuerdo con la expresin anterior se puede escribir

0 1
c c
m m
M =
+
................................................................................................. (2.44)
En la ecuacin (2.44), el trmino
0
es el error que se comete al hacer la
suposicin inicial para iniciar el proceso de solucin por el mtodo de ensayo y
error de Jacobi, el cual ser convergente si despus de m iteraciones, cuando m
es muy grande, el error entre la solucin calculada y la solucin real,
m+1
, tiende a
cero; o sea que usando la ecuacin (2.44) el mtodo de Jacobi ser convergente
si


= =
+
m
M
m m
0 lim lim
0 1
c c
...................................................................... (2.45)

Puesto que
0
no es cero, ya que lo lgico es que el vector solucin inicial
supuesto es diferente del vector solucin real, la nica posibilidad de que se
cumpla la ecuacin (2.45) es que la matriz M multiplicada m veces por si misma se
convierta en la matriz 0, o sea una matriz en la que todos sus elementos sean
cero. Una matriz cumpla con esta condicin en los siguientes casos:

95
1-.)Una condicin suficiente para que

=
m
M
m
0 lim
es que 1 < M donde M se
conoce como la norma de la matriz M.

2-.)
( )

< =
m
M M
m
0 0 lim

donde (M) es el radio espectral de M y est dado por ( ) { } max ) ( = M , donde
son los valores propios de la matriz M.

2-. Mtodo de Gauss Seidel

El mtodo de Gauss Seidel es una mejora del mtodo de Jacobi para acelerar
su convergencia, en donde las valores que se van calculando en una iteracin
dada para las variables se van usando en la misma iteracin para calcular los
valores de las variables siguientes. Se basa en lo siguiente:

La matriz de los coeficientes A de la ecuacin (2.2) para el sistema de ecuaciones
se descompone en tres sumandos, las matrices D, U y L

U L D A = ............................................................................................... (2.46)

donde: D est definida otra vez por la ecuacin (2.35)

(
(
(
(
(
(

n n
n n
a
a
a
a
D
,
, 1
22
11
.... 0 0 0
0 ... 0 0
... ... ... ... ...
0 .... 0 0
0 .... 0 0
...................................................................... (2.35)

L est definida por

96
(
(
(
(
(
(

=
0 ....
.... .... .... .... ....
0 .... 0
0 .... 0 0
0 .... 0 0 0
3 , 2 , 1 ,
2 , 3 1 , 3
21
n n n
a a a
a a
a
L
....................................................................... (2.47)

y U est definida por

.....................................................................................................................
(
(
(
(
(
(

=
0 .... 0 0 0
.... .... .... .... ....
.... 0 0
.... 0 0
.... 0
, 3
, 2 3 , 2
, 1 3 , 1 2 , 1
n
n
n
a o
a a
a a a
U
...................................................................... (2.48)

El sistema de ecuaciones se puede entonces representar matricialmente por

( ) b X U L D = *

o sea
( ) b X U X L D + = * * .................................................................................. (2.49)

La ecuacin (2.49) cumple si X es la solucin exacta, pero al igual que en el caso
de Jacobi tambin se puede plantear un procedimiento de ensayo y error de tal
forma que

( ) b X U X L D
m m
+ =
+
* *
1
............................................................................ (2.50)

donde
m
X es un vector solucin supuesto en una iteracin y
1 + m
X es un nuevo y
mejor solucin calculado en la misma iteracin.

97
De acuerdo con las definiciones de las matrices D, L y U y de la suma y
multiplicacin de matrices, los elementos
1 + m
i
x del vector
1 + m
X se pueden calcular
de la siguiente manera

m
n n
m m m
x a x a x a b x a
, 1 3 3 , 1 2 2 , 1 1
1
1 1 , 1
.. .......... =
+


( )
m
n n
m m m
x a x a x a b
a
x
, 1 3 3 , 1 2 2 , 1 1
1 , 1
1
1
.. ..........
1
=
+


m
n n
m m m m
x a x a x a b x a x a
, 2 4 4 , 2 3 3 , 2 2
1
2 2 , 2
1
1 1 , 2
.. .......... = +
+ +


( )
m
n n
m m m
x a x a x a b
a
x
, 2 3 3 , 2
1
2 1 , 2 2
2 , 2
1
2
.. ..........
1
=
+ +


y en general


+ =

=
+ +
= +
n
i j
m
j j i i
i
j
m
i i i
m
j j i
x a b x a x a
1
,
1
1
1
,
1
,

(

= + =
+ +
1
1 1
,
1
,
,
1
*
1
i
j
n
i j
m
j j i
m
j j i i
i i
m
i
x a x a b
a
x ...................................................... (2.51)

La ecuacin (2.51) es la ecuacin general para calcular el nuevo vector solucin
en la iteracin (m+1) del proceso de Gauss Seidel; obsrvese que los elementos
se van calculando secuencialmente y cuando se est calculando el elemento
1 + m
i
x
del vector
1 + m
X , usando la ecuacin (2.51), se usan los valores de los elementos
1 + m
j
x para todo j<i los cuales ya han sido calculados en pasos anteriores de la
misma iteracin, los elementos x
j
para j>i son los calculados en la iteracin m
porque an no han sido calculados en la iteracin (m+1).

98
La convergencia de este mtodo se prueba de la misma forma que en el mtodo
de Jacobi teniendo en cuenta que en este caso la matriz M est dada por

( ) U L D M
1
= ............................................................................................. (2.52)

y para que haya convergencia la matriz M dada por la ecuacin (2.52) debe
cumplir con una de las dos condiciones que se plantearon para la matriz M del
mtodo de Jacobi.

3-. Mtodo de Relacin Sucesiva Puntual (PSOR)

Es una mejora del mtodo de Gauss Seidel para acelerar la convergencia, como
se puede ver aplicando el siguiente desarrollo

Aplicando la ecuacin (2.46) a la ecuacin (2.2) se tiene

( ) b X U L D = *

Si se introduce un factor W tal que
( ) b w X U L D w = *

La ecuacin anterior se puede modificar de la siguiente forma

( ) X D wb X wU wL wD D + = + *

( ) ( ) | | X wU w D b w X wL D * 1 * * + + = ....................................................... (2.53)

Nuevamente, la ecuacin (2.49) cumple si X es la solucin exacta del sistema de
ecuaciones, pero tambin puede plantearse a partir de ella un procedimiento de
ensayo y error escribindola de la siguiente forma:

99
.....................................................................................................................
( ) ( ) | |
m m
X wU w D b w X wL D * 1 * *
1
+ + =
+
................................................. (2.54)

La ecuacin (2.54) plantea que es posible suponer un vector solucin
m
X para el
sistema de ecuaciones, que no es la solucin exacta, y a partir de el calcular un
nuevo vector
1 + m
X que es mejor solucin. Esto es, al igual que en los casos de
Jacobi y Gauss Seidel, un proceso de ensayo y error y en una iteracin m se
supone un vector solucin
m
X y se calcula un nuevo vector solucin
1 + m
X que es
mejor solucin que el vector solucin supuesto en la misma iteracin; para una
iteracin m el vector supuesto es el calculado en la iteracin anterior. El proceso
termina cuando en una iteracin se cumpla que la diferencia entre el vector
supuesto y el calculado est dentro de una tolerancia establecida, o sea cuando
se cumpla la ecuacin (2.42)

Los elementos del vector
1 + m
X se calculan de la siguiente manera:

Nuevamente recordando la definicin de las matrices D, L, U y de las operaciones
de suma y multiplicacin de matrices se tiene aplicando la ecuacin (2.54)


+ =

=
+ +
+ + = +
n
i j
m
j j i
m
i i i i
i
j
m
j i i
m
j j i
x a w x a w b w x a x a w
1
, ,
1
1
1
,
1
,
* * ) 1 ( * *

y despejando
1 + m
i
x se tiene

( )
(

+ =

= + =
+ +
1
1
1
, ,
,
1
* 1
i
i j
n
i j
m
j j i
m
j j i i
i i
m
i
m
i
x a x a b
a
w
x w x ........................................ (2.55)


100
Observando la ecuacin (2.55) se ve que el segundo sumando es la expresin
para calcular
1 + m
i
x por el mtodo de Gauss Seidel, ecuacin (2.51), multiplicado
por w, o sea que el valor de
1 + m
i
x en este mtodo se obtiene haciendo un promedio
ponderado entre su valor obtenido entre la iteracin anterior y su valor en la
iteracin (m+1) obtenido por el mtodo de Gauss - Seidel.

El trmino w se conoce como factor de relajacin y dependiendo de su valor se
tendr mejor convergencia o no que con el mtodo de Gaus Seidel, y el mtodo
se conoce como de relajacin; cuando el valor de w es menor que 1 se dice que
se tiene subrelajacin, si es mayor que 1 se habla de sobrerelajacin y si es igual
a 1 se tiene el mtodo de Gauss Seidel.

La convergencia de este mtodo tambin se prueba siguiendo el procedimiento
aplicado en el mtodo de Jacobi y partiendo de la ecuacin (50), sea que en este
caso la matriz M est dada por

( ) ( ) | | wU D w wL D M + =

1 *
1
.................................................................... (2.56)

La convergencia se tendr si la matriz M cumple con una de las dos condiciones
que se plantearon en la demostracin de convergencia del mtodo de Jacobi.

2.4- APROXIMACION DE ECUACIONES DIFERENCIALES A DIFERENCIAS
FINITAS.

En el captulo anterior se discutieron algunas ecuaciones diferenciales parciales
denominadas ecuaciones fundamentales de flujo en medios porosos. Las
soluciones a estas ecuaciones permiten simular el comportamiento de un
yacimiento, bajo diferentes condiciones de flujo, en trminos de la posicin y el
tiempo. En algunos casos de limitada aplicacin, estas ecuaciones pueden ser
resueltas analticamente. Sin embargo, en la mayora de las aplicaciones se
requiere utilizar tcnicas numricas para su solucin.

Este captulo se dedica a la discusin de los conceptos bsicos que rigen la
aproximacin de una ecuacin diferencial parcial a diferencias finitas.
Inicialmente, se discuten algunas nociones sobre variables discretas y diferencias
101
finitas. Luego, se consideran las aproximaciones para la primera y segunda
derivada, as como los esquemas de aproximacin implcita y explcita.
Posteriormente, se discuten los conceptos de error de truncamiento, estabilidad,
convergencia y consistencia. En la parte final del captulo, se habla de sistemas
de malla y condiciones de lmite. A travs de todo el captulo se trabaja con una
ecuacin de flujo simple con la finalidad de introducir los conceptos bsicos
referentes a los temas antes mencionados. Esto facilita considerablemente la
aplicacin y extensin de conceptos bsicos a ecuaciones de flujo complejas que
gobiernan el comportamiento real de un yacimiento.

2.4.1-.DISCRETIZACION, VARIABLES DISCRETAS Y DIFERENCIAS FINITAS
Considrese un medio poroso lineal a travs del cual fluye un fluido compresible o
levemente compresible y el cual se encuentra a una presin inicial
i
p (Figura 2.1).
Si el sistema se abre a producir en el punto L x = , la presin disminuye a medida
que se incrementa el volumen de fluido extrado del sistema. Supngase que para
este sistema, la variacin de la presin p , en funcin de la posicin x y el tiempo
t , est descrita por la siguiente ecuacin diferencial:

t
p
x
p
c
c
=
c
c
2
2
......................................................................................... (2.57)

La solucin analtica a la Ecuacin 2.57 permite obtener una funcin continua
( ) t x p p , = para calcular la presin en cada punto a un tiempo determinado. Un
grfico de p en funcin de x y t lucira en forma anloga a como se ilustra en la
Figura 2.1.

Para la solucin numrica de una ecuacin diferencial de flujo, tal como la
Ecuacin 2.57, se requiere tratar el yacimiento como un conjunto de bloques. Esta
solucin permite calcular la presin en cada bloque a diferentes intervalos de
tiempo. A diferencia de la solucin analtica, en lugar de obtenerse una expresin
que relacione la presin como una funcin continua de la posicin y el tiempo, se
obtiene una serie de valores discretos de presin, cada uno de los cuales
corresponde a un bloque determinado del yacimiento. La accin de dividir el
yacimiento en un nmero determinado de bloques e intervalos de tiempo se
denomina discretizacin. La Figura 2.2 presenta el ejemplo hipottico de un
yacimiento real y su correspondiente discretizacin. A las longitudes de cada
bloque,
1
x A ,
2
x A , ... , y a cada intervalo de tiempo,
1
t A ,
2
t A , ... , se les denomina
diferencias finitas. La Figura 2.3 presenta la forma tpica de la distribucin de
presin en un sistema lineal a diferentes tiempos obtenida mediante la solucin
numrica de una ecuacin fundamental de flujo. A manera de ejemplo, se pudiera
considerar que las Figuras 2.1 y 2.3 corresponden a las soluciones analtica y
numrica, respectivamente, de una misma ecuacin de flujo (por ejemplo de la
Ecuacin 2.57) para un mismo sistema.
102

Figura 2.1-.Forma Tpica de la Distribucin de Presin en Sistema Lineal
(Solucin Analtica)


Los siguientes son algunos comentarios importantes relacionados con la
discretizacin de un yacimiento:

a. Los lados de cada bloque, perpendiculares a la direccin de flujo, actan como
barreras de espesor infinitesimal (vase Figura 2.2).

b. La rata de entrada y salida de fluidos en cada bloque est gobernada por las
permeabilidades de las barreras adyacentes a cada bloque y la diferencia de
presin entre ellos.

c. Las propiedades dentro de cada bloque son las mismas en todos los puntos del
bloque. Por ejemplo: a un tiempo determinado, a cada bloque le corresponde
un nico valor de presin y de permeabilidad relativa.

d. La variacin de las propiedades en el yacimiento se representa mediante la
variacin de las propiedades en cada bloque. Por ejemplo: a cada bloque se le
puede asignar un valor diferente de porosidad. Esto implica que pueden existir
variaciones fuertes en las propiedades de bloques adyacentes. El grado de
variacin de las propiedades de diferentes bloques es funcin, en parte, del
tamao de cada bloque.
x x = L
x = 0 x = L x
t
1
t
2
t
3
t
4
p
x = 0
103
Figura 2.2-. Esquema de un Sistema Real y su Discretizacin


e. La precisin con la cual el comportamiento real del yacimiento puede ser
descrito por el modelo depende del nmero de bloques utilizados. En general,
a mayor nmero de bloques mayor precisin. Sin embargo, el nmero de
bloques est limitado por factores tales como el costo del simulador, en lo que
se refiere a eficiencia computacional y capacidad de almacenamiento, y la
disponibilidad de datos de entrada.

f. La vida del yacimiento es discretizada en intervalos de tiempo. Las condiciones
del yacimiento se definen nicamente al principio y al final de cada intervalo de
tiempo, razn por la cual las condiciones en cada bloque pueden cambiar
considerablemente de un intervalo de tiempo a otro. Intervalos de corta
duracin (es decir, alto nmero de intervalos) hacen que estos cambios sean
menos fuertes y aumentan la precisin de los resultados.

g. La discretizacin del yacimiento hace que el comportamiento continuo de las
condiciones en l sea distorsionado. Por ejemplo: la Figura 3.4 presenta un
esquema de la distorsin ocasionada en la saturacin de agua en un punto fijo
de un sistema lineal debido a su discretizacin.



Ax
1
Ax
2
Ax
3 Ax
4
Barrera Impermeable
k
1
p
1
k
2
p
2
k
3
p
3
k
4
p
4
Bloque 1 Bloque 3 Bloque 4 Bloque 2
|
1
|
2
|
3
|
4
104
Figura 2.3-.Forma Tpica de la Distribucin de Presin en un Sistema Lineal
(Solucin Numrica).

2.4.2-. APROXIMACIONES PARA LA PRIMERA DERIVADA
Probablemente, el mtodo ms comnmente utilizado para dar solucin numrica
a la Ecuacin 2.57 es expresar la funcin ( ) t , x p mediante una serie de Taylor.

Considrese una funcin ( ) x f continua. La funcin puede ser expandida
alrededor de un punto a , localizado a una distancia a x del punto x (Figura
2.5), mediante una serie de Taylor (Protter y Morrey, 1980):

( )

=
0
) (
!
) (
) (
n
n n
n
a x a f
x f ............................................................................ (2.58)

En la Ecuacin 2,58,
( )
( ) a f
n
es la n-sima derivada de f evaluada en el punto
a x = . Es decir, la Ecuacin 3.2 puede ser escrita como:



x = 0 x = L x
x = 0 x = L x
p
1
t
2
t
3
t
4
t
105



Figura 2.4-. Distorsin Ocasionada en la Saturacin de Agua en un Punto
Determinado del Yacimiento Debido a la Discretizacin del Tiempo


Figura 2.5-. Esquema Ilustrativo de los Bloques Tenidos en Cuenta en la
Aproximacin Progresiva de la Primera Derivada


x = 0
x = L
x
i x
i + 1
a x
x a
x
Tiempo
: Yacimiento
: Simulador
0
1
S
w
106
( ) ( )
( ) ( ) ( ) ( )
! 3 ! 2 ! 1 ! 0
3
3
3 2
2
2
0
a x
x
f a x
x
f a x
x
f a x
a f x f
a x a x a x

|
|
.
|
c
c
+

|
|
.
|
c
c
+

|
|
.
|
c
c
+

=
= = =

( )
+

|
|
.
|
c
c
+ +
=
! n
a x
x
f
n
a x
n
n
........................................................ (2.59)

2.4.3-. La Aproximacin Progresiva y su Error de Truncamiento. Supngase
un sistema lineal de bloques de igual longitud, tal como se ilustra en la Figura 2.5.
A cada bloque le corresponde una variable discreta
1
x ,
2
x , ... ,
1 i
x ,
i
x ,
1 + i
x , ...,
n
x , de acuerdo a su posicin dentro del sistema. El incremento entre dos puntos
sucesivos,
i
x y
1 + i
x , est dado por la diferencia finita x A . Si se definen los
puntos
i
x a = y
1 +
=
i
x x , entonces
i i
x x a x =
+ 1
. Por ser bloques de igual
longitud, x a x A = .

Si adicionalmente se tiene en cuenta que la funcin ( ) x f que se pretende
expandir mediante la serie de Taylor es la presin en el sistema, P , entonces de
la Ecuacin 2,59 se tiene:
( ) ( )
( ) ( ) ( )
+ +
A

|
|
.
|
c
c
+
A

|
|
.
|
c
c
+
A

|
|
.
|
c
c
+ =
+

! 3 ! 2 ! 1
3
3
3 2
2
2
1
3 2 1
x
x
P x
x
P x
x
P
x P x P
x x x
i i


( )
+
A

|
|
.
|
c
c
! n
x
x
P
n
x
n
n
n
....................................................................... (2.60)
En la Ecuacin 2.60, ( )
1 + i
x P y ( )
i
x P son las presiones en los bloques 1 + i e i ,
respectivamente. Para simplificar la notacin, las presiones ( )
1 + i
x P y ( )
i
x P se
notarn en adelante como
1 + i
P y
i
P , respectivamente. Similarmente, la derivada
xi
x
P
|
|
.
|
c
c
se notar como
i
x
P
|
|
.
|
c
c
. Teniendo en cuenta esta notacin, la Ecuacin 2.60
puede ser escrita como:

( ) ( ) ( ) ( )
+
A

|
|
.
|
c
c
+ +
A

|
|
.
|
c
c
+
A

|
|
.
|
c
c
+
A

|
|
.
|
c
c
+ =
+
! ! 3 ! 2 ! 1
3
3
3 2
2
2
1
n
x
x
P x
x
P x
x
P x
x
P
P P
n
i
n
n
i i i
i
i
..................................................................................................................... (2.61)

Resolviendo para
i
x
P
|
|
.
|
c
c
:

( ) ( ) ( )
+
A

|
|
.
|
c
c

A

|
|
.
|
c
c

|
|
.
|
c
c

=
|
|
.
|
c
c

+
! ! 3 ! 2
1 2
3
3
2
2
1
n
x
x
P x
x
P x
x
P
x
P P
x
P
n
i
n
n
i i
i i
i
....... (2.62)
107

La Ecuacin 2.62 puede ser escrita como:

p
i
i i
i
R
x
P P
x
P
+
A

=
|
|
.
|
c
c
+ 1
................................................................................. (2.63)

En la Ecuacin 2.63, el trmino
p
i
R est dado por la siguiente ecuacin:

( ) ( )
(
(

+
A

|
|
.
|
c
c
+
A

|
|
.
|
c
c
=
! 3 ! 2
2
3
3
2
2
x
x
P x
x
P
R
i i
p
i
.................................................. (2.64)

La Ecuacin 2.63 puede ser escrita en forma aproximada de la siguiente forma:

x
P P
x
P
i i
i
A

~
|
|
.
|
c
c
+ 1
......................................................................................... (2.65)

A la Ecuacin 2.65 se le conoce como la aproximacin progresiva a la primera
derivada y permite evaluar la primera derivada de la presin en el bloque i en
trminos de la presiones en los bloques i e 1 + i .

A la parte truncada de la Ecuacin 2.62,
p
i
R , se le conoce como error de
truncamiento de la aproximacin progresiva. En este caso, se dice que el error de
truncamiento es de primer orden debido a que la menor potencia a la cual se
encuentra elevada la diferencia finita x A , en el error de truncamiento, es uno
(Ecuacin 2.64). En general, se dice que la aproximacin numrica de una
ecuacin diferencial tiene un error de truncamiento de orden n , lo cual se nota
como ( ) | |
n
x O A , cuando la mnima potencia de la diferencia finita, x A , en la parte
truncada es n .

Obsrvese que a mayor orden del error de truncamiento, mejor es la aproximacin
numrica de la ecuacin diferencial. Por ejemplo, considrese dos esquemas de
aproximacin cuyos errores de truncamiento sean:

|
( ) ( ) ( )
(
(

+
A

|
|
.
|
c
c
+
A

|
|
.
|
c
c
+
A

|
|
.
|
c
c
=
! 4 ! 3 ! 2
3
4
4 2
3
3
2
2
1
x
x
P x
x
P x
x
P
R
i i i
i


y

108
|
( ) ( )
(
(

+
A

|
|
.
|
c
c
+
A

|
|
.
|
c
c
=
! 4 ! 3
3
4
4 2
3
3
2
x
x
P x
x
P
R
i i
i


La aproximacin correspondiente al error de truncamiento |
2
i
R , es la de mayor
exactitud debido a que se est truncando nicamente desde el trmino ( )
2
x A , en
tanto que la aproximacin correspondiente al error de truncamiento |
1
i
R est
truncando desde el trmino ( ) x A .

2.4.4-. La Aproximacin Regresiva y su Error de Truncamiento. Considere el
sistema lineal de bloques ilustrado en la Figura 2.6. Si se definen los puntos
i
x a = y
1
=
i
x x , entonces
i i
x x a x =
1
( )
1
=
i i
x x x A = . Teniendo en
cuenta estas definiciones y recordando que ( ) P x f = , de la Ecuacin 2.59, se
tiene:

( ) ( ) ( )
( ) ( )
+ +
A

|
|
.
|
c
c
+
A

|
|
.
|
c
c
+ A
|
|
.
|
c
c
+ =

! 3 ! 2
3
3
3 2
2
2
1
x
x
P x
x
P
x
x
P
x P x P
Xi Xi Xi
i i


( )
+
A

|
|
.
|
c
c
! n
x
x
P
n
Xi
n
n
................................................................. (2.66)

Con la finalidad de simplificar notacin, ( )
1 i
x P se notar como
1 i
P . Teniendo en
cuenta est notacin, y la notacin para ( )
i
x P definida anteriormente, la Ecuacin
2.66 puede ser escrita como:

( ) ( ) ( ) ( )
+

|
|
.
|
c
c
+ +
|
|
.
|
c
c

|
|
.
|
c
c
+ |
.
|
c
c
=

! 3

! 2

! 1

3
3
3 2
2
2
1
n
x
x
P x
x
P x
x
P x
x
P
P P
n
i
n
n
i i
i
i
i
(2.67)

Resolviendo la Ecuacin 2.66 para
i
x
P
|
|
.
|
c
c
, se tiene:

( ) ( ) ( )
+
A

|
|
.
|
c
c
+ +
A

|
|
.
|
c
c

|
|
.
|
c
c
+
A

=
|
|
.
|
c
c

! ! 3 ! 2
1 2
3
3
2
2
1
n
x
x
P x
x
P x
x
P
x
P P
x
P
n
i
n
n
i i
i i
i
.. (2.68)

La Ecuacin 2.68 puede ser escrita como:

r
i
i i
i
R
x
P P
x
P
+
A

=
|
|
.
|
c
c
1
................................................................................. (2.69)

109


Figura 2.6-. Esquema Ilustrativo de los Bloques Tenidos en Cuenta en la
Aproximacin Regresiva de la Primera Derivada



En la Ecuacin 2.69
r
i
R est dado por:

( ) ( ) ( )
+
A

|
|
.
|
c
c
+ +
A

|
|
.
|
c
c

|
|
.
|
c
c
=

! ! 3 ! 2
1 2
3
3
2
2
n
x
x
P x
x
P x
x
P
R
n
i
n
n
i i
r
i
.................... (2.70)
La Ecuacin 2.69 puede ser aproximada de la siguiente forma:

x
P P
x
P
i i
i
A

~
|
|
.
|
c
c
1
......................................................................................... (2.71)

A la Ecuacin 2.71 se le conoce como la aproximacin regresiva de la primera
derivada. Su error de truncamiento est dado por la Ecuacin 2.70. Obsrvese
que el error de truncamiento de la aproximacin regresiva es de primer orden,
( ) | | x O A .

2.4.5-. La Aproximacin Central y su Error de Truncamiento. Si se sustrae la
Ecuacin 2.67 de la Ecuacin 2.61, se obtiene:

x = 0 x = L
a x = - ( x a)
x
i
x
i - 1
x a
110
( ) ( ) ( )
( ) ( ) ( )
(
(

+
A

|
|
.
|
c
c

|
|
.
|
c
c
+
A

|
|
.
|
c
c

(
(

+
A

|
|
.
|
c
c
+
A

|
|
.
|
c
c
+
A

|
|
.
|
c
c
+ =
+

! 3 ! 2 ! 1
! 3 ! 2 ! 1
3
3
3 2
2
2
3
3
3 2
2
2
1 1
x
x
P x
x
P x
x
P
P
x
x
P x
x
P x
x
P
P P P
i i i
i
i i i
i i i


Simplificando,

( ) ( ) ( ) ( )

! ! ! 3
2
! 1
2
3
3
3
1 1
n
x
x
P
n
x
x
P x
x
P x
x
P
P P
n
i
n
n n
i
n
n
i i
i i
A

|
|
.
|
c
c
+
A

|
|
.
|
c
c
+
A

|
|
.
|
c
c
+
A

|
|
.
|
c
c
=
+


Despejando
i
x
P
|
|
.
|
c
c
,

( )
( ) ( ) ( )
+
A

|
|
.
|
c
c

|
|
.
|
c
c

A

|
|
.
|
c
c

=
|
|
.
|
c
c

+
! 2 ! 2 ! 3 2
1 1 2
3
3
1 1
n
x
x
P
n
x
x
P x
x
P
x
P P
x
P
n
i
n
n n
i
n
n
i
i i
i
..................................................................................................................... (2.72)

La Ecuacin 2.72 puede ser escrita como:

( )
c
i
i i
i
R
x
P P
x
P
+
A

=
|
|
.
|
c
c
+
2
1 1
................................................................................ (2.73)

En la Ecuacin 2.73,
c
i
R est dada por:
( ) ( ) ( )

A

|
|
.
|
c
c

|
|
.
|
c
c
+
A

|
|
.
|
c
c
=

! ! ! 3
1 1 2
3
3
n
x
x
P
n
x
x
P x
x
P
R
n
i
n
n n
i
n
n
i
c
i
............ (2.74)

El error el truncamiento,
c
i
R , tiende a cero cuando x A tiende a cero. En este
caso, la Ecuacin 2.73, se simplifica de la siguiente forma:

( ) x
P P
x
P
i i
i
A

~
|
|
.
|
c
c
+
2
1 1
...................................................................................... (2.75)

A la Ecuacin 2.75 se le conoce como la aproximacin central a la primera
derivada. Su error de truncamiento, dado por la Ecuacin 2.74, es de segundo
orden, ( ) | |
2
x O A . Obsrvese que el error de truncamiento de la aproximacin
central es menor que el error de truncamiento de las aproximaciones progresiva y
regresiva.
111

Las Figuras 2.7, 2.8 y 2.9 presentan la interpretacin geomtrica de las
aproximaciones progresiva, regresiva y central, respectivamente. De estas
figuras, se observa que la aproximacin progresiva estima la primera derivada de
la presin en el bloque i con base a los valores de presin en el bloques i y en el
bloque posterior ( 1 + i ). La aproximacin regresiva estima la primera derivada de la
presin en el bloque i con base a los valores de presin en el bloque i y en el
bloque anterior ( 1 i ). La aproximacin central calcula la primera derivada de la
presin en el bloque i con base a los valores de presin en los bloques anterior
( 1 i ) y posterior ( 1 + i ). Este ltimo hecho es justamente la razn por la cual el
error de truncamiento de la aproximacin central es menor que los errores de
truncamiento de las aproximaciones progresiva y regresiva.


2.4.6 -.APROXIMACIN NUMRICA PARA LA SEGUNDA DERIVADA


La aproximacin numrica para la segunda derivada puede ser obtenida sumando
las Ecuaciones 2.61 y 2.67:

( )
( ) ( )
+
A

|
|
.
|
c
c
+
A

|
|
.
|
c
c
+ A
|
|
.
|
c
c
+ = +
+
! 6
2
! 4
2 2
6
6
6 4
4
4
2
2
2
1 1
x
x
P x
x
P
x
x
P
P P P
i i i
i i i


Despejando
i
x
P
|
|
.
|
c
c
2
2
,


Figura 2.7-. Interpretacin Geomtrica de la Aproximacin Progresiva
Ax
P
i + 1
P
i
P
i + 1
P
i
progresiva
x
p
(

c
c
:
verdadera
x
p
(

c
c
:
x

x
i + 1
112



Figura 2.8-. Interpretacin Geomtrica de la Aproximacin Regresiva

Figura 2.9-. Interpretacin Geomtrica de la Aproximacin Central


( )
( ) ( )
+
A

|
|
.
|
c
c

A

|
|
.
|
c
c

A
+
=
|
|
.
|
c
c
+
! 6
2
! 4
2
2
4
6
6 2
4
4
2
1 1
2
2
x
x
P x
x
P
x
P P P
x
P
i i
i i i
i
........... (2.76)
P
i 1
P
i
P
i + 1
x
i
x
i 1
x
i + 1
verdadera
x
p
(

c
c
:
centrada
x
p
(

c
c
:
Ax
P
i 1
P
i
P
i
P
i 1
x
i
x
i 1
regresiva
x
p
(

c
c
:
verdadera
x
p
(

c
c
:
113

La Ecuacin 2.76 puede ser escrita de la siguiente forma:

( )
2
2
1 1
2
2
2
i
i i i
i
R
x
P P P
x
P
+
A
+
=
|
|
.
|
c
c
+
.................................................................. (2.77)

En la Ecuacin 2.77,
2
i
R est dada por:

( ) ( )
(
(

+
A

|
|
.
|
c
c
+
A

|
|
.
|
c
c
=
! 6
2
! 4
2
4
6
6 2
4
4
2
x
x
P x
x
P
R
i i
i
.......................................... (2.78)

2
i
R tiende a cero cuando x A tiende a cero. En este caso, de la Ecuacin 2.77 se
obtiene:

( )
2
1 1
2
2
2
x
P P P
x
P
i i i
i
A
+
~
|
|
.
|
c
c
+
.......................................................................... (2.79)

A la Ecuacin 2.79 se le conoce como la aproximacin numrica para la segunda
derivada. Esta expresin permite estimar el valor de
2
2
x
P
c
c
en el bloque i en
trminos de las presiones en los bloques 1 i , i e 1 + i . El error de truncamiento
de esta ecuacin est dado por la Ecuacin 2.78 de donde se puede observar que
el error de truncamiento para la segunda derivada es de segundo orden.

Anlogamente se pueden plantear expresiones para P
i+2
y P
i-2
de la siguiente
manera

( ) ( ) ( ) ( )
..
!
2
..
! 3
2
! 2
2
! 1
2
3
3
3 2
2
2
2
+
A

|
|
.
|
c
c
+ +
A

|
|
.
|
c
c
+
A

|
|
.
|
c
c
+
A

|
|
.
|
c
c
+ =
+
n
x
x
P x
x
P x
x
P x
x
P
P P
n
i
n
n
i i i
i
i

..................................................................................................................... (2.80)

( ) ( ) ( ) ( )
+
A

|
|
.
|
c
c
+ +
A

|
|
.
|
c
c

|
|
.
|
c
c
+
A

|
|
.
|
c
c
=

!
2
! 3 ! 2
2
! 1
3
3
3 2
2
2
2
n
x
x
P x
x
P x
x
P x
x
P
P P
n
i
n
n
i i i
i
i

..................................................................................................................... (2.81)

Si la ecuacin (2.61) la multiplicamos por 4 y le restamos la ecuacin (2.81) se
puede tener la siguiente expresin para
1
|
.
|

\
|
x
P
o
o


114
( ) ( )
.....
! 4 * 2
12
! 3 * 2
4
2
3 4
4
4 3
3
3 2
2 1
+
|
|
.
|

\
| A

|
|
.
|

\
| A

A

= |
.
|

\
|
+ +
i i
i i i
i
x
P x
x
P x
x
P P P
x
P
o
o
o
o
o
o


La expresin anterior en forma resumida se puede escribir como

( ) ( )
2
2 1
2
3 4
x
x
P P P
x
P
i i i
i
A +
A

= |
.
|

\
|
+ +
u
o
o


y si se aproxima la expresin anterior a

x
P P P
x
P
i i i
i
A

~ |
.
|

\
|
+ +
2
3 4
2 1
o
o
...................................................................... (2.82)

se est cometiendo un error del orden de ( )
2
x A .

Combinando las ecuaciones (2.61) y (2.79)-(2.81) se puede obtener

( ) ( )
2
2 1 1 2
12
8 8
x
x
P P P P
x
P
i i i i
i
A +
A
+ +
= |
.
|

\
|
+ +
u
o
o
................................................... (2.83)

la cual se puede aproximar a

x
P P P P
x
P
i i i i
i
A
+ +
~ |
.
|

\
|
+ +
12
8 8
2 1 1 2
o
o
.................................................................. (2.84)

cometiendo un error del orden de( )
2
x A , y

( )
( ) ( )
4
2
2 1 1 2
2
2
12
16 30 16
x
x
P P P P P
x
P
i i i i i
i
A +
A
+ +
=
|
|
.
|

\
|
+ +
u
o
o


la cual si se aproxima a

( )
2
2 1 1 2
2
2
12
16 30 16
x
P P P P P
x
P
i i i i i
i
A
+ +
~
|
|
.
|

\
|
+ +
o
o
.................................................. (2.85)

se comete un error del orden de ( )
4
x A .

En general se puede decir que es posible obtener expresiones para la primera y
segunda derivada en funcin del valor de la funcin en diferentes puntos; por
ejemplo supongamos se quiere obtener una expresin para la segunda derivada
115
en funcin del valor de la funcin en los puntos i-1, i e i+1; o sea se quiere tener
una expresin para
i
x
P
|
|
.
|

\
|
2
2
o
o
en funcin de P
i-1
, P
i
y P
i+1
de tal forma que

u o o o + + + =
+ 1 3 2 1 1
"
i i i i
P P P P

Para hallar los
i
se usa el mtodo de los coeficientes indeterminados de la
siguiente manera:

- Se reemplazan P
i-1
, P
i
y P
i+1
por su serie y al reemplazarlos en la ecuacin para
"
i
P se tiene

( ) ( )( ) ( )
( )
( )( ) ( ) ! 3 /
! 2
3
3 1
"
2
3 1
'
3 1 3 2 1
"
x P
x
P x P P
i i i i
A +
A
+ + A + + + = o o u o o o o o o o

- Los coeficientes de la izquierda se hacen iguales a los de la derecha para cada
uno de los trminos similares y entonces

0
3 2 1
= + + o o o

( ) 3 0
1 3 1
o o o o = = A x
( )( )
( ) ( )
2
1
2
1
2
3 1
1 2
2 1 ! 2 /
3
x x
x
A
= =
A
= = A + o o o o o
( ) ( )
2
2 2
2
3 2 1
2
0
2
x x A
= = +
A
= + + o o o o o

- Reemplazando ahora los valores de
1
,
2
y
3
en la ecuacin para
"
i
P se tiene

( ) ( ) ( )
( )( ) | |
( )
( ) | |
3
2
1 1
3
3 1 1
2 2
1
2
"
2
! 3 /
1 2 1
x
x
p P P
x P
x
P
x
P
x
P
i i i
i i i i
A +
A
+
= A +
A
+
A

A
=
+
+
u o o u

En la expresin anterior el error en trminos de ( )
3
x A no existira porque
1
=
3
y
por tanto habra que buscar el siguiente trmino del error, o sea que finalmente se
podra escribir

( ) ( ) ( )
( ) | |
( )
( ) | |
4
2
1 1
4
1
2 2
1
2
"
2 1 2 1
x
x
p P P
x P
x
P
x
P
x
P
i i i
i i i i
A +
A
+
= A +
A
+
A

A
=
+
+
u u

La expresin anterior es la ecuacin (2.79) obtenida antes.

116
Dos expansiones que sern muy tiles para expandir la ecuacin de difusividad
son

x
u u
u
i i
i
A

=
+ 2 / 1 2 / 1 '
................................................................................ (2.86)

x
u u
u
i i
i
A

=
+
+
1 '
2 / 1
.................................................................................. (2.87)

Con las expresiones anteriores se puede expandir la expresin siguiente de la
ecuacin de difusividad
|
|
.
|

\
|
x
P kA
x o
o
o
o
para flujo lineal en una dimensin.

Aplicando la expresin (2.86)

( ) ( ) ( ) ( )
x
x P kA x P kA
x
P
A k
x
i i i i
A

=
|
|
.
|

\
|
+ + 2 / 1 2 / 1 2 / 1 2 / 1
/ / / / o o o o
o
o
o
o
............... (2.88)

Aplicando ahora la expresin (2.87)

x
P P
x
P
i i
i
A

= |
.
|

\
|
+
+
1
2 / 1
o
o
............................................................................ (2.89)

x
P P
x
P
i i
i
A

= |
.
|

\
|

1
2 / 1
o
o
............................................................................ (2.90)


Llevando ahora las expresiones (2.89) y (2.90) a la ecuacin (2.88) se tiene


( ) ( ) ( ) ( )
2
1 2 / 1 2 / 1 2 / 1 1 2 / 1 1
) (
/ ) / / ( /
x
P kA P kA kA P kA
x
P
A k
x
i i i i i i
A
+ +
=
|
|
.
|

\
|
+ + +

o
o
o
o
............................................................................................................. (2.91)

La ecuacin (2.91) si k
x,i+1/2
= k
x,i-1/2
se convierte en

2
2
2
1 1
) (
2
x
P kA
x
P P P kA
x
P
A k
x
i i i
o
o
o
o
o
o
=
|
|
.
|

\
|
A
+
=
|
|
.
|

\
|
+


Las ecuaciones anteriores suponen que los pasos (x) son constantes

117

2.5-. ESQUEMAS DE APROXIMACIN EXPLICITO E IMPLICITO
Considrese que la Ecuacin 2.57 es vlida en el intervalo L x 0 < < y 0 > t . Para
dar solucin numrica a esta ecuacin, se requiere dividir el intervalo | | L , 0 en
diferencias finitas en el espacio, x A , y el intervalo | | t , 0 en diferencias finitas en el
tiempo, t A (Figura 2.10). As mismo, es necesario expandir las derivadas
2
2
x
P
c
c
y
t
P
c
c
.
La derivada
t
P
c
c
puede ser expandida mediante una aproximacin progresiva
(Ecuacin 2.65), regresiva (Ecuacin 2.71) central (Ecuacin 2.75). En las
Ecuaciones 2.65, 2.71 y 2.75 se aproxim la variacin de la presin con respecto a
la posicin a un tiempo fijo,
t
x
P
|
|
.
|
c
c
. En este caso se requiere evaluar la variacin

Figura 2.10-. Discretizacin de un Sistema Lineal en el Espacio y el Tiempo


de la presin con respecto al tiempo en una posicin fija,
x
t
P
|
|
.
|
c
c
. Por lo anterior,
en lugar de emplear la diferencia finita x A , tal como aparece en las Ecuaciones
2.65, 2.71 y 2.75 , se debe emplear la diferencia finita t A . As mismo, el subndice
At
Ax
i 1
L 0
t
i i + 1
118
i (el cual indica posicin) en este caso permanece fijo. El ndice que vara debe
ser aquel que indique el nivel de tiempo. Denotaremos por la letra n los ndices
que indiquen nivel de tiempo y sern superndices de la variable P . Tal como se
indic anteriormente, los ndices que indiquen posicin los denotaremos por la
letra i , y sern subndices de la variable P . Por ejemplo,
i
P indica presin en el
bloque i ,
n
P indica presin al tiempo n ,
n
i
P indica la presin en el bloque i al
tiempo n .

Teniendo en cuenta estos comentarios, las ecuaciones para las aproximaciones
progresiva, regresiva y central de
t
P
c
c
son, respectivamente,

( ) | | t O
t
P P
t
P
n
i
n
i
A +
A

=
c
c
+1
........................................................................... (2.92)

( ) | | t O
t
P P
t
P
n
i
n
i
A +
A

=
c
c
1
........................................................................... (2.93)

( ) | |
2
1 1
2
t O
t
P P
t
P
n
i
n
i
A +
A

=
c
c
+
...................................................................... (2.94)

En las Ecuaciones 2.92 a 2.94, los superndices 1 n , n y 1 + n representan los
tiempos
1 n
t ,
n
t y
1 + n
t , respectivamente.

Para la aproximacin de la primera derivada con respecto al tiempo, no
utilizaremos la aproximacin regresiva, Ecuacin 2.93, debido a que nuestro
inters es calcular la presin en cada bloque i del sistema a tiempos futuros (es
decir,
1 + n
i
P ), partiendo de valores actuales (es decir,
n
i
P ) y, contrario a nuestro
inters, la aproximacin regresiva relaciona la presin actual (
n
i
P ) con presiones a
tiempos pasados (
1 n
i
P ). La aproximacin central tampoco ser utilizada debido a
problemas de estabilidad, concepto que ser discutido ms adelante. Por estas
razones, se emplear la aproximacin progresiva, (Ecuacin 2.92), para aproximar
el trmino
t
P
c
c
.

Para aproximar la segunda derivada se aplicar la Ecuacin 2.79:

( )
2
1 1
2
2
2
x
P P P
x
P
i i i
i
A
+
~
|
|
.
|
c
c
+
.......................................................................... (2.79)
119

Esta ecuacin, tal como est escrita, no incluye el nivel de tiempo. La pregunta
que surge sera: A que tiempo hace referencia las presiones en la Ecuacin
2.79?. Dependiendo del nivel de tiempo asignado a cada trmino de presin en la
Ecuacin 2.79, se suele hablar de diferentes esquemas de aproximacin, de los
cuales los ms comunes son: la aproximacin explcita, la aproximacin implcita y
la aproximacin de Crank-Nicolson. Desde el punto de vista numrico, cada uno
de estos esquemas posee caractersticas diferentes.

Esquema de Aproximacin Explcita. Supngase que las presiones incluidas en
la Ecuacin 2.79 son evaluadas al tiempo
n
t . En este caso, se tiene:

( )
2
1 1
2
2
2
x
P P P
x
P
n
i
n
i
n
i
A
+
~
c
c
+
............................................................................ (2.95)

Si se lleva las Ecuaciones 2.92y 2,95 a la Ecuacin 2,57, se obtiene:

( )
t
P P
x
P P P
n
i
n
i
n
i
n
i
n
i
A

=
A
+
+
+
1
2
1 1
2
.................................................................. (2.96)

La Ecuacin 2.96 puede ser escrita como:

( ) ( )
2 2 1
1 1
2 x P x P t P t P t P
n
i
n
i
n
i
n
i
n
i
A A = A + A A
+
+


o bien:

( ) ( ) ( )
2
1
2 2
1
1
1
2
x
t
P
x
t
P
x
t
P P
n
i
n
i
n
i
n
i
A
A
+
|
|
.
|

\
|

A
A

A
A
=
+
+


Definiendo la variable como:

( )
2
x
t
A
A
= .................................................................................................... (2.97)

Se obtiene:

( )
n
i
n
i
n
i
n
i
P P P P
1 1
1
1 2
+
+
+ = ........................................................... (2.98)

A la Ecuacin 2.98 se le conoce como la aproximacin explcita a la Ecuacin
2.57. La aplicacin repetitiva de esta aproximacin numrica permite estimar la
presin en cada bloque i del sistema al tiempo
1 + n
t , con base en los valores de
120
presin en los bloques 1 i , i e 1 + i al tiempo
n
t . Para llevar a cabo esta serie
de clculos, se requiere de una condicin inicial y dos condiciones de frontera. A
continuacin se ilustra un ejemplo de aplicacin de la Ecuacin 2.98.

Ejemplo. Supngase un sistema lineal de siete bloques tal como el presentado en
la Figura 2.3. Si
0
t es tiempo inicial, obtener las expresiones de la aproximacin
explcita para el clculo de las presiones en cada bloque a un tiempo futuro
t t t A + =
0 1
y a un tiempo t t t A + =
1 2
, asumiendo que en cada caso las
condiciones de lmite son:

a. =
0
P constante
cero
P =
b. 0
8
= P

Solucin. Definiendo
n
t t =
0
, entonces
1 0 1 +
= A + =
n
t t t t . Debido a que
0
t es el
tiempo actual (condicin inicial), entonces la presiones a este tiempo,
( ) 0
1
P ,
( ) 0
2
P ,
( ) 0
7
P , son conocidas. De la Ecuacin 2.98 se pueden obtener las expresiones para
el clculo de las presiones en cada nodo al tiempo
1 1 +
=
n
t t :

Para el primer bloque, 1 = i :

( ) ( )
( )
( ) ( ) 0
2
0
1
0
0
1
1
1 2 P P P P + =

Asignando al bloque 0 = i la condicin de lmite izquierdo, se tiene =
0
P constante
cero
P = , por consiguiente:

( ) ( )
( )
( ) ( ) 0
2
0
1
0 1
1
1 2 P P P P
cero
+ =

Para el segundo bloque, 2 = i :

( ) ( )
( )
( ) ( ) 0
3
0
2
0
1
1
2
1 2 P P P P + =



Para el sptimo bloque, 7 = i :

( ) ( )
( )
( ) ( ) 0
8
0
7
0
6
1
7
1 2 P P P P + =

Asignando al bloque 8 = i la condicin de lmite derecho, se tiene 0
8
= P , por
consiguiente:
121

( )
( )
( ) ( ) 0
6
0
7
1
7
1 2 P P P + =

Al tiempo t t t A + =
1 2
, se pueden re-definir la variables
n
t y
1 + n
t de tal forma que
1
t t
n
= y
2 1
t t
n
=
+
. Teniendo en cuenta esta re-definicin de variables, de la
Ecuacin 2,98 se tiene:

Para el primer bloque, 1 = i :

( ) ( )
( )
( ) ( ) 1
2
1
1
1
0
2
1
1 2 P P P P + =

De nuevo, asignando al bloque 0 = i la condicin de lmite izquierdo, se tiene
cero
P P = = constante
0
, luego:

( )
( )
( ) ( ) 1
2
1
1
2
1
1 2 P P P P
cero
+ =

Para el segundo bloque, 2 = i :

( ) ( )
( )
( ) ( ) 1
3
1
2
1
1
2
2
1 2 P P P P + =



Para el sptimo bloque, 7 = i :

( ) ( )
( )
( ) ( ) 1
8
1
7
1
6
2
7
1 2 P P P P + =

Nuevamente, asignando al bloque 8 = i la condicin de lmite derecho, se tiene
0
8
= P , por consiguiente:

( )
( )
( ) ( ) 1
6
1
7
2
7
1 2 P P P + =

En la parte izquierda de la Figura 2.11 se indican los puntos en el espacio y los
niveles de tiempo tenidos en cuenta en la aproximacin explcita.

De la Ecuacin 2.98, y su aplicacin al ejemplo anterior, se infiere el procedimiento
general para el clculo de la distribucin de presiones mediante aplicacin de la
aproximacin explcita:

a. A cada uno de los bloques del sistema se le asigna un valor de presin al
tiempo inicial 0 = t teniendo en cuenta las condiciones inicial y de frontera. De
122
esta forma se obtienen valores para las variables
0
1 + i
P ,
0
i
P y
0
1 i
P ,
x
..., n i , 1 , 0 =
(la variable
x
n hace referencia al nmero variables discretas).

b. Se hace 0 = n y se calcula los valores de
1 + n
i
P para cada uno de los bloques
de la malla (
x
..., n i , 1 , 0 = ).

c. Una vez barrida toda la malla, se repite el procedimiento para el tiempo
2 + n
t ,
asignando a los valores de n los respectivos valores de 1 + n previamente
calculados y aplicando la Ecuacin 2,96 para 1 + n y 2 + n , en lugar de n y
1 + n , respectivamente.

Tal como se discutir ms adelante, para obtener estabilidad en este esquema de
solucin, se requieren incrementos de tiempo, t A , muy pequeos, lo que limita su
aplicacin en un gran nmero de problemas prcticos.



Figura 2.11-. Representacin Esquemtica de las Aproximaciones
Explcita e Implcita



Esquema de Aproximacin Implcita. Supngase que las presiones
1 i
P ,
i
P y
1 + i
P en la Ecuacin 2.79 hacen referencia al tiempo
1 + n
t . En este caso, la
Ecuacin 2.95 puede ser escrita como:

( )
2
1
1
1 1
1
2
2
2
x
P P P
x
P
n
i
n
i
n
i
A
+
~
c
c
+
+
+ +

....................................................................... (2.99)

t
n + 1
t
n
x
i
x
i 1
x
i + 1
Esquema Explcito
t
n + 1
t
n
x
i
x
i 1
x
i + 1
Esquema Implcito
123
Sustituyendo las Ecuaciones 2.92 y 2.99 en la Ecuacin 2.57, se obtiene:

( ) t
P P
x
P P P
n
i
n
i
n
i
n
i
n
i
A

=
A
+
+ +

+ +
+
1
2
1
1
1 1
1
2
.............................................................. (2.100)

La Ecuacin 2.100 puede ser escrita de la siguiente forma:

( )
n
i
n
i
n
i
n
i
P P P P = + +
+

+ +
+
1
1
1 1
1
1 2 ..................................................... (2.101)

A la Ecuacin 2.101 se le conoce como la aproximacin implcita a la Ecuacin
2.57. Si esta ecuacin se aplica a cada uno de los bloques de la malla, se genera
un sistema de ecuaciones cuya solucin simultnea permite obtener la distribucin
de presiones del sistema al tiempo
1 + n
t .

Ejemplo. Obtener el sistema de ecuaciones generado de la aproximacin
implcita, Ecuacin 2.101, para el ejemplo anterior.

Solucin. Definiendo
n
t t =
0
y
1 1 +
=
n
t t . La Ecuacin 2.101 genera una ecuacin
para cada bloque del sistema, tal como se indica a continuacin:

Para el primer bloque, 1 = i ,

( )
( )
( ) ( ) ( ) 0
1
1
2
1
1
1
0
1 2 P P P P = + +

Al igual que en el ejemplo anterior, asignando al bloque 0 i = la condicin de lmite
izquierdo, se tiene =
0
P constante
cero
P = , por consiguiente:

( )
( ) ( ) ( )
cero
P P P P = + +
0
1
1
2
1
1
1 2

Para el segundo bloque, 2 = i :

( )
( )
( ) ( ) ( ) 0
2
1
3
1
2
1
1
1 2 P P P P = + +



( )
( )
( ) ( ) ( ) 0
7
1
8
1
7
1
6
1 2 P P P P = + +

Al igual que en el ejemplo anterior, asignando al bloque 8 = i la condicin de lmite
derecho, se tiene 0
8
= P , por consiguiente:

( )
( )
( ) ( ) 0
7
1
7
1
6
1 2 P P P = +
124

Este sistema de ecuaciones puede ser escrito en forma matricial como
( )
F AP
1
= ,
donde A es la matriz de coeficientes definida como:

( )
( )
( )
( )
( )
( )
( )
(
(
(
(
(
(
(
(
(

+
+
+
+
+
+
+
=
1 2
1 2
1 2
1 2
1 2
1 2
1 2







A

P y F son los vectores de presin y trminos independientes, respectivamente,
definidos de la siguiente forma:

( )
( )
( )
( )
( )
( )
( )
( ) (
(
(
(
(
(
(
(
(

=
1
7
1
6
1
5
1
4
1
3
1
2
1
1
P
P
P
P
P
P
P
1
P ,
( )
( )
( )
( )
( )
( )
( ) (
(
(
(
(
(
(
(
(

=
0
7
0
6
0
5
0
4
0
3
0
2
0
1
P
P
P
P
P
P
P P
cero

F

La solucin del sistema de ecuaciones
( )
F AP
1
= permite obtener los valores de
presin en cada bloque al tiempo
1 1
t t
n
=
+
. Las tcnicas de solucin de este tipo
de sistemas, ms comnmente utilizadas en simulacin de yacimientos, sern
discutidas a lo largo de este texto, a partir del captulo cuatro.

La Figura 2.11 incluye los puntos en el espacio y los niveles de tiempo tenidos en
cuenta en la aproximacin implcita.

Aproximacin de Crank-Nicolson. Se obtiene de una combinacin de las
aproximaciones implcita y explcita; consiste en aproximar la ecuacin diferencial
en el punto |
.
|

\
|
+
2
1
, n i tal como se ilustra en la Figura 2.12.

La aproximacin de Crank-Nicolson aplicada a la Ecuacin 2.57 est dada por la
siguiente ecuacin:
125
Figura 2.12-. Representacin Esquemtica del Mtodo de Crank-Nicolson


2 1
2
2
1
2
2
2
1
2
1
+ +
|
|
.
|
c
c
=
|
|
.
|
c
c
+
|
|
.
|
c
c

n en Central
i
n
i
n
i
t
P
x
P
x
P
..................................................... (2.102)
Reemplazando los trminos diferenciales de la Ecuacin 2.102 por sus
correspondientes aproximaciones numricas, se obtiene:

( ) ( )
|
.
|

\
| A

=
(
(

A
+
+
(
(

A
+

+
+
+

+ +
+
2
2
2
2
1
2
2
1
1
2
1 1
2
1
1
1 1
1
t
P P
x
P P P
x
P P P
n
i
n
i
n
i
n
i
n
i
n
i
n
i
n
i


( ) ( ) t
P P
x
P P P
x
P P P
n
i
n
i
n
i
n
i
n
i
n
i
n
i
n
i
A

=
(
(

A
+
+
(
(

A
+

+
+
+

+ +
+
1
2
1 1
2
1
1
1 1
1
2
2
1
2
2
1
................ (2.103)

Obsrvese que la expansin numrica de la derivada de la presin con respecto al
tiempo, proviene de una aproximacin central en
2
1
+
=
n
t t y, en consecuencia,
puede ser escrita como:

...
6
2
2
2
2
1
2
2
2
1
+
+
|
|
.
|
c
c

|
.
|

\
| A

|
.
|

\
| A

=
c
c
n
i
n
i
n
i
t
P
t
t
P P
t
P
......................................................... (2.104)

Otras aproximaciones de expansin para la ecuacin (2.57) son:

La de Richardson dada por

t
n + 1
t
n
x
i
x
i 1
x
i + 1
( x
i
,

t
n +
)

126
( )
t
P P
x
P P P
n i n i n i n i n i
A

=
A
+
+ +
2
2
1 , 1 ,
2
, 1 , , 1
........................................................ (2.105)

La cual es una expansin implcita pero la expansin de la derivada del tiempo se
ha tomado como centrada; esta expansin implicara manejar informacin a tres
niveles de tiempo lo cual muchas veces puede no ser recomendable por razones
de requerimiento de espacio del computador.

La de Dufort Frankel dada por

( )
( )
t
P P
x
P P P P
n i n i n i n i n i n i
A

=
A
+ +
+ + , 1 ,
2
, 1 , , , 1
.................................................. (2.106)

la cual es tambin una expansin explcita.

Finalmente, tambin podra expandirse la ecuacin (2.57) de la siguiente forma


( )
( )
t
P P
x
P P P P
n i n i n i n i n i n i
A

=
A
+ +
+ + + + , 1 ,
2
1 , 1 1 , , , 1
.............................................. (2.107)

la cual es una forma implcita.


2.6 -. CRITERIOS TENIDOS EN CUENTA PARA LA SELECCION DE UN
ESQUEMA DE APROXIMACIN NUMERICA.
Existen diferentes esquemas de aproximacin numrica a las ecuaciones
diferenciales de flujo, algunos de los cuales han sido discutidos en el numeral 2.5.
La seleccin de uno u otro de estos esquemas se debe realizar teniendo en
cuenta, fundamentalmente, los siguientes criterios: error de truncamiento,
estabilidad, convergencia y consistencia. Esta parte del captulo se dedica a la
discusin del significado de estos conceptos y algunos mtodos para su anlisis.

2.6.1-. Error de Truncamiento. Tal como se menciona en la seccin 2.3, el error
de truncamiento hace referencia al error en que se incurre al utilizar la
aproximacin en diferencias finitas, luego de ser truncada, en lugar de la ecuacin
diferencial misma. En trminos generales el error de truncamiento se define
como:

(

=
Finitas s Diferencia
en Ecuacin
l Diferencia
Ecuacin
T
c .................................................. (2.108)

127
Es importante anotar que este error de truncamiento es diferente al error de
truncamiento en el que se incurre debido a las operaciones repetitivas del
computador. Los computadores tienen una memoria de capacidad finita y, en
consecuencia, pueden retener nicamente nmeros de "tamao finito" (es decir,
hasta cierto nmero de cifras significativas). El error cometido por esta clase de
truncamiento, el cual se propaga a medida que se incrementa el nmero de
operaciones, es diferente al error de truncamiento numrico a que hace referencia
la Ecuacin 2.108.

Ejemplo. Considrese la Ecuacin 2.57, cuya aproximacin explcita est dada
por la Ecuacin (2.96):

( )
t
P P
x
P P P
n
i
n
i
n
i
n
i
n
i
A

=
A
+
+
+
1
2
1 1
2
.................................................................. (2.96)

Aplicando la Ecuacin 2.108, se tiene:

( )
(
(

A
+

c
c

c
c
=
+
+
t
P P
x
P P P
t
p
x
p
n
i
n
i
n
i
n
i
n
i
T
1
2
1 1
2
2
2
c

Llevando las Ecuaciones 2.79 y 2.92, en reemplazo de
2
2
x
p
c
c
y
t
p
c
c
, a la anterior
ecuacin se tiene:

( )
( ) ( )

+
A

|
|
.
|
c
c

A

|
|
.
|
c
c

A
+
=
+
+
t
P P
x
x
P x
x
P
x
P P P
n
i
n
i
i i
n
i
n
i
n
i
T
1 4
6
6 2
4
4
2
1 1
! 6
2
! 4
2
2
c
( )
(
(

A
+

(
(

|
|
.
|
c
c
+
+
+
t
P P
x
P P P
t
t
P
n
i
n
i
n
i
n
i
n
i
i
1
2
1 1
2
2
2
2


Simplificando,

( )
( ) ( ) | | t x O
t
P t
x
P x
i i
T
A + A = +
|
|
.
|
c
c

|
|
.
|
c
c

A
=
2
2
2
4
4 2
2 12
c ............................ (2.109)

La ecuacin (2.109) indica que el error de truncamiento al hacer la expansin
explcita depende de (x)
2
y (t), pero es ms afectado por el tamao del paso en
t que por el tamao del paso en x; adems el truncamiento tiende a cero mientras
ms pequeos sean los pasos, o sea es vlida la expansin explcita.

Encontremos el error de truncamiento en la expansin de Crank Nicolson:
128

Llevando las Ecuaciones 2.57 y 2.103 a la definicin de error de truncamiento,
Ecuacin 2.108, se tiene:

( ) ( )
(
(

A
+
+
A
+

(

c
c

c
c
=
+
+
+

+ +
+
t
P P
x
P P P
x
P P P
t
p
x
p
n
i
n
i
n
i
n
i
n
i
n
i
n
i
n
i
T
1
2
1 1
2
1
1
1 1
1
2
2
2
2
1
2
2
1
c
..................................................................................................................... (2.110)

Combinando las Ecuaciones (2.76), (2.92), (2.95), y (2.99), se obtiene:

( )
( )
( )
2
1 1
4
4 2
2
1
1
1 1
1
2
2
1
12 2
1
2
2
1
x
P P P
x
P x
x
P P P
n
i
n
i
n
i
i
n
i
n
i
n
i
T
A
+
+

|
|
.
|
c
c

A

A
+
=
+
+

+ +
+
c

( )
(
(
(
(
(

|
|
.
|
c
c

|
.
|

\
| A
+
|
.
|

\
| A

|
|
.
|
c
c

A

+
+

2
1
3
3
2
1
4
4 2
6
2
2
2
12 2
1
n
i
n
i
n
i
i
t
P
t
t
P P
x
P x

( ) ( )
(
(

(
A

A
+
+
A
+

+
+
+

+ +
+
t
P P
x
P P P
x
P P P
n
i
n
i
n
i
n
i
n
i
n
i
n
i
n
i
1
2
1 1
2
1
1
1 1
1
2
2
1
2
2
1


Simplificando,

( ) ( )
+
|
|
.
|
c
c

|
.
|

\
| A
+
|
|
.
|
c
c

|
|
.
|
c
c

A
=
+ +
2
1
3
3
2
4
4 2
1
4
4 2
6
2
24 24
n
i
n
i
n
i
T
t
P
t
x
P x
x
P x
c

Es decir, el orden del error de truncamiento de la aproximacin de Crack-Nicolson
est dado por:

( ) ( ) | |
2 2
t x O
T
A + A = c ................................................................................ (2.111)

Al comparar las Ecuaciones (2.109) y (2.111), se infiere que la aproximacin de
Crank-Nicolson presenta un menor error de truncamiento que las aproximaciones
implcita y explcita.

2.6.2-. Estabilidad. Anlisis de estabilidad es un procedimiento a travs del cual
es posible determinar si el error obtenido en un punto (bloque) de la malla, en un
momento determinado, disminuye o aumenta al incrementar el tiempo. Si el error
disminuye, se dice que la aproximacin es estable. Si aumenta, se dice que es
inestable. Si el error disminuye nicamente bajo ciertas condiciones, se dice que
es condicionalmente estable. Para que una aproximacin numrica conlleve a
129
resultados tiles se requiere que dicha aproximacin sea estable. La Figura 2.13
ilustra la variacin de presin en un punto determinado del sistema en funcin del
tiempo para un esquema estable y un esquema inestable. A continuacin se
presentan dos criterios muy tiles para determinar estabilidad: El criterio de
Karplus y el Mtodo de Von Neumann.

El Criterio de Karplus. El mtodo de Karplus no considera el efecto de
condiciones de lmite en el anlisis de estabilidad. El procedimiento a seguir para
la aplicacin del mtodo puede ser resumido de la siguiente forma:

a. La ecuacin en diferencias finitas es reordenada de manera tal que tome la
siguiente forma:

( ) ( ) ( ) ( ) 0
* 1 * 1 *
1
*
1
= + + + + +
+
+

i
n
i i
n
i i
n
i i
n
i
P P d P P c P P b P P a ... (2.112)
En la Ecuacin 2.112,
*
i
P indica la presin en el punto i y al tiempo en torno al
cual se efecta la aproximacin numrica.


Figura 2.14-. Comportamiento en el Cambio de Presin en un Punto
Determinado de un Esquema Estable y un Esquema Inestable


b. La aproximacin ser estable cuando todos los coeficientes a, b, c, d, sean
positivos o si alguno de ellos es negativo, cuando la suma de los coeficientes
a , b , c , ..., d , ... sea menor o igual a cero.

Nivel de tiempo, t
n
n 1 n
P P P =
+
A
n
P A
Esquema Estable
Esquema Inestable
130
Debido a que los coeficientes a , b , c , ..., d pueden cambiar de un nodo a otro, el
criterio debe ser aplicado a cada nodo. Un esquema determinado puede ser
estable en un nodo pero inestable en otros.

Ejemplo. Con la finalidad de aplicar el criterio de Karplus a la aproximacin
explcita, la Ecuacin 2.96 puede ser escrita de la siguiente forma:

( ) ( ) ( ) 0
1
1 1
= +
+
+
n
i
n
i
n
i
n
i
n
i
n
i
P P P P P P ............................................. (2.113)

Esta aproximacin ser estable siempre y cuando se cumpla que:

( )
( )
0 1 2 1 2 1
2
s
A
A
= = +
x
t


Despejando para ( )
2
x t A A , se obtiene:
( )
( )
2
1
2
s
A
A
x
t
................................................................................................... (2.114)

Por consiguiente, la aproximacin explcita es condicionalmente estable y la
condicin de estabilidad est dada por la Ecuacin 2.114 .

Siguiendo un procedimiento completamente anlogo es posible determinar que la
aproximacin implcita es estable. La expresin para la expansin implcita es la
ecuacin (2.100)

( )
n
i
n
i
n
i
n
i
P P P P = + +
+

+ +
+
1
1
1 1
1
1 2 ........................................... (2.101)

La cual se puede llevar a

( ) ( )( ) ( ) 0 1 2
1
1
1 1
1
= + +
+

+ +
+
n
i
n
i
n
i
n
i
n
i
n
i
P P P P P P

En la expresin anterior el coeficiente ( ) 1 2 + es menor que cero y por tanto
para que se tenga estabilidad segn el criterio de Karplus se requiere que

( ) 0 1 0 1 2 s s + + ................................................................ (2.115)

La expresin (2.115) siempre se va a cumplir y por tanto la expansin implicita es
incondicionalmente estable.

Por su parte la estabilidad de la expansin de Crank Nicolson se puede probar
por el criterio de Karplus de la siguiente manera:

La expresin para esta expansin es la ecuacin (2.103)
131

( ) ( ) t
P P
x
P P P
x
P P P
n
i
n
i
n
i
n
i
n
i
n
i
n
i
n
i
A

=
(
(

A
+
+
(
(

A
+

+
+
+

+ +
+
1
2
1 1
2
1
1
1 1
1
2
2
1
2
2
1
................. (2.103)

La cual tambin se puede escribir como

n
i
n
i
n
i
n
i
n
i
n
i
n
i
n
i
P P P P P P P P
=
(
(

+
+
(
(

+
+
+
+

+ +
+
1
1 1
1
1
1 1
1
2
2 2


Y llevndola a una forma similar a la ecuacin (2.103) se tiene

( ) ( )( ) ( ) ( ) ( ) 0 1 2
1 1
1
1
1 1
1
= + + + +
+
+
+
+ +
+
n
i
n
i
n
i
n
i
n
i
n
i
n
i
n
i
n
i
n
i
P P P P P P P P P P

Y como en la expresin anterior hay un coeficiente menor que cero, segn Karplus
para que la expansin dada por la ecuacin (2.103) se requiere que

( )
2
1
0 3 1 2 s s + + ............................................................... (2.116)
La condicin dada por la ecuacin (2.116) es la misma dada para la expansin
explcita, ecuacin (2.115).

El Mtodo de Von Neumann o Anlisis Armnico. Una de las suposiciones
bsicas de este mtodo se fundamenta en el hecho de que el error de una
aproximacin numrica satisface la aproximacin numrica. El valor exacto de
presin en un punto
i
x a un tiempo
n
t , ( )
n i
t , x p , y el valor calculado al aplicar una
aproximacin numrica,
n
i
P ,estn relacionados mediante la siguiente expresin:

( )
n
i
n
i n i
P t , x p c + = ....................................................................................... (2.117)

En la Ecuacin 3.58
n
i
c es el error de truncamiento en el cual se incurre al
aproximar el valor exacto ( )
n i
t , x p al valor obtenido de la aproximacin numrica,
n
i
P .

Los valores de presin obtenidos de la solucin exacta satisfacen la aproximacin
numrica, pues la solucin exacta tiene el mismo valor de la aproximacin
numrica cuando el error tiende a cero. A manera de ejemplo, considrese la
aproximacin explcita a la Ecuacin 2.57. Puesto que la solucin exacta satisface
la aproximacin numrica, de la Ecuacin 2.96 se obtiene:

132
( ) ( ) ( )
( )
( ) ( )
t
x p x p
x
x p x p x p
i i i i i
A