MATEMATICAS AVANZADAS

por

Antonio Flores Tlacuahuac

Departmento de Ingenierı́a Quı́mica

Universidad Iberoamericana
México, DF

October 22, 2008

Chapter 1

Ecuaciones Diferenciales
Ordinarias: Problemas de Valores
en la Frontera

1.1 Introducción
En este capitulo abordaremos el problema de la solución numerica de ODE’s donde
el valor de la variable dependiente se especifica en algun(os) lugar(es) del espacio de
solución. Ejemplo de este tipo de problemas es la solución de la ecuación de difusión
d2 c
, 06x61 (1.1.1)
dx2
sobre la cual se imponen las siguientes restricciones
c(0) = 0 y c(1) = 1
(1.1.2)
las restricciones se conocen como las condiciones frontera del problema (dado que
estan definidas sobre los lı́mites del dominio de solución), y significa que la solución
determinada c(x) deberá tomar los valores c(0) cuando x = 0, y c(1) = 1 cuando x = 1.
Note que el problema planteado por las ecuaciones 1.1.1 y 1.1.2 no puede ser resuelto
de manera directa, por tecnicas de valores iniciales; por lo cual se deben usar un nuevo
tipo de tecnicas de solución.
Las condiciones frontera pueden clasificarse de la siguiente forma
Dirichlet : y(x) = f1 (x)
dy(x)
N eumann : = f2 (x)
dx
dy(x)
Robin : ay(x) + b = f3 (x)
dx

1
2CHAPTER 1. ECUACIONES DIFERENCIALES ORDINARIAS: PROBLEMAS DE VALORES EN L

como puede notarse de la definición anterior, las condiciones frontera de Dirichlet

simplemente definen o especifican el valor de la solución en alguna frontera del dominio
de solución. Las condiciones frontera de Neumann involucran la especificación de la
derivada de la solución sobre el mismo espacio de solución. Finalmente, las condiciones
frontera de Robin son una mezcla del tipo de condiciones frontera de Dirichlet y de
Neumann. Las condiciones frontera también se denominan de primera clase si son
condiciones de Dirichlet, de segunda clase si son condiciones de Neumann, y de tercera
clase si son condiciones de Robin.

1.2 Colocación Ortogonal.

El método de colocación ortogonal presenta algunas ventajas sobre el método de colo-
cación convencional: a) la función de aproximación a la solución involucra el empleo de
polinomios ortogonales, b) los puntos de colocación son simplemente las raices de uno
de estos polinomios ortoganales (evitando asi una posible selección pobre de los puntos
de colocación), y c) el método se plantea de tal modo que la variable dependiente son
los valores de la solución, en los puntos de colocación, en vez de los coeficientes de la
solución aproximada.
Procediendo de manera semejante a la formulación del MRP, la solución y(x) se
aproxima mediante la siguiente expansión finita
N
X
y(x) ≈ ŷ(x) = ai yo (x) (1.2.3)
i=1

donde N representa los puntos de colocación, yi (x) denota el conjunto de funciones

base, ai son los coeficientes de la expansión que deseamos determinar. Es claro que,
una vez que los coeficientes ai sean determinados, podemos evaluar la solución, en
cualquier punto xj del espacio de solución, usando la ecuación de interpolación
N
X
y(xj ) ≈ ŷ(xj ) = ai yi (xj ) (1.2.4)
i=1

ahora la ecuación 1.2.4 se evalua en N puntos xj , podemos usar el N xN sistema de

ecuaciones resultante, para determinar el conjunto de coeficientes ai
N
X
ai = [yi (xj )]−1 [y(xj )] (1.2.5)
i=1

esto, por supuesto, supone que yi (xj ) forma un conjunto de ecuaciones linealmente
independientes, tal que la inversa de yi (xj ) existe. La ecuación anterior podrı́a per-
mitirnos usar el conjunto de soluciones en los puntos de colocación y(xj ) como las
variables desconocidas, en vez del conjunto de coeficientes ai .
Dado que en los problemas acerca de la solución numerica de ODE’s aparecen
derivadas de primer orden (al menos) y de segundo y ordenes superiores, podemos
Colocación Ortogonal. 3

expresar dichas derivadas también en terminos de la expansión dada por la ecuación

1.2.4
N
dy(xj ) dŷ(xj ) X dyi (xj )
≈ = ai (1.2.6)
dx dx i=1
dx
de manera semejante
N
d2 y(xj ) d2 ŷ(xj ) X d2 yi (xj )
≈ = ai (1.2.7)
dx2 dx2 i=1
dx2

las ecuaciones 1.2.6 y 1.2.7 tambiı́n pueden plantearse en terminos de la solución en

los puntos de colocación y(xj ) en vez del conjunto de coeficientes ai . Substituyendo la
ecuación 1.2.5 en las ecuaciones 1.2.6 y 1.2.7
N N
dy(xj ) dŷ(xj ) X X dyi (xj )
≈ = [yi (xk )]−1 [y(xk )] (1.2.8)
dx dx i=1 k=1
dx

y
N N
d2 y(xj ) d2 ŷ(xj ) X X −1 d2 yi (xj )
≈ = [y i (x k )] [y(x k )] (1.2.9)
dx2 dx2 i=1 k=1
dx2
las ecuaciones 1.2.8 y 1.2.9 pueden reescribirse como
N
dy(xj ) dŷ(xj ) X
≈ = Ajk y(xk ) (1.2.10)
dx dx i=1

y
N
d2 y(xj ) d2 ŷ(xj ) X
≈ = Bjk y(xk ) (1.2.11)
dx2 dx2 i=1

donde
XN
dyi (xj )
Ajk = [yi (xk )]−1 (1.2.12)
k=1
dx
N
X d2 yi (xj )
Bjk = [yi (xk )]−1 (1.2.13)
k=1
dx2
A continuación las funciones base se seleccionan como polinomios ortogonales, donde
el polinomio Pm se define como una combinación lineal de potencias de x siendo m el
exponente de x de mayor grado
m
X
Pm (x) = c j xj (1.2.14)
j=0

por ejemplo si m = k, entonces el conjunto de coeficientes cj se selecccionan de manera

tal que el polinomio P1 sea ortogonal al polinomio Po ; del mismo modo P2 debe ser
ortogonal a Po y P1 . Continuando de esta manera finalmente tenemos que el polinomio
4CHAPTER 1. ECUACIONES DIFERENCIALES ORDINARIAS: PROBLEMAS DE VALORES EN L

Pk deberá ser ortogonal a todos los polinomios de grado menor que k. En algunos casos
la condición de ortogonalidad puede incluir alguna función peso W (x)
Z b
Pm (x)Pk (x)W (x)dx = 0, k = 0, 1, ..., m − 1 (1.2.15)
a

1.2.1 Determinación de raices de un polinomio ortogonal

A manera de ejemplificar la manera como se determinan los puntos de colocación (los
cuales corresponden a las raices de un polinomio ortogonal) consideremos el siguiente
ejemplo.
Ejemplo 1 Determine los coeficientes b,c y d del conjunto de polinomios dados por Po = 1, P1 =
1 + bx, P2 = 1 + cx + dx2 requiriendo que P1 sea ortogonal a Po y P2 sea ortogonal P1 y Po . Suponga
W (x) = 1 y a = 0, b = 1.
Solución.
Dado que P1 debe ser ortogonal a Po
Z 1 Z 1
P1 Po dx = (1 + bx)dx = 0
0 0

cuya solución es b = −2. De manera semejante, ya que P2 debe ser ortogonal a Po y P1

Z 1 Z 1
P2 Po dx = (1 + cx + dx2 )dx = 0
0 0

Z 1 Z 1
P2 P1 dx = (1 + cx + dx2 )(1 − 2x)dx = 0
0 0

resolviendo simultaneamente este conjunto de ecuaciones se tiene c = −6, d = 6. Por lo tanto los
polinomios ortogonales Po , P1 y P2 estan dados por
Po = 1
P1 = 1 − 2x
P2 = 1 − 6x − 6x2

las raices de P2 están ubicadas en: 0.2113 y 0.7887.

Debe notarse que la elección del intérvalo de integración (0,1) es completamente arbi-
traria. Uno podrı́a seleccionar cualquier otro intérvalo de interes (0,tf ) y la ubicación
de las raices serı́a la misma empleando el correspondiente cambio de coordenadas tal
como se muestra en el siguiente ejemplo.
Ejemplo 2 Repetir el ejemplo 1 donde supongamos, por ejemplo, que seleccionamos arbitraria-
mente tf = 12.5.
•
Z tf Z tf
P1 Po dx = (1 + bx)dx = 0
0 0

integrando,
µ ¶¯12.5
bx2 ¯¯
x+ = 12.5 + 78.125b = 0 (1.2.16)
2 ¯0
Colocación Ortogonal. 5

• similarmente,
Z tf Z tf
P2 Po dx = (1 + cx + dx2 )dx = 0
0 0

integrando,
µ ¶¯12.5
cx2 dx3 ¯¯
x+ + = 12.5 + 78.125c + 651.0417d = 0 (1.2.17)
2 3 ¯0

si formamos el siguiente polinomio ortogonal,

Z tf Z tf
P1 P2 dx = (1 + bx)(1 + cx + dx2 )dx = 0
0 0

integrando,
µ ¶¯12.5
bdx4 (d + bc)x3 (c + b)x2 ¯
+ + + x ¯¯ =
4 3 2 0
6103.5bd + 651.0417(d + bc) + 78.125(c + b) + 12.5 = 0 (1.2.18)

resolviendo las ecuaciones 1.2.16, 1.2.17 y 1.2.18:

b = −0.16, c = −0.48, d = 0.0348

por lo tanto el polinomio ortogonal de segundo orden estará dado por:

P2 (x) = 1 − 0.48x + 0.0384x2

cuyas raices son:

x1 = 9.8584, x2 = 2.6416

si dividimos x1 y x2 entre tf obtendremos las mismas raices que el correspondiente polinomio ortogonal
del ejemplo 1:
9.8554
= 0.7887
12.5
2.6416
= 0.2113
12.5

En el siguiente ejemplo se muestra la aplicación del método de colocación ortogonal

para resolver un problema de valores iniciales.

Ejemplo 3 Simular la operación dinámica de un reactor batch isotérmico donde ocurre la reacción:
2A → B

el model matemático está dado por:

dCA 2
= −kCA , CA (0) = 3 mol/lt
dt
donde k = 0.1 lt/(mol-min); simular la operación del reactor durante los primeros 60 minutos de
operación.
El problema se resolverá usando dos puntos internos de colocación tal como se mustra en la Figura
1.1.
6CHAPTER 1. ECUACIONES DIFERENCIALES ORDINARIAS: PROBLEMAS DE VALORES EN L

u u

x0 x1 x2 x3
0 1

Figure 1.1: El espacio de solución se divide en dos puntos de colocación internos.

el punto x0 representa la solución dada por las condiciones iniciales del problema; o sea en este punto
la concentración del compuesto A es igual a 3 mol/lt. Los puntos x1 y x2 representan la concentración
del mismo compuesto A en los 2 puntos internos de colocación, mientras que el punto x3 representa
la concentración cuando se ha alcanzado el lı́mite superior de integración (1 o 60 dependiendo de si
se usa el enfoque adimensional o no adimensional, respectivamente).
Para verificar que la solución que se obtiene es la misma, independientemente del intervalo de
integración seleccionado, el problema se resuelve usando el intérvalo adimensional (0,1) de operación.
Las soluciones numéricas se comparan contra la solución analı́tica del problema dada por:

CAo
CA =
1 + ktCAo

Usando el intérvalo adimensional (0,1) las raices estarán úbicadas en:

xo = 0
x1 = 0.2113
x2 = 0.7887
x3 = 1

la matriz de primeras derivadas parciales de polinomios de Lagrange lineales está dada por:
 
−2.7321 1.7321 1.7321 −0.7321
A =  0.7321 −1.7321 −1.7321 2.7321 
−1.0 2.1962 −8.1962 7.0

el correspondiente modelo discretizado está dado por:

2
A11 CAo + A12 CA1 + A13 CA2 + A14 CA3 = −kCA1 tf
2
A21 CAo + A22 CA1 + A23 CA2 + A24 CA3 = −kCA2 tf
2
A31 CAo + A32 CA1 + A33 CA2 + A34 CA3 = −kCA3 tf

donde tf = 60 mins representa el factor de escala de tiempo para que el intérvalo de integración esté
comprendido en el rango adimensional (0,1). Reemplazado valores:
2
−8.1963 + 1.7321CA1 + 1.7321CA2 − 0.7321CA3 = −6CA1
2
2.1963 − 1.7321CA1 + 1.7321CA2 + 2.7321CA3 = −6CA2
2
−3 + 2.1962CA1 − 8.1962CA2 + 7CA3 = −6CA3
Colocación Ortogonal. 7

N xj Wj
1 .5000000000 .6666666667
2 .2113248654 .5000000000
.7886751346 .5000000000
3 .1127016654 .2777777778
.5000000000 .4444444444
.8872983346 .2777777778
4 .0694318442 .1739274226
.3300094783 .3260725774
.6699905218 .3260725774
.9305681558 .1739274226
5 .0469100771 .1184634425
.2307653450 .2393143353
.5000000000 .2844444444
.7692346551 .2393143353
.9530899230 .1184634425
6 .0337652429 .0856622462
.1693953068 .1803807865
.3806904070 .2339569678
.6193095931 .2338569678
.8306046933 .1803807865
.9662347571 .0856622462

Table 1.1: Raices de polinomios ortogonales y funciones de peso definidos por la

ecuación ...W=1,a=0,b=1

En el método de colocación ortogonal, los m puntos de colocación corresponden

a las m raices del polinomio de mayor grado Pm (x). Esto define completamente la
localización de los puntos de colocación una vez que se han seleccionado la función
peso W (x), y el intervalo de integración a, b. La tabla 1.1 muestra la localización de
los puntos de colocación, junto con las funciones peso, definidos por la ecuación 1.2.15.
Notese que los puntos de colocación depende de la elección especifica de la función
peso, y del intervalo de colocación.

1.2.2 Sistemas sin simetria

El siguiente paso en la aplicación del método de colocación ortogonal consiste en pro-

poner una función de aproximación que tenga una relación estrecha con el problema
partı́cular que nos interesa resolver. Por ejemplo, podriamos proponer una función de
aproximación donde el primer termino satisface las condiciones frontera, y los siguientes
terminos (que son parte de una serie) satisfacen las condiciones frontera homogeneas.
8CHAPTER 1. ECUACIONES DIFERENCIALES ORDINARIAS: PROBLEMAS DE VALORES EN L

Entonces, una posible elección serı́a

N
X
y ≈ ŷ = x + x(1 − x) ai Pi−1 (x) (1.2.19)
i=1

que puede reescribirse como

N
X +2
y ≈ ŷ = bi Pi−1 (x) (1.2.20)
i=1

suponiendo que Pi−1 (x) representa un polinomio de orden i − 1, la ecuación anterior

estará dada por
N
X +2
y ≈ ŷ = di xi−1 (1.2.21)
i=1
la ecuación 1.2.21 define la forma básica de la función de aproximación en terminos
de una familia de polinomios ortogonales. Esta ecuación debe entonces escribirse para
cada uno de los puntos de colocación. Como se dijo antes, los N puntos de colocación
xj estarán dados por las N raices del polinomio ortogonal de mayor grado PN (x).
Normalmente x1 y xN +2 corresponden a los puntos, del espacio de solución, donde se
especifican las condiciones frontera (por ejemplo, si el problema está adimensionalizado
entre 0 y 1, x1 = 0, xN +2 = 1); los restantes puntos de colocación xj (2 6 j 6 N + 1) se
denominan puntos internos de colocación. Entonces, en cada punto de colocación xj ,
la función de aproximación y sus primera y segunda derivadas están dadas por
N
X +2
y(xj ) ≈ ŷ(xj ) = di xi−1
j (1.2.22)
i=1
N
X +2
dy(xj ) dŷ
≈ = di (i − 1)xi−2
j (1.2.23)
dx dxj i=1
N
X +2
d2 y(xj ) d2 ŷ
≈ 2 = di (i − 1)(i − 2)xji−3 (1.2.24)
dx2 dxj i=1
o bien en notación matricial
y = Qd (1.2.25)
dy
= Cd (1.2.26)
dx
d2 y
= Dd (1.2.27)
dx2
(1.2.28)
donde
Qji = xi−1
j (1.2.29)
Cji = (i − 1)xi−1
j (1.2.30)
Dji = (i − 1)(i − 2)xji−3 (1.2.31)
(1.2.32)
Colocación Ortogonal. 9

N W   A   B 
1
6
-3 4 -1 4 -8 4
1  2   -1 0 1   4 -8 4 
3
1
 6   1 -4 3   4 -8 4 
0 -7 8.196 -2.196 1 24 -37.18 25.18 -12
 1   -2.732 1.732 1.732 -.7321   16.39 -24 12 -4.392 
2  2     
 1   .7321 -1.7321 -1.7321 2.732   -4.392 12 -24 16.39 
2
1
2
-1 2.196 -8.196 7 -12 25.18 -37.18 24

Table 1.2:

por facilidad computacional se acostumbra escribir la primera y segunda derivada, de

la función de aproximación, en terminos de la solución en los puntos de colocación
y(xj ) en vez del conjunto de coeficientes di . Entonces

dy
= CQ−1 y ≡ Ay (1.2.33)
dx
d2 y
= DQ−1 y ≡ By (1.2.34)
dx2
(1.2.35)

la tabla 1.2 muestra las matrices A y B para diversos puntos de colocación. Las
suposiciones hechas para obtener estas matrices son las mismas que aquellas hechas
para determinar la localización de los puntos de colocación mostrados en la tabla 1.1.
Finalmente las funciones peso W (x), necesarias para la evaluación de la localización
de los puntos de colocación (ver ecuación 1.2.15) pueden evaluarse de la formula de
cuadratura
Z 1 N
X +2
f (x)dx = Wj f (xj ) (1.2.36)
0 j=1

por ejemplo si f (x) = xi−1 entonces

Z 1 N
X +2
1
x i−1
dx = Wj xji−1 = (1.2.37)
0 j=1
i

o en notación matricial
WQ = f (1.2.38)
de donde
W = f Q−1 (1.2.39)
La aplicación del método de colocación ortogonal se ilustra en el siguiente ejemplo.

Ejemplo 4 Reactor catalı́tico tubular con dispersión axial.

Considere el problema de determinar los perfiles de concentración y temperatura en
un reactor tubular empacado con material catalı́tico. Los reactivos fluyen a través del
10CHAPTER 1. ECUACIONES DIFERENCIALES ORDINARIAS: PROBLEMAS DE VALORES EN

catalizador, reaccionan sobre la superficie de este, y los productos fluyen sobre el resto
del lecho catalı́tico. Normalmente, debido al movimiento de los materiales alrededor del
catalizador, ocurre dispersión del material a lo largo de las direcciones radial y axial.
En este ejemplo suponemos que la dispersión radial es muy pequen̄a, en comparación
con la dispersión axial, de tal forma que podemos despreciar su efecto. El modelo
estacionario no-isotérmico del reactor catalitico tubular con dispersión axial está dado
por:
d2 c dc
2
− P eM − P eM R(c, T ) = 0 (1.2.40a)
dz dz
d2 T dT
2
− P eH − P eH βR(c, T ) = 0 (1.2.40b)
dz dz
sujeto a las condiciones frontera:
dc
= P eM (c − 1), en z = 0 (1.2.41a)
dz
dc
= 0, en z = 1 (1.2.41b)
dz
dT
= P eH (T − 1), en z = 0 (1.2.41c)
dz
dT
= 0, en z = 1 (1.2.41d)
dz
donde, γ
R = kc2 eγ− T (1.2.42)

Aplicando el método de colocación ortogonal al conjunto de ecuaciones 1.2.40a-1.3.55

tenemos:1
N
X +2 N
X +2
Bij cj − P eM Aij cj − P eM R(ci , Ti ) = 0, i = 2, ..., N + 1 (1.2.43a)
j=1 j=1

N
X +2 N
X +2
Bij Tj − P eH Aij Tj − P eH βR(ci , Ti ) = 0, i = 2, ..., N + 1 (1.2.43b)
j=1 j=1

N
X +2
A1j cj = P eM (c1 − 1), en z = 0 (1.2.43c)
j=1

N
X +2
AN +2,j cj = 0, en z = 1 (1.2.43d)
j=1

N
X +2
A1j Tj = P eH (T1 − 1), en z = 0 (1.2.43e)
j=1

1
en la formulación mostrada N representa el número de puntos internos de colocación.
Colocación Ortogonal. 11

u u u

x1 x2 x3 x4 x5
0 1

Figure 1.2: Aplicación del método de colocación ortogonal.

N
X +2
AN +2,j Tj = 0 en z = 1 (1.2.43f)
j=1

como puede notarse, el método de colocación ortogonal ha permitido transformar un

conjunto de ecuaciones diferenciales ordinarias no-lineales, ecuaciones 1.2.40a-1.3.55
y sus condiciones frontera, en un conjunto de ecuaciones algebraicas simultaneas no-
lineales dado por las ecuaciones 1.2.43a-1.2.43f. La naturaleza no-lineal de estas ecua-
ciones es debida a la presencia del termino velocidad de reacción R(c, T ) que acopla
a los balances de materia y de energı́a. El problema se resuelve ahora para 3 puntos
internos de colocación.

• Solución empleando 3 puntos de colocación. Expandiendo el grupo de

ecuaciones 1.2.43a-1.2.43f (ver figura 1.2),

A11 c1 + A12 c2 + A13 c3 + A14 c4 + A15 c5 − P eM (c1 − 1) = 0 (1.2.44a)

B21 c1 + B22 c2 + B23 c3 + B24 c4 + B25 c5 − P eM [A21 c1 + A22 c2

γ− Tγ
+A23 c3 + A24 c4 + A25 c5 ] − P eM kc22 e 2 =0 (1.2.44b)

B31 c1 + B32 c2 + B33 c3 + B34 c4 + B35 c5 − P eM [A31 c1 + A32 c2

γ− Tγ
+A33 c3 + A34 c4 + A35 c5 ] − P eM kc23 e 3 =0 (1.2.44c)

B41 c1 + B42 c2 + B43 c3 + B44 c4 + B45 c5 − P eM [A41 c1 + A42 c2

γ− Tγ
+A43 c3 + A44 c4 + A45 c5 ] − P eM kc24 e 4 =0 (1.2.44d)

A51 c1 + A52 c2 + A53 c3 + A54 c4 + A55 c5 = 0 (1.2.44e)

A11 T1 + A12 T2 + A13 T3 + A14 T4 + A15 T5 − P eH (T1 − 1) = 0 (1.2.44f)
12CHAPTER 1. ECUACIONES DIFERENCIALES ORDINARIAS: PROBLEMAS DE VALORES EN

B21 T1 + B22 T2 + B23 T3 + B24 T4 + B25 T5 − P eH [A21 T1 + A22 T2

γ− Tγ
+A23 T3 + A24 T4 + A25 T5 ] − P eH βkc22 e 2 = 0 (1.2.44g)

B31 T1 + B32 T2 + B33 T3 + B34 T4 + B35 T5 − P eH [A31 T1 + A32 T2

γ− Tγ
+A33 T3 + A34 T4 + A35 T5 ] − P eH βkc23 e 3 = 0 (1.2.44h)

B41 T1 + B42 T2 + B43 T3 + B44 T4 + B45 T5 − P eH [A41 T1 + A42 T2

γ− Tγ
+A43 T3 + A44 T4 + A45 T5 ] − P eH βkc24 e 4 =0 (1.2.44i)

A51 T1 + A52 T2 + A53 T3 + A54 T4 + A55 T5 = 0 (1.2.44j)

Las matrices de colocación estan dadas por

 
-13.00 14.79 -2.67 1.88 -1.00
 -5.32 3.87 2.07 -1.29 0.68 
 
A=
 1.50 -3.23 0.00 3.23 -1.50 
 (1.2.45)
 -0.68 1.29 -2.07 -3.87 5.32 
1.00 -1.88 2.67 -14.79 13.00

 
84.00 -122.06 58.67 -44.60 24.00
 53.24 -73.33 26.67 -13.33 6.76 
 
B=
 -6.00 16.67 -21.33 16.67 -6.00 
 (1.2.46)
 6.76 -13.33 26.67 -73.33 53.24 
24.00 -44.60 58.67 -122.06 84.00

dos casos, referentes al mismo problema, fueron resueltos.

γ
1. P eM = P eH = 2, R = 3.36c2 eγ− T , γ = 17.6, β = −.056

Nodo 1 2 3 4 5
Concentracción .58031396 .49549975 .31429452 .23989384 .23537842
Temperatura 1.0234997 1.0282494 1.0383976 1.0425647 1.0428180
γ
2. P eM = P eH = 96, R = 3.817037c2 eγ− T , γ = 17.6, β = −.056

Nodo 1 2 3 4 5
Concentracción .96294732 .65499445 .23656508 .13525980 .12590711
Temperatura 1.0020749 1.0193202 1.0427526 1.0484249 1.0489490
Colocación Ortogonal. 13

El sistema de ecuaciones algebraicas no-lineales, dado por las ecuaciones 1.2.44a-

1.2.44j, fue resuelto usando la rutina fsolve del modulo de optimización de
Matlab. La rutina fsolve es un programa para la solución numerica de sistemas
de ecuaciones algebraicas simultaneas no-lineales, el cual emplea un método de
optimización de minimos cuadrados no-lineales para resolver el problema; es de-
cir, el problema se resuelve como un problema de optimización. El programa
requiere de una estimación inicial xo de las variables desconocidas. Por ejemplo,
para el primer caso resuelto

xo=[1,1,1,1,1,1,1,1,1,1]’

la rutina se invoca, desde Matlab, de la siguiente forma

x=fsolve(’ecs’,xo)

donde ecs es el nombre de la rutina donde el usuario define las ecuaciones a

resolver, x es el vector de resultados. El programa normalmente funciona con una
serie de valores nominales, en caso de que el usuario no desee valores diferentes,
tales como el numero maximo de iteracciones, la precision requerida, etc. 2 A
continuación se muestra el listado de la rutina ecs para el primer caso resuelto
anteriormente.

function [f]=ecs(x)
a=[
-0.1300000E+02 0.1478830E+02 -0.2666666E+01 0.1878361E+01 -0.1000000E+01;
-0.5323791E+01 0.3872984E+01 0.2065591E+01 -0.1290995E+01 0.6762101E+00;
0.1500000E+01 -0.3227486E+01 0.0000000E+00 0.3227487E+01 -0.1500000E+01;
-0.6762100E+00 0.1290994E+01 -0.2065591E+01 -0.3872984E+01 0.5323791E+01;
0.1000000E+01 -0.1878361E+01 0.2666666E+01 -0.1478831E+02 0.1300000E+02];
b=[
0.8400001E+02 -0.1220632E+03 0.5866666E+02 -0.4460350E+02 0.2400000E+02;
0.5323790E+02 -0.7333334E+02 0.2666667E+02 -0.1333334E+02 0.6762102E+01;
-0.6000001E+01 0.1666667E+02 -0.2133333E+02 0.1666667E+02 -0.6000001E+01;
0.6762101E+01 -0.1333333E+02 0.2666666E+02 -0.7333334E+02 0.5323790E+02;
0.2400000E+02 -0.4460350E+02 0.5866666E+02 -0.1220632E+03 0.8400001E+02];

c1=x(1);
c2=x(2);
c3=x(3);
c4=x(4);
c5=x(5);
t1=x(6);
t2=x(7);
t3=x(8);
t4=x(9);
t5=x(10);

pe=2;
2
desde Matlab teclee help foptions para obtener mayor información sobre dichos parametros.
14CHAPTER 1. ECUACIONES DIFERENCIALES ORDINARIAS: PROBLEMAS DE VALORES EN

gama=17.6;
beta=-.056;
k=3.36;

f(1)= a(1,1)*c1+a(1,2)*c2+a(1,3)*c3+a(1,4)*c4+a(1,5)*c5-pe*(c1-1)
f(2)= b(2,1)*c1+b(2,2)*c2+b(2,3)*c3+b(2,4)*c4+b(2,5)*c5-pe*(a(2,1)*c1...
+a(2,2)*c2+a(2,3)*c3+a(2,4)*c4+a(2,5)*c5)-pe*k*c2*c2*exp(gama-(gama/t2))
f(3)= b(3,1)*c1+b(3,2)*c2+b(3,3)*c3+b(3,4)*c4+b(3,5)*c5-pe*(a(3,1)*c1...
+a(3,2)*c2+a(3,3)*c3+a(3,4)*c4+a(3,5)*c5)-pe*k*c3*c3*exp(gama-(gama/t3))
f(4)= b(4,1)*c1+b(4,2)*c2+b(4,3)*c3+b(4,4)*c4+b(4,5)*c5-pe*(a(4,1)*c1...
+a(4,2)*c2+a(4,3)*c3+a(4,4)*c4+a(4,5)*c5)-pe*k*c4*c4*exp(gama-(gama/t4))
f(5)= a(5,1)*c1+a(5,2)*c2+a(5,3)*c3+a(5,4)*c4+a(5,5)*c5
f(6)= a(1,1)*t1+a(1,2)*t2+a(1,3)*t3+a(1,4)*t4+a(1,5)*t5-pe*(t1-1)
f(7)= b(2,1)*t1+b(2,2)*t2+b(2,3)*t3+b(2,4)*t4+b(2,5)*t5-pe*(a(2,1)*t1...
+a(2,2)*t2+a(2,3)*t3+a(2,4)*t4+a(2,5)*t5)-pe*beta*k*c2*c2*exp(gama-(gama/t2))
f(8)= b(3,1)*t1+b(3,2)*t2+b(3,3)*t3+b(3,4)*t4+b(3,5)*t5-pe*(a(3,1)*t1...
+a(3,2)*t2+a(3,3)*t3+a(3,4)*t4+a(3,5)*t5)-pe*beta*k*c3*c3*exp(gama-(gama/t3))
f(9)= b(4,1)*t1+b(4,2)*t2+b(4,3)*t3+b(4,4)*t4+b(4,5)*t5-pe*(a(4,1)*t1...
+a(4,2)*t2+a(4,3)*t3+a(4,4)*t4+a(4,5)*t5)-pe*beta*k*c4*c4*exp(gama-(gama/t4))
f(10)= a(5,1)*t1+a(5,2)*t2+a(5,3)*t3+a(5,4)*t4+a(5,5)*t5

Ejemplo 5 Resolver el siguiente problema de valores a la frontera:

3py
y 00 + =0
(p + t2 )2

en el intérvalo t ∈ [−0.1, 0.1]. Las condiciones frontera dadas por

−0.1
y(−0.1) = √
p + 0.01
0.1
y(0.1) = √
p + 0.01

usar p =10−5 . Comparar la solución numérica contra la solución analı́tica dada por
t
y(t) = p
p + t2

Antes de proceder a aplicar el método de colocación ortogonal al problema en cuestión

resulta conveniente escalar la variable independiente t para que, la nueva variable escal-
ada, esté comprendida en el intérvalo [0,1]. Esto se puede hacer de manera muy sencilla
definiendo al nuevo tiempo escalado τ como,
t + 0.1
τ=
0.2
de esta forma cuando t = −0.1 tendremos τ = 0 y cuando t = 0.1 entonces τ = 1.
El empleo de este escalamiento implica que la ecuación diferencial original a resolver debe
reescribirse para que la variable independiente sea ahora τ . Entonces de la ecuación anterior,

t = 0.2τ − 0.1
Colocación Ortogonal. 15

por lo tanto
d2 τ = 0.2d2 t
reemplazando los términos t y dt2 en la ecuación diferencial original y rearreglando un poco
tenemos que la ecuación diferencial a resolver, con la nueva variable independiente τ , está
dada por,
d2 y 0.6py
2
+ =0
dτ [p + (0.2τ − 0.1)2 ]2
En la Figura 1.3 se muestran el perfil de variación del estado del sistema y con respecto
al tiempo evaluado empleando 10 puntos internos de colocación. Los listados Matlab para
resolver este ejemplo se muestran a continuación.

clear all; clc;

global npc bmat p tau

npc = input(’Number of internal collocation points ==> ?’);

[amat,bmat,qinv,roots,wx] = planar(npc);
p = 1e-5;
tau = roots;

y0 = -ones(npc+2,1);
options = optimset(’display’,’iter’);

x = fsolve(’problem4’,y0,options);
time = 0.2*tau-0.1;

tanalytic = linspace(-0.1,0.1);
xanalytic =tanalytic./sqrt(p+tanalytic.^2);

plot(tanalytic,xanalytic,’r-’), hold plot(time,x,’o-’)

legend (’Analytic solution’,’OC approximation’,0)

%-- End of the run.m file --

function fx = problem4 (y) global npc bmat p tau

%
% Approximation of the model within the domain solution
% (since the boundary conditions are first kind BCs, no
% discretization of the BCs is needed)
%
k = 1;
fx (k) = sqrt(p+0.01)*y(1)+0.1 ; % left boundary condition
for i = 2:npc+1,
sum = 0 ;
k = k+1 ;
16CHAPTER 1. ECUACIONES DIFERENCIALES ORDINARIAS: PROBLEMAS DE VALORES EN

for j = 1:npc+2,
sum = sum+bmat(i,j)*y(j);
end
fx (k) = sum+0.6*p*y(i)/(p+(0.2*tau(i)-0.1)^2)^2 ;
end
fx (k+1) = sqrt(p+0.01)*y(k+1)-0.1 ; % right boundary condition

%-- End of the problem4.m file --

0.8
Analytic solution
OC approximation
0.6

0.4

0.2

−0.2

−0.4

−0.6

−0.8

−1
−0.1 −0.08 −0.06 −0.04 −0.02 0 0.02 0.04 0.06 0.08 0.1

Figure 1.3: Comparación entre la solución analı́tica y la aproximación numérica por

colocación ortogonal para el ejemplo 5. Se emplearon 10 puntos internos de colocación.
El sı́mbolo “◦” indica los puntos donde el algoritmo de colocación ortogonal evalua la
solución.

1.3 Aplicación a un sistema distribuido

A continuación considereramos la manera de emplear el método de colocación ortogonal
cuando el modelo matemático del sistema en cuestión está planteado en términos de
ecuaciones diferenciales parciales.

Ejemplo 6 Considere el mismo sistema de reacción catalı́tico en un reactor tubular

con dispersión axial descrito anteriormente en el ejemplo 4. Sin embargo en este caso
nos interesa analizar la respuesta dinámica del reactor. Por lo tanto en esta situación el
sistema distribuido estará expresado en términos del siguiente conjunto de ecuaciones
diferenciales parciales no lineales.

∂c ∂ 2c ∂c
= 2 − P eM − P eM R(c, T ) (1.3.47)
∂t ∂z ∂z
Aplicación a un sistema distribuido 17

∂T ∂ 2T ∂T
= 2
− P eH − P eH βR(c, T ) (1.3.48)
∂t ∂z ∂z
sujeto a las condiciones iniciales:
c(z, 0) = 1 (1.3.49)
T (z, 0) = 1 (1.3.50)
y frontera:
∂c
= P eM (c − 1), en z = 0 (1.3.51)
∂z
∂c
= 0, en z = 1 (1.3.52)
∂z
∂T
= P eH (T − 1), en z = 0 (1.3.53)
∂z
∂T
= 0, en z = 1 (1.3.54)
∂z
(1.3.55)
donde, γ
R = kc2 eγ− T (1.3.56)
Si de nueva cuenta usamos tres puntos de colocación internos tendremos 10 incog-
nitas: c1 , c2 , c3 , c4 , c5 , T1 , T2 , T3 , T4 y T5 (ver figura 1.2) . Esto significa que necesitamos
plantear 10 ecuaciones para determinar el valor de dichas incognitas. En este caso
el modelo matemático sólo debe discretizarse sobre los puntos de colocación internos.
Esto nos dará el siguiente total de 6 ecuaciones diferenciales ordinarias discretizadas.

dc2
= B21 c1 + B22 c2 + B23 c3 + B24 c4 + B25 c5 − P eM [A21 c1 + A22 c2 + A23 c3 + A24 c4 + A25 c5 ]
dt
γ− γ
−P eM kc22 e T2
dc3
= B31 c1 + B32 c2 + B33 c3 + B34 c4 + B35 c5 − P eM [A31 c1 + A32 c2 + A33 c3 + A34 c4 + A35 c5 ]
dt
γ− γ
−P eM kc23 e T3
dc4
= B41 c1 + B42 c2 + B43 c3 + B44 c4 + B45 c5 − P eM [A41 c1 + A42 c2 + A43 c3 + A44 c4 + A45 c5 ]
dt
γ− γ
−P eM kc24 e T4
dT2
= B21 T1 + B22 T2 + B23 T3 + B24 T4 + B25 T5 − P eH [A21 T1 + A22 T2 + A23 T3 + A24 T4 + A25 T5 ]
dt
γ− γ
−P eH βkc22 e T2
dT3
= B31 T1 + B32 T2 + B33 T3 + B34 T4 + B35 T5 − P eH [A31 T1 + A32 T2 + A33 T3 + A34 T4 + A35 T5 ]
dt
γ− γ
−P eH βkc23 e T3
dT4
= B41 T1 + B42 T2 + B43 T3 + B44 T4 + B45 T5 − P eH [A41 T1 + A42 T2 + A43 T3 + A44 T4 + A45 T5 ]
dt
γ− γ
−P eH βkc24 e T4
18CHAPTER 1. ECUACIONES DIFERENCIALES ORDINARIAS: PROBLEMAS DE VALORES EN

mientras que la discretización de cada una de las condiciones frontera nos dará el
restante conjunto de 4 ecuaciones.

A11 c1 + A12 c2 + A13 c3 + A14 c4 + A15 c5 − P eM (c1 − 1) = 0

A51 c1 + A52 c2 + A53 c3 + A54 c4 + A55 c5 = 0
A11 T1 + A12 T2 + A13 T3 + A14 T4 + A15 T5 − P eH (T1 − 1) = 0
A51 T1 + A52 T2 + A53 T3 + A54 T4 + A55 T5 = 0

Para el conjunto de parámetros del caso 1 la respuesta dinámica del sistema, desde que
el conjunto de condiciones iniciales hasta que el reactor alcanza el estado estacionario se
muestra en la figuras 1.4 y 1.5. De manera semejante en las figuras 1.6 y 1.7 se muestra
la correspondiente respuesta dinámica para los parámetros del caso 2. A continuación
se muestra el programa Matlab empleado para realizar los cálculos reportados en dichas
figuras.

1
c2
c3
0.9 c4

0.8
Dimensionless Concentration

0.7

0.6

0.5

0.4

0.3

0.2
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Dimensionless time

Figure 1.4: Caso 1: Respuesta dinámica del reactor catalitico.

clear all; clc;

%
% Main driver program for performing the dynamic simulation of a
% catalytic tubular reactor. The PDE system is discretized using the
% orthogonal collocation approach.
%
Aplicación a un sistema distribuido 19

1.045

1.04

1.035
Dimensionless Temperature

1.03

1.025

1.02

1.015

1.01

T2
1.005
T3
T4
1
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Dimensionless time

Figure 1.5: Caso 1: Respuesta dinámica del reactor catalitico.

% Written by Antonio Flores T./ 6 Nov,2003

%
global pe gama beta k global a b
%
% Parameter values
%
pe = 2;
gama = 17.6;
beta = -.056;
k = 3.36;
%
% The collocation matrices
%
a=[
-0.1300000E+02 0.1478830E+02 -0.2666666E+01 0.1878361E+01 -0.1000000E+01;
-0.5323791E+01 0.3872984E+01 0.2065591E+01 -0.1290995E+01 0.6762101E+00;
0.1500000E+01 -0.3227486E+01 0.0000000E+00 0.3227487E+01 -0.1500000E+01;
-0.6762100E+00 0.1290994E+01 -0.2065591E+01 -0.3872984E+01 0.5323791E+01;
0.1000000E+01 -0.1878361E+01 0.2666666E+01 -0.1478831E+02 0.1300000E+02];
b=[
0.8400001E+02 -0.1220632E+03 0.5866666E+02 -0.4460350E+02 0.2400000E+02;
0.5323790E+02 -0.7333334E+02 0.2666667E+02 -0.1333334E+02 0.6762102E+01;
20CHAPTER 1. ECUACIONES DIFERENCIALES ORDINARIAS: PROBLEMAS DE VALORES EN

1
c2
c
3
0.9 c
4

0.8
Dimensionless Concentration

0.7

0.6

0.5

0.4

0.3

0.2

0.1
0 0.002 0.004 0.006 0.008 0.01 0.012 0.014 0.016 0.018 0.02
Dimensionless time

Figure 1.6: Caso 2: Respuesta dinámica del reactor catalitico.

-0.6000001E+01 0.1666667E+02 -0.2133333E+02 0.1666667E+02 -0.6000001E+01;

0.6762101E+01 -0.1333333E+02 0.2666666E+02 -0.7333334E+02 0.5323790E+02;
0.2400000E+02 -0.4460350E+02 0.5866666E+02 -0.1220632E+03 0.8400001E+02];
%
% Define initial conditions
%
xinitial = [1 1 1 1 1 1];
tspan = linspace(0,1,40);
[time,x] = ode15s(’ex5ode’,tspan,xinitial);
figure(1)
plot(time,x(:,1:3)), xlabel(’Dimensionless
time’),ylabel(’Dimensionless Concentration’)
legend(’c_{2}’,’c_{3}’,’c_{4}’,0)
figure(2)
plot(time,x(:,4:6)),xlabel(’Dimensionless
time’),ylabel(’Dimensionless Temperature’)
legend(’T_{2}’,’T_{3}’,’T_{4}’,0)

%-- End of the ex5.m file --

function [f]=ex5ode(time,x)
global pe gama beta k
Aplicación a un sistema distribuido 21

1.06
T2
T3
T4
1.05

Dimensionless Temperature

1.04

1.03

1.02

1.01

1
0 0.002 0.004 0.006 0.008 0.01 0.012 0.014 0.016 0.018 0.02
Dimensionless time

Figure 1.7: Caso 2: Respuesta dinámica del reactor catalitico.

global a b
global c2 c3 c4 t2 t3 t4
%
% Get the states values
%
c2=x(1);
c3=x(2);
c4=x(3);
t2=x(4);
t3=x(5);
t4=x(6);
%
% Solve the associated algebraic equations
%
options = optimset(’MaxIter’,200);
x0 = [1 1 1 1];

[xa,fval,flag] = fsolve(’ex5ae’,x0,options);
if flag <= 0
display(’ Warning: Convergence problems...’);
end
c1 = xa(1);
22CHAPTER 1. ECUACIONES DIFERENCIALES ORDINARIAS: PROBLEMAS DE VALORES EN

c5 = xa(2);
t1 = xa(3);
t5 = xa(4);

f(1)=b(2,1)*c1+b(2,2)*c2+b(2,3)*c3+b(2,4)*c4+b(2,5)*c5-pe*(a(2,1)*c1...
+a(2,2)*c2+a(2,3)*c3+a(2,4)*c4+a(2,5)*c5)-pe*k*c2*c2*exp(gama-(gama/t2));
f(2)=b(3,1)*c1+b(3,2)*c2+b(3,3)*c3+b(3,4)*c4+b(3,5)*c5-pe*(a(3,1)*c1...
+a(3,2)*c2+a(3,3)*c3+a(3,4)*c4+a(3,5)*c5)-pe*k*c3*c3*exp(gama-(gama/t3));
f(3)=b(4,1)*c1+b(4,2)*c2+b(4,3)*c3+b(4,4)*c4+b(4,5)*c5-pe*(a(4,1)*c1...
+a(4,2)*c2+a(4,3)*c3+a(4,4)*c4+a(4,5)*c5)-pe*k*c4*c4*exp(gama-(gama/t4));
f(4)=b(2,1)*t1+b(2,2)*t2+b(2,3)*t3+b(2,4)*t4+b(2,5)*t5-pe*(a(2,1)*t1...
+a(2,2)*t2+a(2,3)*t3+a(2,4)*t4+a(2,5)*t5)-pe*beta*k*c2*c2*exp(gama-(gama/t2));
f(5)=b(3,1)*t1+b(3,2)*t2+b(3,3)*t3+b(3,4)*t4+b(3,5)*t5-pe*(a(3,1)*t1...
+a(3,2)*t2+a(3,3)*t3+a(3,4)*t4+a(3,5)*t5)-pe*beta*k*c3*c3*exp(gama-(gama/t3));
f(6)=b(4,1)*t1+b(4,2)*t2+b(4,3)*t3+b(4,4)*t4+b(4,5)*t5-pe*(a(4,1)*t1...
+a(4,2)*t2+a(4,3)*t3+a(4,4)*t4+a(4,5)*t5)-pe*beta*k*c4*c4*exp(gama-(gama/t4));

f = f’;

%-- End of the ex5ode.m file --

function [f] = ex5ae(x)

global pe gama beta k
global a b
global c2 c3 c4 t2 t3 t4
%
% This file computes the values of both concentration and temperature
% on the system boundaries.
%
% Get the value of the algebraic variables
%
c1=x(1);
c5=x(2);
t1=x(3);
t5=x(4);

f(1)= a(1,1)*c1+a(1,2)*c2+a(1,3)*c3+a(1,4)*c4+a(1,5)*c5-pe*(c1-1);
f(2)= a(5,1)*c1+a(5,2)*c2+a(5,3)*c3+a(5,4)*c4+a(5,5)*c5;
f(3)=a(1,1)*t1+a(1,2)*t2+a(1,3)*t3+a(1,4)*t4+a(1,5)*t5-pe*(t1-1);
f(4) = a(5,1)*t1+a(5,2)*t2+a(5,3)*t3+a(5,4)*t4+a(5,5)*t5;

f = f’;

%-- End of the ex5ae.m file --

Problemas en varias dimensiones 23

1.4 Problemas en varias dimensiones

El método de colocación ortogonal también puede ser aplicado a sistemas distribuidos
en dos y tres dimensiones espaciales, tanto en cooordenadas cartesianas, cilindrı́cas y
esféricas.
Supongamos que nuestro modelo de parámetros distribuidos, en dos dimensiones
x, y, se puede representar por la siguiente ecuación

f (u, x, y, t, ut , ux , uxx , uy , uyy ) = 0

donde u representa a la variable dependiente. La discretización de las derivadas par-

ciales, a lo largo de cada una de las dimensiones, está dada por el siguiente conjunto
de ecuaciones.
¯ Nx +2
∂ 2 u ¯¯ X
x
= Bjk uik (1.4.57)
∂x2 ¯ij k=1
¯ N x +2
∂u ¯¯ X
= Axjk uik (1.4.58)
∂x ¯ij k=1
¯ Ny +2
∂ 2 u ¯¯ X y
= Bik ukj (1.4.59)
∂y 2 ¯ij k=1
¯ Ny +2
X
∂u ¯¯
= Ayik ukj (1.4.60)
∂y ¯ij k=1

donde i = 1, ..., Nx + 2, j = 1, ..., Ny + 2. Nx y Ny representan los puntos de colocación

internos (sin incluir las fronteras) a lo largo de las coordenadas x e y, respectivamente.
Nótese que Ax y B x representan las matrices de colocación de las primeras y segundas
derivadas parciales a lo largo de la dimensión x. De manera semejante Ay y B y repre-
sentan las correspondientes matrices de colocación de primeras y segundas derivadas
parciales a lo largo de la dimensión y. Las matrices de discretización a lo largo de cada
eje sólo serán iguales entre si el número de puntos de colocación sobre el eje x es igual
al número de puntos de colocación sobre el eje y.

Ejemplo 7 Considere el problema dinámico de transferencia de calor en dos dimensiones

tal como se ilustra en la figura 1.8. Se supone que toda la lamina se encuentra inicialmente a
la temperatura To y que dos de los extremos se aislan mientras que los otros dos restantes
extremos se encuentran expuestos al medio ambiente. Debido a las caracterı́sticas de
simetrı́a del problema en cuestión sólo se considerara una porción de la barra. El modelo
dinámico bidimensional de transferencia de calor por conducción está dado por la siguiente
ecuación diferencial parcial lineal
µ ¶
∂T ∂ 2T ∂2T
=D +
∂τ ∂x2 ∂y 2
24CHAPTER 1. ECUACIONES DIFERENCIALES ORDINARIAS: PROBLEMAS DE VALORES EN

y
6

¢¢
B ¢
¢¢
¢¢ -
¢¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢
x

A A

Figure 1.8: Transferencia de calor por conducción, en dos dimensiones, en una barra.

sujeto a las siguientes condiciones iniciales y frontera:

T (x, y, 0) = To
∂T (0, y, τ )
= 0
∂x
∂T (x, 0, τ )
= 0
∂y
µ ¶
∂T (A, y, τ )
−k = hT
∂x
µ ¶
∂T (x, B, τ )
−k = hT
∂y

en términos de las siguientes variables adimensionales

x y T Dτ
x̂ = , ŷ = , θ = , t = 2
A B To A
donde h = 1,k = 1,A = 1 y B = 1. Tal como se muestra en la figura 1.9 en este caso
usaremos 3 puntos colocación internos sobre el eje x y dos puntos de colocación internos
sobre el eje y.
El modelo matemático original, junto con sus correspondientes condiciones iniciales y
frontera, se puede reescribir en términos de las nuevas variables adimensionales como se
Problemas en varias dimensiones 25

1 y4 u u u u u

y3 u u u u u

y2 u u u u u

0 y1 u u u u u

x1 x2 x3 x4 x5
0 1

Figure 1.9: Aplicación del método de colocación ortogonal en dos dimensiones .

muestra a continuación.

∂θ ∂2θ 2
2∂ θ
= + L
∂t ∂ x̂2 ∂ ŷ 2
θ(x̂, ŷ, 0) = 1
∂θ(0, ŷ, t)
= 0
∂ x̂
∂θ(x̂, 0, t)
= 0
∂ ŷ
∂θ(1, ŷ, t) hAθ
− =
∂ x̂ k
∂θ(x̂, 1, t) hBθ
− =
∂ ŷ k

A
donde L = B
es la relación de aspecto.
26CHAPTER 1. ECUACIONES DIFERENCIALES ORDINARIAS: PROBLEMAS DE VALORES EN

1. Aproximación del modelo matemático

dθ11 x x x x x y y y y
= B11 θ11 + B12 θ12 + B13 θ13 + B14 θ14 + B15 θ15 + L2 (B11 θ11 + B12 θ21 + B13 θ31 + B14 θ41 )
dt
dθ12 x x x x x y y y y
= B21 θ11 + B22 θ12 + B23 θ13 + B24 θ14 + B25 θ15 + L2 (B11 θ12 + B12 θ22 + B13 θ32 + B14 θ42 )
dt
dθ13 x x x x x y y y y
= B31 θ11 + B32 θ12 + B33 θ13 + B34 θ14 + B35 θ15 + L2 (B11 θ13 + B12 θ23 + B13 θ33 + B14 θ43 )
dt
dθ14 x x x x x y y y y
= B41 θ11 + B42 θ12 + B43 θ13 + B44 θ14 + B45 θ15 + L2 (B11 θ14 + B12 θ24 + B13 θ34 + B14 θ44 )
dt
dθ15 x x x x x y y y y
= B51 θ11 + B52 θ12 + B53 θ13 + B54 θ14 + B55 θ15 + L2 (B11 θ15 + B12 θ25 + B13 θ35 + B14 θ45 )
dt
dθ21 x x x x x y y y y
= B11 θ21 + B12 θ22 + B13 θ23 + B14 θ24 + B15 θ25 + L2 (B21 θ11 + B22 θ21 + B23 θ31 + B24 θ41 )
dt
dθ22 x x x x x y y y y
= B21 θ21 + B22 θ22 + B23 θ23 + B24 θ24 + B25 θ25 + L2 (B21 θ12 + B22 θ22 + B23 θ32 + B24 θ42 )
dt
dθ23 x x x x x y y y y
= B31 θ21 + B32 θ22 + B33 θ23 + B34 θ24 + B35 θ25 + L2 (B21 θ13 + B22 θ23 + B23 θ33 + B24 θ43 )
dt
dθ24 x x x x x y y y y
= B41 θ21 + B42 θ22 + B43 θ23 + B44 θ24 + B45 θ25 + L2 (B21 θ14 + B22 θ24 + B23 θ34 + B24 θ44 )
dt
dθ25 x x x x x y y y y
= B51 θ21 + B52 θ22 + B53 θ23 + B54 θ24 + B55 θ25 + L2 (B21 θ15 + B22 θ25 + B23 θ35 + B24 θ45 )
dt
dθ31 x x x x x y y y y
= B11 θ31 + B12 θ32 + B13 θ33 + B14 θ34 + B15 θ35 + L2 (B31 θ11 + B32 θ21 + B33 θ31 + B34 θ41 )
dt
dθ32 x x x x x y y y y
= B21 θ31 + B22 θ32 + B23 θ33 + B24 θ34 + B25 θ35 + L2 (B31 θ12 + B32 θ22 + B33 θ32 + B34 θ42 )
dt
dθ33 x x x x x y y y y
= B31 θ31 + B32 θ32 + B33 θ33 + B34 θ34 + B35 θ35 + L2 (B31 θ13 + B32 θ23 + B33 θ33 + B34 θ43 )
dt
dθ34 x x x x x y y y y
= B41 θ31 + B42 θ32 + B43 θ33 + B44 θ34 + B45 θ35 + L2 (B31 θ14 + B32 θ24 + B33 θ34 + B34 θ44 )
dt
dθ35 x x x x x y y y y
= B51 θ31 + B52 θ32 + B53 θ33 + B54 θ34 + B55 θ35 + L2 (B31 θ15 + B32 θ25 + B33 θ35 + B34 θ45 )
dt
dθ41 x x x x x y y y y
= B11 θ41 + B12 θ42 + B13 θ43 + B14 θ44 + B15 θ45 + L2 (B41 θ11 + B42 θ21 + B43 θ31 + B44 θ41 )
dt
dθ42 x x x x x y y y y
= B21 θ41 + B22 θ42 + B23 θ43 + B24 θ44 + B25 θ45 + L2 (B41 θ12 + B42 θ22 + B43 θ32 + B44 θ42 )
dt
dθ43 x x x x x y y y y
= B31 θ41 + B32 θ42 + B33 θ43 + B34 θ44 + B35 θ45 + L2 (B41 θ13 + B42 θ23 + B43 θ33 + B44 θ43 )
dt
dθ44 x x x x x y y y y
= B41 θ41 + B42 θ42 + B43 θ43 + B44 θ44 + B45 θ45 + L2 (B41 θ14 + B42 θ24 + B43 θ34 + B44 θ44 )
dt
dθ45 x x x x x y y y y
= B51 θ41 + B52 θ42 + B53 θ43 + B54 θ44 + B55 θ45 + L2 (B41 θ15 + B42 θ25 + B43 θ35 + B44 θ45 )
dt

2. Condiciones frontera
Problemas en varias dimensiones 27
¯
∂θ ¯
(a) ∂x x=0
=0

Ax11 θ11 + Ax12 θ12 + Ax13 θ13 + Ax14 θ14 + Ax15 θ15 = 0
Ax11 θ21 + Ax12 θ22 + Ax13 θ23 + Ax14 θ24 + Ax15 θ25 = 0
Ax11 θ31 + Ax12 θ32 + Ax13 θ33 + Ax14 θ34 + Ax15 θ35 = 0
Ax11 θ41 + Ax12 θ42 + Ax13 θ43 + Ax14 θ44 + Ax15 θ45 = 0
¯
∂θ ¯ hAθ
(b) − ∂x x=1
= k

hAθ15
Ax51 θ11 + Ax52 θ12 + Ax53 θ13 + Ax54 θ14 + Ax55 θ15 = −
k
hAθ25
Ax51 θ21 + Ax52 θ22 + Ax53 θ23 + Ax54 θ24 + Ax55 θ25 = −
k
hAθ35
Ax51 θ31 + Ax52 θ32 + Ax53 θ33 + Ax54 θ34 + Ax55 θ35 = −
k
hAθ45
Ax51 θ41 + Ax52 θ42 + Ax53 θ43 + Ax54 θ44 + Ax55 θ45 = −
k
¯
∂θ ¯
(c) ∂y ¯
=0
y=0

Ay11 θ11 + Ay12 θ21 + Ay13 θ31 + Ay14 θ41 = 0

Ay11 θ12 + Ay12 θ22 + Ay13 θ32 + Ay14 θ42 = 0
Ay11 θ13 + Ay12 θ23 + Ay13 θ33 + Ay14 θ43 = 0
Ay11 θ14 + Ay12 θ24 + Ay13 θ34 + Ay14 θ44 = 0
Ay11 θ15 + Ay12 θ25 + Ay13 θ35 + Ay14 θ45 = 0
¯
∂θ ¯ hBθ
(d) − ∂y ¯ = k
y=1

hBθ41
Ay41 θ11 + Ay42 θ21 + Ay43 θ31 + Ay44 θ41 = −
k
hBθ42
Ay41 θ12 + Ay42 θ22 + Ay43 θ32 + Ay44 θ42 = −
k
hBθ43
Ay41 θ13 + Ay42 θ23 + Ay43 θ33 + Ay44 θ43 = −
k
hBθ44
Ay41 θ14 + Ay42 θ24 + Ay43 θ34 + Ay44 θ44 = −
k
hBθ45
Ay41 θ15 + Ay42 θ25 + Ay43 θ35 + Ay44 θ45 = −
k
Como puede observarse el número total de ecuaciones resultantes de la discretización (mod-
elo y condiciones frontera) es de 38 pero sólo tenemos 20 incognitas (todas las temperaturas
en cada uno de los puntos de colocación de la figura 1.9). En realidad el modelo matemático
sólo debe aproximarse en los puntos de colocación internos, mientras que las condiciones
frontera se usan para definir el valor de las temperatura en los puntos de discretización
sobre las fronteras del sistema.
28CHAPTER 1. ECUACIONES DIFERENCIALES ORDINARIAS: PROBLEMAS DE VALORES EN

Debe notarse también que algunas de las condiciones frontera han sido aproximadas dos
veces. Esto es cierto en el caso de los puntos de discretización que están colocados sobre
los vertices del dominio de solución de la figura 1.9 (θ11 , θ15 , θ41 ¯y θ45 ). En¯ este problema
arbitrariamente decidimos aproximar las condiciones frontera ∂∂θx̂ ¯x̂=0 y ∂∂θx̂ ¯x̂=1 empleando
todos los puntos
¯ de colocación sobre la coordenada y, pero aproximar las condiciones
∂θ ¯
¯
frontera ∂ ŷ ¯ = y ∂∂θx̂ ¯ŷ=1 sólo sobre los puntos internos de colocación definidos sobre la
ŷ=0
coordenada x. De esta forma se obtiene el mismo número de ecuaciones y de incognitas.

Procediendo de la manera explicada anteriormente el sistema de ecuaciones discretizadas

que modela la dinámica de transferencia de calor en dos dimensiones estará dado como se
muestra a continuación.

dθ22 x x x x x y y y y
= B21 θ21 + B22 θ22 + B23 θ23 + B24 θ24 + B25 θ25 + L2 (B21 θ12 + B22 θ22 + B23 θ32 + B24 θ42 )
dt
dθ23 x x x x x y y y y
= B31 θ21 + B32 θ22 + B33 θ23 + B34 θ24 + B35 θ25 + L2 (B21 θ13 + B22 θ23 + B23 θ33 + B24 θ43 )
dt
dθ24 x x x x x y y y y
= B41 θ21 + B42 θ22 + B43 θ23 + B44 θ24 + B45 θ25 + L2 (B21 θ14 + B22 θ24 + B23 θ34 + B24 θ44 )
dt
dθ32 x x x x x y y y y
= B21 θ31 + B22 θ32 + B23 θ33 + B24 θ34 + B25 θ35 + L2 (B31 θ12 + B32 θ22 + B33 θ32 + B34 θ42 )
dt
dθ33 x x x x x y y y y
= B31 θ31 + B32 θ32 + B33 θ33 + B34 θ34 + B35 θ35 + L2 (B31 θ13 + B32 θ23 + B33 θ33 + B34 θ43 )
dt
dθ34 x x x x x y y y y
= B41 θ31 + B42 θ32 + B43 θ33 + B44 θ34 + B45 θ35 + L2 (B31 θ14 + B32 θ24 + B33 θ34 + B34 θ44 )
dt
Problemas en varias dimensiones 29

sujetas a las siguientes condiciones frontera

Ax11 θ11 + Ax12 θ12 + Ax13 θ13 + Ax14 θ14 + Ax15 θ15 = 0
Ax11 θ21 + Ax12 θ22 + Ax13 θ23 + Ax14 θ24 + Ax15 θ25 = 0
Ax11 θ31 + Ax12 θ32 + Ax13 θ33 + Ax14 θ34 + Ax15 θ35 = 0
Ax11 θ41 + Ax12 θ42 + Ax13 θ43 + Ax14 θ44 + Ax15 θ45 = 0
hAθ15
Ax51 θ11 + Ax52 θ12 + Ax53 θ13 + Ax54 θ14 + Ax55 θ15 = −
k
hAθ25
Ax51 θ21 + Ax52 θ22 + Ax53 θ23 + Ax54 θ24 + Ax55 θ25 = −
k
hAθ35
Ax51 θ31 + Ax52 θ32 + Ax53 θ33 + Ax54 θ34 + Ax55 θ35 = −
k
hAθ45
Ax51 θ41 + Ax52 θ42 + Ax53 θ43 + Ax54 θ44 + Ax55 θ45 = −
k
Ay11 θ12 + Ay12 θ22 + Ay13 θ32 + Ay14 θ42 = 0
Ay11 θ13 + Ay12 θ23 + Ay13 θ33 + Ay14 θ43 = 0
Ay11 θ14 + Ay12 θ24 + Ay13 θ34 + Ay14 θ44 = 0
hBθ42
Ay41 θ12 + Ay42 θ22 + Ay43 θ32 + Ay44 θ42 = −
k
hBθ43
Ay41 θ13 + Ay42 θ23 + Ay43 θ33 + Ay44 θ43 = −
k
hBθ44
Ay41 θ14 + Ay42 θ24 + Ay43 θ34 + Ay44 θ44 = −
k

En la figura 1.10 se muestran las gráficas de respuesta dinámica del proceso de transfer-
encia de calor por conducción en dos dimensiones para cada uno de los puntos de colocación
internos. Los resultados mostrados fueron obtenidos usando el siguiente programa Matlab.

clc; clear all; warning off;

%
% This file is the main driver program for performing the dynamic simulation of the
% conduction heat process in an initially heated bar. The heat conduction is assumed
% to occur along the X and Y cartesian dimensions. Orthogonal collocation is used to
% get a discretized model consisting of a set of odes.
%
%
% Made by Antonio Flores T./ 4 Nov,2003
%
global nx ny bx ax by ay
global l alpha beta

nix = 3; % number of internal collocation points (x-axis)

niy = 2; % number of internal collocation points (y-axis)
k = 1; % thermal conductivity
h = 1; % heat transfer coefficient
a = 1; % bar length
30CHAPTER 1. ECUACIONES DIFERENCIALES ORDINARIAS: PROBLEMAS DE VALORES EN

1.4

1.2
θ
22
θ23
1 θ
24
θ
32
θ
33
θ
34
0.8
θ

0.6

0.4

0.2

0
0 0.5 1 1.5 2 2.5 3
Dimensionless time

Figure 1.10: Respuesta dinámica del proceso de transferencia de calor por conducción
en dos dimensiones.

b = 1; % bar width
nx = nix+2; % total number of points along x-axis
ny = niy+2; % total number of pints along y-axis
l = (a/b)^2; % squared aspect ratio
alpha = h*a/k;
beta = h*b/k;
%
% Collocation matrices along x-direction
%
ax = [
-0.130000010E+02 0.147883043E+02 -0.266666627E+01 0.187836111E+01 -0.100000000E+01
-0.532379103E+01 0.387298441E+01 0.206559086E+01 -0.129099452E+01 0.676210105E+00
0.150000024E+01 -0.322748613E+01 0.000000000E+00 0.322748661E+01 -0.150000024E+01
-0.676209986E+00 0.129099429E+01 -0.206559062E+01 -0.387298441E+01 0.532379055E+01
0.100000000E+01 -0.187836087E+01 0.266666627E+01 -0.147883062E+02 0.130000010E+02];
bx = [
0.840000076E+02 -0.122063164E+03 0.586666641E+02 -0.446035042E+02 0.240000019E+02
0.532379036E+02 -0.733333359E+02 0.266666679E+02 -0.133333368E+02 0.676210260E+01
-0.600000095E+01 0.166666660E+02 -0.213333321E+02 0.166666679E+02 -0.600000095E+01
0.676210117E+01 -0.133333340E+02 0.266666641E+02 -0.733333359E+02 0.532378998E+02
0.240000019E+02 -0.446034966E+02 0.586666641E+02 -0.122063187E+03 0.840000076E+02];
%
% Collocation matrices along y-direction
%
ay = [
Problemas en varias dimensiones 31

-0.700000000E+01 0.819615269E+01 -0.219615269E+01 0.100000000E+01

-0.273205090E+01 0.173205090E+01 0.173205101E+01 -0.732050836E+00
0.732050776E+00 -0.173205078E+01 -0.173205090E+01 0.273205066E+01
-0.100000000E+01 0.219615245E+01 -0.819615269E+01 0.700000000E+01];

by =[
0.240000000E+02 -0.371769142E+02 0.251769180E+02 -0.120000000E+02
0.163923054E+02 -0.240000019E+02 0.120000019E+02 -0.439230537E+01
-0.439230490E+01 0.120000000E+02 -0.240000000E+02 0.163923035E+02
-0.120000000E+02 0.251769142E+02 -0.371769142E+02 0.240000000E+02];
%
% Initial conditions
%
tspan = linspace(0,3);
xinitial = [1 1 1 1 1 1]’;
[time,x] = ode15s(’ex5ode’,tspan,xinitial);
plot(time,x), xlabel(’Dimensionless time’), ylabel(’\theta’)
legend(’\theta_{22}’,’\theta_{23}’,’\theta_{24}’,’\theta_{32}’,’\theta_{33}’,’\theta_{34}’)

%-- End of the ex5.m file --

function fx = ex5ode (t,x)

global nx ny bx ax by ay
global l alpha beta
global theta22 theta23 theta24 theta32 theta33 theta34
%
% Dynamic mathematical of the discretized 2-dimensional
% heat conduction conduction equation.
%
%
% Get states (only internal collocation points are true states)
%
theta22 = x(1);
theta23 = x(2);
theta24 = x(3);
theta32 = x(4);
theta33 = x(5);
theta34 = x(6);
%
% Solve for the variables defined on the system boundaries
%
x0 = [1 1 1 1 1 1 1 1 1 1 1 1 1 1];
[thetab] = fsolve(’ex5ae’,x0);
theta11 = thetab(1); theta12 = thetab(2); theta13 = thetab(3);
theta14 = thetab(4); theta15 = thetab(5); theta41 = thetab(6);
theta42 = thetab(7); theta43 = thetab(8); theta44 = thetab(9);
theta45 = thetab(10); theta21 = thetab(11); theta31 = thetab(12);
theta25 = thetab(13); theta35 = thetab(14);
%
% Discretized dynamic math model...
%
dt22=bx(2,1)*theta21+bx(2,2)*theta22+bx(2,3)*theta23+bx(2,4)*theta24+bx(2,5)*theta25+...
l*(by(2,1)*theta12+by(2,2)*theta22+by(2,3)*theta32+by(2,4)*theta42);
32CHAPTER 1. ECUACIONES DIFERENCIALES ORDINARIAS: PROBLEMAS DE VALORES EN

dt23=bx(3,1)*theta21+bx(3,2)*theta22+bx(3,3)*theta23+bx(3,4)*theta24+bx(3,5)*theta25+...
l*(by(2,1)*theta13+by(2,2)*theta23+by(2,3)*theta33+by(2,4)*theta43);
dt24=bx(4,1)*theta21+bx(4,2)*theta22+bx(4,3)*theta23+bx(4,4)*theta24+bx(4,5)*theta25+...
l*(by(2,1)*theta14+by(2,2)*theta24+by(2,3)*theta34+by(2,4)*theta44);
dt32=bx(2,1)*theta31+bx(2,2)*theta32+bx(2,3)*theta33+bx(2,4)*theta34+bx(2,5)*theta35+...
l*(by(3,1)*theta12+by(3,2)*theta22+by(3,3)*theta32+by(3,4)*theta42);
dt33=bx(3,1)*theta31+bx(3,2)*theta32+bx(3,3)*theta33+bx(3,4)*theta34+bx(3,5)*theta35+...
l*(by(3,1)*theta13+by(3,2)*theta23+by(3,3)*theta33+by(3,4)*theta43);
dt34=bx(4,1)*theta31+bx(4,2)*theta32+bx(4,3)*theta33+bx(4,4)*theta34+bx(4,5)*theta35+...
l*(by(3,1)*theta14+by(3,2)*theta24+by(3,3)*theta34+by(3,4)*theta44);

fx = [dt22 dt23 dt24 dt32 dt33 dt34]’;

%-- End of the ex5ode.m file --

function fx = ex5ae (x)

%
% This file computes the values of the temperatures defined on the boundaries
%
global nx ny bx ax by ay
global l alpha beta
global theta22 theta23 theta24 theta32 theta33 theta34
%
% Get the states
%
theta11 = x(1);
theta12 = x(2);
theta13 = x(3);
theta14 = x(4);
theta15 = x(5);
theta41 = x(6);
theta42 = x(7);
theta43 = x(8);
theta44 = x(9);
theta45 = x(10);
theta21 = x(11);
theta31 = x(12);
theta25 = x(13);
theta35 = x(14);

fx(1)=ax(1,1)*theta11+ax(1,2)*theta12+ax(1,3)*theta13+ax(1,4)*theta14+ax(1,5)*theta15;
fx(2)=ay(1,1)*theta12+ay(1,2)*theta22+ay(1,3)*theta32+ay(1,4)*theta42;
fx(3)=ay(1,1)*theta13+ay(1,2)*theta23+ay(1,3)*theta33+ay(1,4)*theta43;
fx(4)=ay(1,1)*theta14+ay(1,2)*theta24+ay(1,3)*theta34+ay(1,4)*theta44;
fx(5)=ax(5,1)*theta11+ax(5,2)*theta12+ax(5,3)*theta13+ax(5,4)*theta14+ax(5,5)*theta15;
fx(6)=ax(1,1)*theta41+ax(1,2)*theta42+ax(1,3)*theta43+ax(1,4)*theta44+ax(1,5)*theta45;
fx(7)=ay(4,1)*theta12+ay(4,2)*theta22+ay(4,3)*theta32+ay(4,4)*theta42+beta*theta42;
fx(8)=ay(4,1)*theta13+ay(4,2)*theta23+ay(4,3)*theta33+ay(4,4)*theta43+beta*theta43;
fx(9)=ay(4,1)*theta14+ay(4,2)*theta24+ay(4,3)*theta34+ay(4,4)*theta44+beta*theta44;
fx(10)=ax(5,1)*theta41+ax(5,2)*theta42+ax(5,3)*theta43+ax(5,4)*theta44+(ax(5,5)+alpha)...
*theta45;
fx(11)=ax(1,1)*theta21+ax(1,2)*theta22+ax(1,3)*theta23+ax(1,4)*theta24+ax(1,5)*theta25;
fx(12)=ax(1,1)*theta31+ax(1,2)*theta32+ax(1,3)*theta33+ax(1,4)*theta34+ax(1,5)*theta35;
fx(13)=ax(5,1)*theta21+ax(5,2)*theta22+ax(5,3)*theta23+ax(5,4)*theta24+(ax(5,5)+alpha)...
Problemas en varias dimensiones 33

*theta25;
fx(14)=ax(5,1)*theta31+ax(5,2)*theta32+ax(5,3)*theta33+ax(5,4)*theta34+(ax(5,5)+alpha)...
*theta35;

fx = fx’;

%-- End of the ex5ae.m file --

El problema discretizado conduce a un sistema simultaneo de ecuaciones diferenciales algebraicas el

cual puede ser integrado directamente por Matlab empleando la rutina ode15s tal como se muestra en el
siguiente listado.

function test_ex5
clc; clear all;
warning off;
%
% This file is the main driver program for performing the dynamic simulation of the
% conduction heat process in an initially heated bar. The heat conduction is assumed
% to occur along the X and Y cartesian dimensions. Orthogonal collocation is used to
% get a discretized model consisting of a set of odes.
%
%
% Made by Antonio Flores T./ 14 Nov,2003
%
global nx ny bx ax by ay
global l alpha beta

nix = 3; % number of internal collocation points (x-axis)

niy = 2; % number of internal collocation points (y-axis)
k = 1; % thermal conductivity
h = 1; % heat transfer coefficient
a = 1; % bar length
b = 1; % bar width
nx = nix+2; % total number of points along x-axis
ny = niy+2; % totakl number of pints along y-axis
l = (a/b)^2; % squared aspect ratio
alpha = h*a/k;
beta = h*b/k;
%
% Collocation matrices along x-direction
%
ax = [
-0.130000010E+02 0.147883043E+02 -0.266666627E+01 0.187836111E+01 -0.100000000E+01
-0.532379103E+01 0.387298441E+01 0.206559086E+01 -0.129099452E+01 0.676210105E+00
0.150000024E+01 -0.322748613E+01 0.000000000E+00 0.322748661E+01 -0.150000024E+01
-0.676209986E+00 0.129099429E+01 -0.206559062E+01 -0.387298441E+01 0.532379055E+01
0.100000000E+01 -0.187836087E+01 0.266666627E+01 -0.147883062E+02 0.130000010E+02];
bx = [
0.840000076E+02 -0.122063164E+03 0.586666641E+02 -0.446035042E+02 0.240000019E+02
0.532379036E+02 -0.733333359E+02 0.266666679E+02 -0.133333368E+02 0.676210260E+01
-0.600000095E+01 0.166666660E+02 -0.213333321E+02 0.166666679E+02 -0.600000095E+01
0.676210117E+01 -0.133333340E+02 0.266666641E+02 -0.733333359E+02 0.532378998E+02
34CHAPTER 1. ECUACIONES DIFERENCIALES ORDINARIAS: PROBLEMAS DE VALORES EN

0.240000019E+02 -0.446034966E+02 0.586666641E+02 -0.122063187E+03 0.840000076E+02];

%
% Collocation matrices along y-direction
%
ay = [
-0.700000000E+01 0.819615269E+01 -0.219615269E+01 0.100000000E+01
-0.273205090E+01 0.173205090E+01 0.173205101E+01 -0.732050836E+00
0.732050776E+00 -0.173205078E+01 -0.173205090E+01 0.273205066E+01
-0.100000000E+01 0.219615245E+01 -0.819615269E+01 0.700000000E+01];

by =[
0.240000000E+02 -0.371769142E+02 0.251769180E+02 -0.120000000E+02
0.163923054E+02 -0.240000019E+02 0.120000019E+02 -0.439230537E+01
-0.439230490E+01 0.120000000E+02 -0.240000000E+02 0.163923035E+02
-0.120000000E+02 0.251769142E+02 -0.371769142E+02 0.240000000E+02];
%
% Initial conditions
%
tspan = linspace(0,3);
xinitial = ones(20,1)’;
mass_matrix =zeros(20,20);
mass_matrix(1,1)=1;
mass_matrix(2,2)=1;
mass_matrix(3,3)=1;
mass_matrix(4,4)=1;
mass_matrix(5,5)=1;
mass_matrix(6,6)=1;
options = odeset(’Mass’,mass_matrix);
[time,x] = ode15s(@ex5ode,tspan,xinitial,options);
plot(time,x(:,1:6)) xlabel(’Dimensionless time’), ylabel(’\theta’)
legend(’\theta_{22}’,’\theta_{23}’,’\theta_{24}’,’\theta_{32}’,’\theta_{33}’,’\theta_{34}’)

%-- End of the test_ex5.m file --

function fx = ex5ode (t,x)

global nx ny bx ax by ay
global l alpha beta

% Dynamic mathematical of the discretized 2-dimensional

% heat conduction conduction equation.
%
%
% Get states (only internal collocation points are true states)
%
theta22 = x(1);
theta23 = x(2);
theta24 = x(3);
theta32 = x(4);
theta33 = x(5);
theta34 = x(6);
%
% Get the value of the algebraic variables
%
theta11 = x(7);
Problemas en varias dimensiones 35

theta12 = x(8);
theta13 = x(9);
theta14 = x(10);
theta15 = x(11);
theta41 = x(12);
theta42 = x(13);
theta43 = x(14);
theta44 = x(15);
theta45 = x(16);
theta21 = x(17);
theta31 = x(18);
theta25 = x(19);
theta35 = x(20);
%
% Discretized ordinary differential equations dynamic math model...
%
%
fx(1)=bx(2,1)*theta21+bx(2,2)*theta22+bx(2,3)*theta23+bx(2,4)*theta24+bx(2,5)*theta25+...
l*(by(2,1)*theta12+by(2,2)*theta22+by(2,3)*theta32+by(2,4)*theta42);
fx(2)=bx(3,1)*theta21+bx(3,2)*theta22+bx(3,3)*theta23+bx(3,4)*theta24+bx(3,5)*theta25+...
l*(by(2,1)*theta13+by(2,2)*theta23+by(2,3)*theta33+by(2,4)*theta43);
fx(3)=bx(4,1)*theta21+bx(4,2)*theta22+bx(4,3)*theta23+bx(4,4)*theta24+bx(4,5)*theta25+...
l*(by(2,1)*theta14+by(2,2)*theta24+by(2,3)*theta34+by(2,4)*theta44);
fx(4)=bx(2,1)*theta31+bx(2,2)*theta32+bx(2,3)*theta33+bx(2,4)*theta34+bx(2,5)*theta35+...
l*(by(3,1)*theta12+by(3,2)*theta22+by(3,3)*theta32+by(3,4)*theta42);
fx(5)=bx(3,1)*theta31+bx(3,2)*theta32+bx(3,3)*theta33+bx(3,4)*theta34+bx(3,5)*theta35+...
l*(by(3,1)*theta13+by(3,2)*theta23+by(3,3)*theta33+by(3,4)*theta43);
fx(6)=bx(4,1)*theta31+bx(4,2)*theta32+bx(4,3)*theta33+bx(4,4)*theta34+bx(4,5)*theta35+...
l*(by(3,1)*theta14+by(3,2)*theta24+by(3,3)*theta34+by(3,4)*theta44);
%
% Algebraic equations
%
fx(7)=ax(1,1)*theta11+ax(1,2)*theta12+ax(1,3)*theta13+ax(1,4)*theta14+ax(1,5)*theta15;
fx(8)=ay(1,1)*theta12+ay(1,2)*theta22+ay(1,3)*theta32+ay(1,4)*theta42;
fx(9)=ay(1,1)*theta13+ay(1,2)*theta23+ay(1,3)*theta33+ay(1,4)*theta43;
fx(10)=ay(1,1)*theta14+ay(1,2)*theta24+ay(1,3)*theta34+ay(1,4)*theta44;
fx(11)=ax(5,1)*theta11+ax(5,2)*theta12+ax(5,3)*theta13+ax(5,4)*theta14+ax(5,5)*theta15;
fx(12)=ax(1,1)*theta41+ax(1,2)*theta42+ax(1,3)*theta43+ax(1,4)*theta44+ax(1,5)*theta45;
fx(13)=ay(4,1)*theta12+ay(4,2)*theta22+ay(4,3)*theta32+ay(4,4)*theta42+beta*theta42;
fx(14)=ay(4,1)*theta13+ay(4,2)*theta23+ay(4,3)*theta33+ay(4,4)*theta43+beta*theta43;
fx(15)=ay(4,1)*theta14+ay(4,2)*theta24+ay(4,3)*theta34+ay(4,4)*theta44+beta*theta44;
fx(16)=ax(5,1)*theta41+ax(5,2)*theta42+ax(5,3)*theta43+ax(5,4)*theta44+(ax(5,5)+alpha)*theta45;
fx(17)=ax(1,1)*theta21+ax(1,2)*theta22+ax(1,3)*theta23+ax(1,4)*theta24+ax(1,5)*theta25;
fx(18)=ax(1,1)*theta31+ax(1,2)*theta32+ax(1,3)*theta33+ax(1,4)*theta34+ax(1,5)*theta35;
fx(19)=ax(5,1)*theta21+ax(5,2)*theta22+ax(5,3)*theta23+ax(5,4)*theta24+(ax(5,5)+alpha)*theta25;
fx(20)=ax(5,1)*theta31+ax(5,2)*theta32+ax(5,3)*theta33+ax(5,4)*theta34+(ax(5,5)+alpha)*theta35;

fx = fx’;

%-- End of the ex5ode.m file --

36CHAPTER 1. ECUACIONES DIFERENCIALES ORDINARIAS: PROBLEMAS DE VALORES EN

Ejemplo 8 3 Una cavidad cilı́ndrica, de 10 pies de diametro y 10 pies de altura, está

enterrada en el suelo que se mantiene inicialmente a una temperatura uniforme de 60 o F.
Al tiempo t = 0, la cavidad se llena con gas natural licuado, manteniendo sus paredes a -260
o
F. El gas licuado que se evapora se recondensa y se regresa a la cavidad. La temperatura
en el suelo que rodea a la cavidad se puede determinar resolviendo la siguiente ecuación,
µ 2 ¶
∂ T 1 ∂T ∂ 2T ∂T
k 2
+ + 2 = ρCp
∂r r ∂r ∂z ∂t

donde r y z son la coordenada radial y vertical, respectivamente. El suelo está seco y tiene
las siguientes propiedades: conductividad térmica k=2.2 BTU/(hr-pie-o F), densidad ρ=125
lb/pie3 y capacidad especifica Cp =0.19 BTU/(lb-o F). Graficar los perfiles de temperatura
sobre ambos ejes usando un periodo largo de tiempo, asi como el flujo de calor Q desde el
suelo hasta la cavidad.

1.5 Colocación Ortogonal sobre elementos finitos

1.5.1 Polinomios de Lagrange
Cuando el perfil de las variables dependientes varia fuertemente, el método de colo-
cación ortogonal puede producir resultados que no reflejan el perfil verdadero. Esto se
debe a que se ha supuesto en el método de colocación ortogonal que las funciones base
(polinomios ortogonales) son continuas sobre todo el espacio de solución.
Una manera de resolver este problema consiste en dividir el espacio de solución en
regiones o elementos no necesariamente de la misma longitud. De hecho los elementos
tienden a ser mas pequen̄os en las regiones donde ocurren cambios pronunciados en
los perfiles de las variables dependientes. Una vez que los perfiles son mas “suaves”, o
muestran mayor continuidad, el taman̄o de los elementos puede incrementarse. Obser-
verse que la forma de discretizar el espacio de solución es muy parecida a la forma como
se elige el taman̄o de paso en problemas de valores iniciales de ecuaciones diferenciales
ordinarias: el taman̄o de paso se hace mas pequen̄o alrededor de la región donde la
conducta dinámica cambia rapidamente; el taman̄o del paso se hace mas grande cuando
la dinámica no cambia demasiado.
Una vez realizada la discretización, o elección de la localización de los elementos,
el método de colocación ortogonal se aplica sobre cada elemento en la misma forma
realizada antes. Cuando se aplica el método de colocación en esta forma, existe el
riezgo de que la solución obtenida en el interior del elemento pueda ser discontinua
cuando se realiza el cambio hacia el siguiente elemento. Normalmente en casos como
este se acostumbra a imponer una restricción extra para asegurar la continuidad de
la solución entre los elementos. Esto es particularmente cierto si las funciones base
usadas son polinomios discontinuos. 4
3
Tomado de Carnahan, Luther, Wilkes, Applied Numerical Methods, John Wile, 1969, Problema
7.25.
4
Existe una forma simple de garantizar que la solución sea continua entre los elementos. Esto
puede lograrse utilizando polinomios de Hermite como funciones base.
Colocación Ortogonal sobre elementos finitos 37

primer segundo
elemento elemento
u u u u

x1 x2 x3 x4 x5 x6 x7
0 .5 1

Figure 1.11: El espacio de solución se divide en dos elementos finitos equidistantes y

cada uno de ellos tiene dos puntos internos de colocación.

Ejemplo 9 Supongamos que deseamos resolver el siguiente problema,

d2 C
= φ2 C (1.5.61)
dx2
sujeto a las condiciones frontera,
dC
= 0 en x = 0 (1.5.62)
dx
C = 1 en x = 1 (1.5.63)
el espacio de solución se divide (arbitrariamente) en dos elementos equidistantes (es decir la longitud
de cada elemento finito es de h1 = h2 = .5), y se emplean dos puntos internos de colocación en
el interior de cada elemento tal como se muestra en la figura 1.11. En esta misma figura x2 , x3
representan los puntos de colocación internos del primer elemento, x5 , x6 son los puntos de colocación
del segundo elemento. x1 , x4 son los puntos de colocación en los extremos o fronteras del primer
elemento, mientras que x4 , x7 son los respectivos puntos externos de colocación del segundo elemento.

Si tuvieramos solo un elemento los puntos de colocación internos estarı́an úbicados en .21132 y
.78868. Sin embargo, si usamos dos elementos, los puntos de colocación referidos a estos dos elementos
estarı́an localizados en:

x1 =0
x2 = 12 (.21132) = .10566
x3 = 12 (.78868) = .39434
x4 = .5
x5 = .5 + 12 (.21132) = .60566
x6 = .5 + 12 (.78868) = .89434
x7 = 1.

El método de colocación ortogonal sobre los dos elementos finitos puede ahora aplicarse en la
siguiente secuencia.5
1. Primer elemento.

5
Notese que uno puede plantear las ecuaciones en cualquier orden, pero el orden sugerido puede
conducir a menores errores en el planteamiento del método.
38CHAPTER 1. ECUACIONES DIFERENCIALES ORDINARIAS: PROBLEMAS DE VALORES EN

a).- Condición frontera en x = 0.

dC
=0
dx
¯ X4
dC ¯¯
= A1i Ci = A11 C1 + A12 C2 + A12 C3 + A14 C4 = 0 (1.5.64)
dx ¯x=x1 i=1

b).- Modelo en los puntos internos de colocación.

d2 C
= φ2 C
dx2
antes de discretizar esta ecuación notese claramente que la variable independiente x toma
valores entre 0 y 1. Sin embargo, la longitud de este elemento fue definida como h1 = .5,
necesitamos escalar a la variable x para referirla a esta longitud. Esto puede hacerse
facilmente definiendo,
x − x1 x − x1
u= =
x4 − x1 h1

entonces el modelo se reescribirı́a como,

1 d2 C
= φ2 C
h1 du2

• primer punto (x = x2 ).
4
1 d2 C 1 X 1
2
= B2i Ci = (B21 C1 + B22 C2 + B23 C3 + B24 C4 )
h1 du h1 i=1 h1

φ2 C = φ2 C2

entonces,

1 d2 C 1
2
− φ2 C = (B21 C1 + B22 C2 + B23 C3 + B24 C4 ) − φ2 C2 = 0 (1.5.65)
h1 du h1

• segundo punto (x = x3 ).

1 d2 C
= φ2 C
h1 du2

4
1 d2 C 1 X 1
2
= B3i Ci = (B31 C1 + B32 C2 + B33 C3 + B34 C4 )
h1 du h1 i=1 h1

φ2 C = φ2 C3

entonces,

1 d2 C 1
− φ2 C = (B31 C1 + B32 C2 + B33 C3 + B34 C4 ) − φ2 C3 = 0 (1.5.66)
h1 du2 h1
Colocación Ortogonal sobre elementos finitos 39

c).- Condición para asegurar continuidad entre los elementos 1 y 2.

Para asegurar continuidad de la solución entre los elementos basta con plantear la conser-
vación de algún término en el primer elemento. Si esta continuidad se conserva, digamos
para los puntos de colocación 1,2,3 y 4, entonces resulta que se asegura continuidad de
la solución ya que el punto de colocación 4 es también el primer punto de colocación del
segundo elemento.
En este problema, ya que el sistema opera en estado estacionario, se puede postular la
conservación del flujo másico,

N1 = N2 = N3 = N4

siendo,
dCi
Ni = Di
du
donde i representa puntos de colocación. Entonces planteando flujo constante de masa
entre el primer y último punto de colocación del primer elemento,

N1 = N4

o sea,
à ¯ ! à ¯ !
D4 dC ¯¯ D1 dC ¯¯
=
h1 du ¯x=x4 h2 du ¯x=x1
1 2

donde los subindices mas externos 1 y 2 se refieren al primer y segundo elementos respec-
tivamente. Si suponemos que las difusividades Di son constantes,
¯ 4
dC ¯¯ X
= A4i Ci = A41 C1 + A42 C2 + A43 C3 + A44 C4
du ¯x=x4 i=1

¯ 4
dC ¯¯ X
= A1i Ci+3 = A11 C4 + A12 C5 + A13 C6 + A14 C7
du ¯x=x1 i=1

finalmente,

N1 −N4 = A41 C1 +A42 C2 +A43 C3 +(A44 −A11 )C4 +A12 C5 +A13 C6 +A17 C7 = 0 (1.5.67)

2. Segundo elemento.

a).- Modelo en los puntos de colocación internos.

• primer punto (x = x5 ).
4
1 d2 C 1 X 1
= B2i Ci+3 = (B21 C4 + B22 C5 + B23 C6 + B24 C7 )
h2 du2 h2 i=1 h2

φ2 C = φ2 C5

entonces,

1 d2 C 1
− φ2 C = (B21 C4 + B22 C5 + B23 C6 + B24 C7 ) − φ2 C5 = 0 (1.5.68)
h2 du2 h2
40CHAPTER 1. ECUACIONES DIFERENCIALES ORDINARIAS: PROBLEMAS DE VALORES EN

• segundo punto (x = x6 ).

1 d2 C
= φ2 C
h2 du2

4
1 d2 C 1 X 1
= B3i Ci+3 = (B31 C4 + B32 C5 + B33 C6 + B34 C7 )
h2 du2 h2 i=1 h2

φ2 C = φ2 C6

entonces,

1 d2 C 1
2
− φ2 C = (B31 C4 + B32 C5 + B33 C6 + B34 C7 ) − φ2 C6 = 0 (1.5.69)
h2 du h2

b).- Condición frontera en x = 1.

C=1

o bien,
C7 = 1 (1.5.70)

las ecuaciones 1.5.64-1.5.70 pueden reescribirse en forma matricial como se muestra a continuación.
    
A11 A12 A13 A14 C1 0
 B21 B22 B23 B24   C2   h1 φ2 C2 
    
 B31 B32 B33 B34   C3   h1 φ2 C3 
    
 A41 A42 A43 (A44 − A11 ) −A12 −A13 −A14   C4  =  0 
    
 B21 B22 B23    
B24   C5   h2 φ C5  2
 
 B31 B32 B33 B34   C6   h2 φ2 C6 
1 C7 1

haciendo φ = 6 y utilizando las siguientes matrices de colocación,

 
-7 8.196 -2.196 1
 -2.732 1.732 1.732 -.7321 
A=  .7321 -1.7321 -1.7321 2.732 


-1 2.96 -8.196 7

 
24 -37.18 25.18 -12
 16.39 -24 12 -4.392 
B=
 -4.392

12 -24 16.39 
-12 25.18 -37.18 24

obtenemos las soluciones en los puntos de colocación,

C1 = 0.00597
C2 = 0.00608
C3 = 0.02863
C4 = 0.05478
C5 = 0.08308
C6 = 0.51965
C7 =1
Colocación Ortogonal sobre elementos finitos 41

1
2
3
0.9
4
5
6
0.8

0.7

0.6
C

0.5

0.4

0.3

0.2

0.1

0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
x

Figure 1.12: Perfil de concentración usando diferente número de puntos internos de

colocación y dos elementos finitos.

En la Figura 1.12 se muestran la forma como la solución se mejora a medida que se aumenta el
número de puntos de colocación internos empleando dos elementos finitos. Por supuesto, la calidad
de la solución puede mejorarse de manera más fácil si usamos varios elementos finitos aun empleando
relativamente pocos puntos internos de colocación. En la figura 1.13 se muestra la solución del prob-
lema en cuestióon usando varios elementos finitos y siempre dos puntos internos de colocación. Como
puede notarse, el uso de elementos finitos coadyuva a obtener perfiles de concentración más cercanos
a la solución esperada. A continuación se meustra el programa Matlab empleado para realizar los
cálculos mostrados en la figura anterior.

clear all; clc;

global nfe npc a b roots h phi

npc = 2; % Number of internal collocation pints

nfe = 12; % Number of finite elements

phi = 6;

[a,b,q,roots,w] = planar(npc);

sum = 0;
for j = 1:nfe,
h(j) = 1/nfe; % finite element length
sum = sum+h(j);
hx(j) = sum;
42CHAPTER 1. ECUACIONES DIFERENCIALES ORDINARIAS: PROBLEMAS DE VALORES EN

1.2 1

0.9
1

0.8

0.8 0.7

0.6
0.6
C

C
0.5

0.4
0.4

0.3
0.2

0.2

0
0.1

−0.2 0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
x x

(a) (b)
1 1

0.9 0.9

0.8 0.8

0.7 0.7

0.6 0.6
C

0.5 0.5

0.4 0.4

0.3 0.3

0.2 0.2

0.1 0.1

0 0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
x x

(c) (d)
1 1

0.9 0.9

0.8 0.8

0.7 0.7

0.6 0.6
C

0.5 0.5

0.4 0.4

0.3 0.3

0.2 0.2

0.1 0.1

0 0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
x x

(e) (f)

Figure 1.13: Resultados del ejemplo 9 empleando: (a) 2, (b) 4, (c) 6, (d) 8, (e) 10 y
(f) 12 elementos finitos. En todos los casos siempre se usaron 2 puntos internos de
colocación.
Colocación Ortogonal sobre elementos finitos 43

end

index = 0;
for k = 1:nfe,
if k == 1
for j = 1:npc+2,
index = index+1;
u(j) = roots(j)*h(k);
end
else
for j = 1:npc+1,
index = index+1;
u(index) = u(index-(npc+1))+h(k);
end
end
end

x0 = ones(nfe*(npc+2)-(nfe-1),1);
options = optimset(’display’,’iter’);

x = fsolve(’difussion_fe’,x0,options)

plot(u,x,’o-’),xlabel(’x’),ylabel(’C’),hold
for j = 1:nfe-1,
plot([hx(j) hx(j)],[0 1],’r--’)
end

%-- End of the run.m file --

function fx = difussion_fe(c)

global nfe npc a b roots h phi

npc1 = npc+1;
npc2 = npc+2;
%
% Boundary condition at x=0
%
suma = 0;
for i = 1:npc2,
suma = suma+a(1,i)*c(i);
end
fx (1) = suma;
%
% Model approximation within the first element
%
index = 1;
for k = 1:nfe,
for j = 2:npc1,
index = index+1;
sumb = 0;
for i = 1:npc2,
sumb = sumb+b(j,i)*c(i+(k-1)*npc1);
end
44CHAPTER 1. ECUACIONES DIFERENCIALES ORDINARIAS: PROBLEMAS DE VALORES EN

fx (index) = sumb/(h(k)*h(k))-phi*phi*c(j+(k-1)*npc+1);
end
%
% Continuity condition between finite elements
%
if k ~= nfe,
sumlhs = 0 ;
sumrhs = 0 ;
for i = 1:npc2,
sumlhs = sumlhs+a(npc2,i)*c(i+(k-1)*npc1);
sumrhs = sumrhs+a(1,i)*c(i+(k-1)*npc1+npc1);
end
index = index+1;
fx (index) = sumlhs-sumrhs;
end
end
%
% Boundary condition at x = 1;
%
index = index+1;
fx (index) = c(index)-1;

%-- End of the difussion_fe.m file

Ejemplo 10 Resolver el siguiente problema de valores a la frontera empleando el método

de Colocación Ortogonal sobre elementos finitos,

2d2 y dy
(² + x ) 2 + 4x + 2y = 0
dx dx
sujeto a las siguientes condiciones frontera:
1
y(−1) = y(1) =
1+²
resolver el problema para valores de ² = 0.05, 0.02 y 0.01. Graficar la respuesta del estado
y como función de x.

1.5.2 Polinomios de Hermite

Existe una manera de garantizar continuidad de la solución entre los diferentes ele-
mentos de que puede constar un sistema. Esta forma consiste en emplear polinomios
cúbicos de Hermite, los cuales son, por su propia definición, continuos entre los ele-
mentos. Por ejemplo, para un elemento k determinado tenemos 4 polinomios cúbicos,
donde los polinomios Hi , i = 1, 4 se definen como (ver figura 1.14),

H1 = (1 − u)2 (1 + 2u)
H2 = u(1 − u)2 hk
H3 = u2 (3 − 2u)
H4 = u2 (u − 1)hk
Colocación Ortogonal sobre elementos finitos 45

1.2

0.8 H
1
H
3

0.6
i
H

0.4

0.2
H2

0
H
4

−0.2
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
u

Figure 1.14: Polinomios cúbicos de Hermite.

entonces la solución c(u) en cualquier punto del dominio se obtiene mediante la siguiente
combinación lineal,
X4
C(u) = ai Hi (u) (1.5.71)
i=1
donde la primera derivada de la solución está dada por,
4
dC 1 dC 1 X dHi
= = ai
dx hk du hk i=1 du

mientras la segunda derivada,

4
d2 C 1 X d2 Hi
= ai
dx2 hk i=1 du2

la aplicación del método de colocación ortogonal ocurre cuando evaluamos la solución

y sus primera y segunda derivadas en los puntos internos de colocación uj . Entonces,
4
X
C(uj ) = Hi (uj )ai
i=1

¯ 4 ¯
dC ¯¯ X dHi ¯¯
= ai
du ¯uj i=1
du ¯uj
46CHAPTER 1. ECUACIONES DIFERENCIALES ORDINARIAS: PROBLEMAS DE VALORES EN
¯ X4 ¯
d2 C ¯¯ d2 Hi ¯¯
= ai
du2 ¯uj i=1
du 2 ¯
uj

o bien,
4
X
C(uj ) = Hji ai
i=1

X 4
dC
(uj ) = Aji ai
du i=1

X 4
d2 C
(u j ) = Bji ai
du2 i=1

donde, por ejemplo, si usamos 2 puntos internos de colocación ( u1 = 0, u2 = .2113248654, u3 =

.7886751346, u4 = 0) las matrices H, A y B se evaluan de la siguiente manera.
• Matriz H.

 
H1 (u1 ) H2 (u1 ) H3 (u1 ) H4 (u1 )
 H1 (u2 ) H2 (u2 ) H3 (u2 ) H4 (u2 ) 
H=
 H1 (u3 )

H2 (u3 ) H3 (u3 ) H4 (u3 ) 
H1 (u4 ) H2 (u4 ) H3 (u4 ) H4 (u4 )
 
1 0 0 0
 .8849018 .131445856hk .11509982 -.035220811hk 
H=
 .11509982 .035220811hk

.88490018 -.131445856hk 
0 0 1 0

• Matriz A.
 dH1 (u1 ) dH2 (u1 ) dH3 (u1 ) dH4 (u1 )

du du du du
 dH1 (u2 ) dH2 (u2 ) dH3 (u2 ) dH4 (u2 ) 
 
A= du
dH1 (u3 )
du
dH2 (u3 )
du
dH3 (u3 )
du
dH4 (u3 ) 
 du du du du

dH1 (u4 ) dH2 (u4 ) dH3 (u4 ) dH4 (u4 )
du du du du

donde,
dH1 (uj )
= −2(1 − uj )(1 + 2uj ) + 2(1 − uj )2
du
dH2 (uj )
= (1 − uj )2 hk − 2uj (1 − uj )hk
du
dH3 (uj )
= 2uj (3 − 2uj ) − 2u2j
du
dH4 (uj )
= 2uj (uj − 1)hk + u2j hk
du
Colocación Ortogonal sobre elementos finitos 47

entonces,
 
0 hk 0 0
 -1 0.288675316hk 1 -0.288675316hk 
A=
 -1 -0.288675316hk

1 0.288675316hk 
0 0 0 hk

• Matriz B.
 d2 H1 (u1 ) d2 H2 (u1 ) d2 H3 (u1 ) d2 H4 (u1 )

du2 du2 du2 du2
 d2 H1 (u2 ) d2 H2 (u2 ) d2 H3 (u2 ) d2 H4 (u2 ) 
 du2 du2 du2 du2 
B= d2 H1 (u3 ) d2 H2 (u3 ) d2 H3 (u3 ) d2 H4 (u3 ) 
 du2 du2 du2 du2

d2 H1 (u4 ) d2 H2 (u4 ) d2 H3 (u4 ) d2 H4 (u4 )
du2 du2 du2 du2

donde,
d2 H1 (uj )
= −6 + 12uj
du2
d2 H2 (uj )
= −4(1 − uj )hk + 2uj hk
du2
2
d H3 (uj )
= 6 − 12uj
du2
2
d H4 (uj )
= 2(uj − 1)hk + 4uj hk
du2
entonces,
 
-6 -4hk 6 -2hk
 -3.46410162 -2.73205081hk 3.46410162 -0.73206081hk 
B=
 3.46410162 0.73205081hk

-3.4641062 2.73205081hk 
6 2hk -6 4hk

Para demostrar que el empleo de polinomios de Hermite, como funciones base,

permite garantizar que la solución es continua entre los elementos se debe cumplir,
µ ¯ ¶ µ ¯ ¶
dC ¯¯ dC ¯¯
= (1.5.72)
dx ¯x=1 k dx ¯x=0 k+1

donde el subindice k se refiere al número del elemento. Entonces,

• en x = 1 (elemento k)
µ ¶
dC 1 dh1 dh2 dh3 dh4
= a1 + a2 + a3 + a4
dx hk du du du du

• en x = 0 (elemento k + 1)
µ ¶
dC 1 dh1 dh2 dh3 dh4
= a3 + a4 + a5 + a6
dx hk du du du du
48CHAPTER 1. ECUACIONES DIFERENCIALES ORDINARIAS: PROBLEMAS DE VALORES EN

elemento elemento
k k+1
u u u u

1, 2 3, 4 5, 6

Figure 1.15: El espacio de solución se divide en dos elementos finitos equidistantes y

cada uno de ellos tiene dos puntos internos de colocación.

observese que los coeficientes a3 , a4 que estan localizados en el extremo derecho del
elemento k son los mismos que los elementos a3 , a4 localizados en el extremo izquierdo
del elemento k + 1 (ver figura 1.15). Entonces,
dh1
= 2(1 − u)2 − 2(1 + 2u)(1 − u)
du
dh2
= −2u(1 − u)hk + (1 − u)2 hk
du
dh3
= −2u2 + 2u(3 − 2u)
du
dh4
= u2 hk + 2u(u − 1)hk
du
• en u = 1
dh1 dh2 dh3 dh4
= 0, = 0, = 0, = hk
du du dt dt
• en u = 0
dh1 dh2 dh3 dh4
= 0, = hk+1 , = 0, =0
du du du du
por lo tanto,
1 1
a4 hk = a4 hk+1
hk hk+1
de esta ecuación puede deducirse que la primera derivada de la solución es continua
entre los elementos.
Ejemplo 11 Resuelva el mismo problema de valores a la frontera planteado en el ejemplo 9 usando
polinomios de Hermite como funciones base. Utilize nuevamente dos elementos de igual taman̄o
h1 = h2 = .5 y dos puntos de colocación internos en cada elemento.
Como en el caso del ejemplo 9 planteamos las ecuaciones comenzando desde el extremo izquiedo
del dominio de solución y progresando asi hasta alcanzar el extremo derecho del espacio de solución.
Colocación Ortogonal sobre elementos finitos 49

1. Primer elemento.

a).- Condición frontera (x=0)

dC
=0
dx
o sea,
A11 a1 + A12 a2 + A13 a3 + A14 a4 = 0 (1.5.73)

b).- Modelo matemático.

d2 C
= φ2 C
dx2
el cual en términos de la coordenada u del elemento se escribe como,

d2 C 1 d2 C
= = φ2 C
dx2 h1 du2

• primer punto de colocación.

1 d2 C 1 ¡ ¢
−φ2 C = 2 B21 a1 + B22 a2 + B23 a3 + B24 a4 ) − φ2 (H21 a1 + H22 a2 + H23 a3 + H24 a4
h21 du2 h1
(1.5.74)
• segundo punto de colocación.

1 d2 C 1 ¡ ¢
2 2
−φ2 C = 2 B31 a1 + B32 a2 + B33 a3 + B34 a4 ) − φ2 (H31 a1 + H32 a2 + H33 a3 + H34 a4
h1 du h1
(1.5.75)

2. Segundo elemento.

a).- Modelo matemático.

d2 C
= φ2 C
dx2
• primer punto de colocación.

1 d2 C 1 ¡ ¢
2 2
−φ2 C = 2 B21 a3 + B22 a4 + B23 a5 + B24 a6 ) − φ2 (H21 a3 + H22 a4 + H23 a5 + H24 a6
h2 du h2
(1.5.76)
• segundo punto de colocación.

1 d2 C 1 ¡ ¢
3 2
−φ2 C = 2 B31 a3 + B32 a4 + B33 a5 + B34 a6 ) − φ2 (H31 a3 + H32 a4 + H33 a5 + H34 a6
h2 du h2
(1.5.77)
b).- Condición frontera (x = 1).

C(1) = 1

o sea,
H41 a3 + H42 a4 + H43 a5 + H44 a6 = 1 (1.5.78)
haciendo nuevamente φ = 6 y resolviendo el sistema de ecuaciones 1.5.73-1.5.78 obten-
emos los coeficientes ai de la expansión representada por la ecuación 1.5.71,
50CHAPTER 1. ECUACIONES DIFERENCIALES ORDINARIAS: PROBLEMAS DE VALORES EN

a1 = 0.006
a2 =0
a3 = 0.0548
a4 = 0.3125
a5 =1
a6 = 5.737

notese que una vez que se han determinado los coeficientes ai , se puede determinar la
solución a lo largo del todo el espacio de solución x empleando la ecuación 1.5.71.

c1 = H11 a1 + H12 a2 + H13 a3 + H14 a4 = 0.0060

c2 = H21 a1 + H22 a2 + H23 a3 + H24 a4 = 0.0061
c3 = H31 a1 + H32 a2 + H33 a3 + H34 a4 = 0.0286
c4 = H41 a1 + H42 a2 + H43 a3 + H44 a4 = 0.0548
c5 = H21 a3 + H22 a4 + H23 a5 + H24 a6 = 0.0831
c6 = H31 a3 + H32 a4 + H33 a5 + H34 a6 = 0.5197
c7 = H41 a3 + H42 a4 + H43 a5 + H44 a6 =1

En la Figura 1.16 se muestra el perfil de concentracción contra la posición en cada uno

de los 2 elementos finitos.

A continuación se muestra el programa Matlab empleado para realizar los cálculos.

0.9

0.8

0.7

0.6
C

0.5

0.4

0.3

0.2

0.1

0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
x

Figure 1.16: Perfil de concentración usando dos puntos internos de colocación y dos
elementos finitos.
Colocación Ortogonal sobre elementos finitos 51

clear all; clc;

global h a b h1 h2 phi

n = 2; % internal collocation points

h1 = 0.5; % Length of first finite element
h2 = 0.5; % Length of second finite element
phi = 6; % Thiele module

[h,a,b,roots]=hermite(n,h1);

x0 = ones(6,1); options = optimset(’display’,’iter’);

a = fsolve(’example11’,x0,options);

% Concentration at collocation points

index = 0; for j = 1:n+2,

index = index+1;
sum = 0;
for i = 1:4,
sum = sum+h(j,i)*a(i);
end
c (index) = sum;
x (index) = roots(index)*h1;
end

for j = 2:n+2,
index = index+1;
sum = 0;
for i = 1:4,
sum = sum+h(j,i)*a(i+2);
end
c (index) = sum;
x (index) = h1+roots(j)*h2;
end

plot(x,c,’o-’),hold,xlabel(’x’),ylabel(’C’)
plot([h1 h1],[0 1],’--’)

%-- end of the run.m file --

function fx = example11(x)

global h a b h1 h2 phi

fx(1) = a(1,1)*x(1) + a(1,2)*x(2) + a(1,3)*x(3) + a(1,4)*x(4);

fx(2) = (1/(h1*h1))*(b(2,1)*x(1)+b(2,2)*x(2)+b(2,3)*x(3)+b(2,4)*x(4))...
-phi*phi*(h(2,1)*x(1) + h(2,2)*x(2) + h(2,3)*x(3) + h(2,4)*x(4));

fx(3) = (1/(h1*h1))*(b(3,1)*x(1) + b(3,2)*x(2) + b(3,3)*x(3)+b(3,4)*x(4))...

-phi*phi*(h(3,1)*x(1) + h(3,2)*x(2) + h(3,3)*x(3) + h(3,4)*x(4));
52CHAPTER 1. ECUACIONES DIFERENCIALES ORDINARIAS: PROBLEMAS DE VALORES EN

fx(4) = (1/(h2*h2))*(b(2,1)*x(3) + b(2,2)*x(4) + b(2,3)*x(5) +b(2,4)*x(6))...

-phi*phi*(h(2,1)*x(3) + h(2,2)*x(4) + h(2,3)*x(5) + h(2,4)*x(6));

fx(5) = (1/(h2*h2))*(b(3,1)*x(3) + b(3,2)*x(4) + b(3,3)*x(5) +b(3,4)*x(6))...

-phi*phi*(h(3,1)*x(3) + h(3,2)*x(4) + h(3,3)*x(5) + h(3,4)*x(6));

fx(6) = h(4,1)*x(3) + h(4,2)*x(4) + h(4,3)*x(5) + h(4,4)*x(6) - 1;

%---end of the example11.m file---

function [h,a,b,roots] = hermite (n,hk)

%
% M file to compute a fourth degree Hermite Polynomial
%
% Inputs:
% -------
% n = number of internal collocation points (without including boundary points)
% hk = length of the finite element
%
% Outputs:
% --------
% h[n,4] = Hermite polynomial values
% a[n,4] = First order derivative of Hermite Polynomial
% b[n,4] = Second order derivative of Hermite Polynomial
% roots[n+2] = Roots of an orthogonal polynomial
%
% Written by Antonio Flores T./ 19,Nov, 2005
%
[aoc,boc,qinv,roots,woc] = planar(n);

for j = 1:n+2,
u = roots(j);
%
% Polynomial values
%
h1 = (1-u)^2*(1 + 2*u);
h2 = u*(1 -u)^2*hk;
h3 = u^2*(3-2*u);
h4 = u^2*(u-1)*hk;
%
% First order derivatives
%
h1p =-2*(1-u)*(1+2*u)+2*(1-u)^2;
h2p = (1-u)^2*hk-2*u*(1-u)*hk;
h3p = 2*u*(3-2*u)-2*u^2;
h4p = 2*u*(u-1)*hk+u^2*hk;
%
% Second order derivatives
%
h1pp = -6+12*u;
h2pp = -4*(1-u)*hk+2*u*hk;
h3pp = 6-12*u;
Colocación Ortogonal sobre elementos finitos 53

h4pp = 2*(u-1)*hk+4*u*hk;
%
% Fill the matrices
%
h(j,1) = h1 ; h(j,2) = h2 ; h(j,3) = h3 ; h(j,4) = h4;
a(j,1) = h1p ; a(j,2) = h2p ; a(j,3) = h3p ; a(j,4) = h4p;
b(j,1) = h1pp; b(j,2) = h2pp; b(j,3) = h3pp; b(j,4) = h4pp;
end

%-- End of the hermite.m file

54CHAPTER 1. ECUACIONES DIFERENCIALES ORDINARIAS: PROBLEMAS DE VALORES EN

1.6 Diferencias Finitas

Uno de los metodos frecuentemente usados para la solucion numerica de ecuaciones
diferenciales (ya sean problemas de valores iniciales o problemas de valores en la fron-
tera) se basa en la aproximacion de la solucion continua por una solucion discreta
determinada en un conjunto de puntos dados del espacio de solucion. Dicha aproxi-
macion reduce la dimensionalidad del sistema, ya que se consigue transformar asi un
sistema de dimension infinita en un sistema de dimension finita. Ademas la aproxi-
macion tambien permite transformar la ecuacion diferencial original en una ecuacion
algebraica (lineal si el operador diferencial es lineal, y no-lineal si el operador diferen-
cial es no-lineal) cuya solucion es mas simple que el problema inicial. Bajo el titulo
diferencias finitas se pueden englobar muchos metodos de aproximacion (incluyendo
aquellos metodos empleados para la solucion de problemas de valores iniciales); sin em-
bargo, bajo este nombre, en este parte unicamente nos referiremos al tipo de metodos
de aproximacion que substituyen las derivadas, de cualquier orden, por su equivalente
discreto en una serie de puntos seleccionados a priori. Esta serie de puntos se conocen
como nodos. Esto significa que la solucion del problema unicamente sera conocida en
los nodos seleccionados; sin embargo, es posible determinar la solucion en cualquier
punto de interes del espacio de solucion usando algun esquema de interpolacion. En
principio podemos elegir arbitrariamente nodos igualmente espaciados; sin embargo,
en algunas aplicaciones, el uso de espaciamiento variable podria ser importante como
una forma de disminuir el error de aproximacion. En terminos globales la aplicacion
del metodo de diferencias finitas involucra los siguientes pasos.

1. Discretizacion del espacio de solucion. El espacio de solucion se divide en

una serie de nodos.

2. Discretizacion de las derivadas. En cada uno de los nodos seleccionados,

las derivadas continuas del sistema original se substituyen por sus equivalentes
discretos.

3. Solucion del sistema resultante. Las ecuaciones resultantes del paso anterior
se resuelven numericamente empleando el metodo apropiado.

En esta seccion aplicaremos el metodo de diferencias finitas para la aproximacion nu-

merica de ecuaciones diferenciales ordinarias, tanto para problemas de valores iniciales
como de problemas de valores en la frontera.

1.6.1 Aproximacion de derivadas en una dimension.

La aproximacion discreta de la derivada de una funcion continua y(x) puede obtenerse
empleando la definicion clasica de la derivada

dy y(x + h) − y(x) y(x + h) − y(x)

≡ lim = lim (1.6.79)
dx (x + h) − x h
h→0 h→0
Diferencias Finitas 55

dy
usando un valor finito de h obtenemos el valor aproximado de la derivada dx en un
nodo determinado xi . La aproximacion numerica de la derivada estara dada entonces
por ¯
dy ¯¯ y(xi + h) − y(xi ) yi+1 − yi
¯ ≈ = (1.6.80)
dx xi (xi + h) − xi xi+1 − xi
donde yi ≡ y(xi ), yi+1 ≡ y(xi + h), xi+1 ≡ xi + h, xi es un nodo. La aproximacion
obtenida usando la ecuacion 1.6.80 permite transformar la derivada por un termino
puramente algebraico. Por supuesto la aproximacion anterior no es la unica que se
puede emplear para discretizar una derivada. Otras opciones son las siguientes
¯
dy ¯¯ yi − yi−1
¯ ≈ (1.6.81)
dx xi xi − xi−1
¯
dy ¯¯ yi+1 − yi−1
¯ ≈ (1.6.82)
dx xi xi+1 − xi−1
donde xi−1 es el nodo anterior al nodo xi .
Las derivadas de ordenes superiores se obtienen de manera anologa al caso de primer
orden. Por ejemplo para una derivada de segundo orden
yi+1 −yi yi −yi−1
d2 y xi+1 −xi
− xi −xi−1
≈ (1.6.83)
dx2 xi+1 − xi−1
si suponemos que el espaciamiento entre los nodos es constante: xi+1 − xi = xi − xi−1 =
∆x, la ecuacion anterior puede reescribirse como
d2 y yi+1 − 2yi + yi−1
≈ (1.6.84)
dx2 (∆x)2

1.6.2 Analisis del error de truncamiento y consistencia

Dado que el metodo de diferencias finitas constituye tan solo una aproximacion a la
solucion verdadera, resulta importante obtener una idea del error que se introduce
cuando las derivadas de la funcion original se aproximan usando el enfoque de diferen-
cias finitas. Tal error se conoce como error de truncamiento dado que, como veremos
mas adelante, la aproximacion de las derivadas se ha realizado usando solo el primer
termino de una expansion que contiene un numero infinito de terminos. La suma de
estos terminos despreciados constituye por lo tanto el error de truncamiento.
Para obtener una expresion analitica del error de truncamiento empleamos el teo-
rema de Taylor de aproximacion de una funcion continua en un punto.
Teorema 1 Considere la funcion y(x) que es continua en el intervalo [ω1 , ω2 ] y que
tiene derivadas continuas hasta de orden N tambien en el mismo intervalo [ω1 , ω2 ].
Dando un punto x∗ ∈ [ω1 , ω2 ] entonces para cualquier punto x ∈ [ω1 , ω2 ],
¯ ¯ ¯
dy ¯¯ d2 y ¯¯ (x − x∗ )2 dN −1 y ¯¯ (x − x∗ )N −1
y(x) = y(x∗ ) + (x − x∗ ) + 2 ¯ + ... + + RN
dx ¯x∗ dx x∗ 2! dxN −1 ¯x∗ (N − 1)!
(1.6.85)
56CHAPTER 1. ECUACIONES DIFERENCIALES ORDINARIAS: PROBLEMAS DE VALORES EN

donde
¯
N dN y ¯¯ (x − x∗ )N
R ≡ , ξ ∈ [ω1 , ω2 ]
dxN ¯ξ N!

notese que si la funcion y(x) es infinitamente diferenciable en el intervalo [ω1 , ω2 ] en-

tonces la serie dada por la ecuacion 1.6.85 se escribe como
X ∞ ¯
dn u ¯¯ (x − x∗ )n
y(x) = (1.6.86)
n=0
dxn ¯x∗ n!

El teorema de Taylor se aplica ahora para obtener expansiones, alrededor del nodo
xi , de la solucion y(x) en los nodos xi+1 y xi−1
¯ ¯
dy ¯¯ d2 y ¯¯ (xi+1 − xi )2
y(xi+1 ) ≡ yi+1 = y(xi ) + (xi+1 − xi ) + 2 ¯ + ... + RN (1.6.87)
dx ¯xi dx xi 2!
¯ ¯
dy ¯¯ d2 y ¯¯ (xi−1 − xi )2
y(xi−1 ) ≡ yi−1 = y(xi ) + (xi−1 − xi ) + 2 ¯ + ... + RN (1.6.88)
dx ¯xi dx xi 2!
dy
despejando el termino dx
de la ecuacion 1.6.87
¯ ¯
dy ¯¯ yi+1 − yi xi+1 − xi d2 y ¯¯
= − + ... (1.6.89)
dx ¯xi xi+1 − xi 2 dx2 ¯xi

el lado derecho de la ecuacion anterior proporciona una medida del error de trun-
camiento. Comparando la ecuacion 1.6.80 con el lado derecho de la ecuacion 1.6.89
resulta que el error de truncamiento (ET) es el error en el que se incurre cuando se
desprecian todos los terminos restantes, a partir del segundo termino, del lado derecho
de la ecuacion 1.6.89
¯ ¯ ¯
dy ¯¯ yi+1 − yi xi+1 − xi d2 y ¯¯ (xi+1 − xi )2 d3 y ¯¯
ET ≡ − =− − + ... (1.6.90)
dx ¯xi xi+1 − xi 2 dx2 ¯xi 3! dx3 ¯xi

de la misma forma pueden obtenerse expresiones semejantes para el error de trun-

camiento para las otras formas de discretizacion de las derivadas dadas por las ecua-
ciones 1.6.81 y 1.6.82. Analizando la ecuacion 1.6.90 resulta claro que el error de
truncamiento puede minimizarse aumentando el numero de terminos en la expansion
de Taylor y/o reduciendo el espaciamiento entre los nodos.
Ejemplo 12 Difusión con reacción de primer orden-I.
La ecuación que describe la difusión en estado estacionario de una sustancia disuelta en el interior
de un fluido en reposo en el cual ocurre una reacción de primer orden, está dada por:
d2 c
D − Kc = 0 (1.6.91)
dx2
sujeto a las condiciones frontera

c(0) = 0
c(1) = c1
Diferencias Finitas 57

donde
c = concentraccion de la sustancia disuelta g/cm3
D = coeficiente de difusion = .01 cm2 /s
K = constante de velocidad de reaccion = .1 s−1
c1 = concentraccion en la frontera donde x = 1,=1 g/cm3

empleando el metodo de diferencias finitas determine el perfil de concentracción c(x) usando las sigu-
ientes discretizaciones sugeridas del espacio de solución.
i) • • •
1
x1 = 0 x2 = 2 x3 = 1

ii) • • • •
1 2
x1 = 0 x2 = 3 x3 = 3 x4 = 1

Solución.

• ∆x = .5
Como se comentó al inicio de esta sección, la aplicación del metodo de diferencias finitas para
la solución numerica de ecuaciones diferenciales involucra los siguientes pasos.

1. Discretización. Usaremos los nodos dados en el enunciado del problema; notese que la
localizacion de los nodos de los extremos coincide con las fronteras del problema.
2. Aproximación. Dado que desconocemos el valor de la solución en 3 nodos, debemos
plantear un sistema de 3 ecuaciones para determinar los valores de c(x1 ), c(x2 ) y c(x3 ).
Sin embargo, dos de los nodos están ubicados sobre las fronteras del sistema, y ya que
las condiciones frontera son de tipo Direchlet, esto fija automaticamente el valor de la
solución en c(x1 ) y c(x3 )
c(x1 ) = 0 (1.6.92)
c(x3 ) = 1 (1.6.93)
esta también es una forma de forzar a que el metodo de diferencias finitas satisfaga las
condiciones frontera. De las ecuaciones 1.6.92 y 1.6.93 la única incognita es c(x2 ) que
es el nodo donde aproximamos la derivada del problema, ecuación 1.6.91, en terminos de
diferencias finitas
¯ · ¸
d2 c ¯ c1 − 2c2 + c3
D 2 ¯¯ − Kc|x2 = D − Kc2 = 0 (1.6.94)
dx x2 (∆x)2

empleando las ecuaciones 1.6.92, 1.6.93 y la ecuación 1.6.94 podemos escribir la siguiente
ecuación matricial para determinar c(x2 )
    
1 0 0 c1 0
 D 2D D   c2  =  0 
(∆x)2 −k − (∆x)2 (∆x)2
0 0 1 c3 1

3. Solución. Usando ∆x = .5 y resolviendo el sistema anterior de ecuaciones

c1 = 0, c2 = .2222, c3 = 1

• ∆x = 13
Usando nuevamente el mismo procediemto que para el caso anterior:
58CHAPTER 1. ECUACIONES DIFERENCIALES ORDINARIAS: PROBLEMAS DE VALORES EN

1. Discretización. La solución empleando la segunda discretización es totalmente analoga

al caso anterior. Los nodos externos estan localizados sobre la frontera del sistema, y se
han definido dos nodos internos equidistantes.
2. Aproximación. De las condiciones frontera, el valor de la solución den los nodos ex-
ternos está fija

c(x1 ) =
c(x4 ) = 1

por lo tanto la ecuación diferencial original únicamente debe aproximarse por diferencias
finitas en los nodos internos x2 y x3
¯ · ¸
d2 c ¯¯ c1 − 2c2 + c3
D 2 ¯ − Kc|x2 = D − Kc2 = 0 (1.6.95)
dx x2 (∆x)2
¯ · ¸
d2 c ¯¯ c2 − 2c3 + c4
D 2 ¯ − Kc|x3 = D − Kc3 = 0 (1.6.96)
dx x3 (∆x)2
o en forma matricial
    
1 0 0 0 c1 0
 D 2D D   c2   0 
−k − (∆x) 0
 (∆x)2 2 (∆x)2    
 0 D
k 2D
− (∆x) D   c3  =  0 
(∆x)2 2 (∆x)2
0 0 0 1 c4 1

1
3. Solución. Resolviendo el sistema anterior de ecuaciones para ∆x = 3

c1 = 0, c2 = .1152, c3 = .3585, c4 = 1

El siguiente ejemplo muestra la forma de resolver problemas de valores en la frontera

cuando alguna de ellas involucra una condición de tipo Neumann.
Ejemplo 13 Difusión con reacción de primer orden-II.
Determine el perfil de concentracción c(x) para el problema enunciado en el ejemplo 12 cuando la
condición frontera en x = 1 es del tipo Neumann.

d2 c
D − Kc = 0 (1.6.97)
dx2
sujeto a las condiciones frontera

c(0) = 0
¯
dc ¯¯
D ¯ = γ
dx x=1

donde γ es una constante. Emplee la discretización sugerida en el inciso i del ejemplo 12.
Solución.

1. Discretización. De acuerdo al enunciado ∆x = .5 y se emplean tres nodos.

2. Aproximación. Este ejemplo es identico al ejemplo 12, excepto por la condición frontera en
x = 1. Dado qur la ubicación de uno de los nodos coincide con el lugar del espacio de solución
donde se especifica la condición frontera de primer orden, esto determina la solución en este
nodo
c1 = 0 (1.6.98)
Diferencias Finitas 59

la aproximación de la ecuación 1.6.97 en el esgundo nodo está dada por

¯ · ¸
d2 c ¯¯ c1 − 2c2 + c3
D 2 ¯ − Kc|x2 = D − Kc2 = 0 (1.6.99)
dx x2 (∆x)2

ya que en el tercer nodo (donde x = 1) está especificada una condición frontera de Neumann,
pueden emplearse diversas formas de obtener la solución c3 en este nodo. Una de estas formas
involucra el plantear la aproximación por diferencias finitas para el tercer nodo empleando un
hipotetico cuarto nodo situado a una distancia ∆x del tercer nodo
¯ · ¸
d2 c ¯¯ c2 − 2c3 + c4
D 2 ¯ − Kc|x3 = D − Kc3 = 0 (1.6.100)
dx x3 2∆x

la condición frontera de tipo Neumann se plantea entonces, en terminos de diferencias finitas,

para el hipotetico cuarto nodo · ¸
c4 − c2
D =γ (1.6.101)
2∆x
despejando c4 de la ecuación 1.6.101
2∆xγ
c4 = + c2
D
substituyendo c4 en la ecuación 1.6.100 eliminamos la ”solución” en el hipotetico cuarto nodo,
y cualquier referencia a este nodo desaparece del problema
· ¸
2c2 − 2c3 2γ
D 2
− Kc3 = − (1.6.102)
(∆x) ∆x

la ecuación 1.6.102 corresponde entonces a la aproximación en terminos de diferencias finitas

en el tercer nodo. De las ecuaciones 1.6.98, 1.6.99 y 1.6.102 se plantea el problema en forma
matricial en los tres nodos
    
1 0 0 c1 0
D 2D D
 (∆x) 2 −K − (∆x) 2 (∆x)2
  c2  =  0 
2D 2D 2γ
0 (∆x)2 −K − (∆x)2
c3 − ∆x

eliminando la primera fila de la ecuación anterior

" #· ¸ · ¸
2D D
−K − (∆x) 2 (∆x)2 c2 0
2D 2D = 2γ
(∆x)2 −K − (∆x) 2 c3 − ∆x

3. Solución. Si γ = .05 la solución del anterior sistema de ecuaciones resulta en los siguientes
valores en los nodos

c2 = .2749, c3 = 1.2329