Captulo III

Solucin de Ecuaciones Diferenciales

Ordinarias

e llama ecuacin diferencial ordinaria a una ecuacin que relaciona la variable

(n)
independiente x, la funcin incgnita y = y(x ) y sus derivadas y, y, ..., y , es decir,

una ecuacin de la forma:

F ( x , y , y, y, ..., y ( n ) ) = 0

(23)

En otras palabras, se llama ecuacin diferencial ordinaria a una ecuacin en la que figura la
derivada o diferencial de una funcin incgnita.
Si la funcin incgnita y = y(x ) depende de una sola variable independiente x, la ecuacin
diferencial se llama ordinaria. Por ejemplo:

dy
+ xy = 0
dx
2 ) y '' + y ' + y = cos x

1)

3) ( x 2 + y 2 )dx + ( x + y )dy = 0

(24)

14 Hidrulica Computacional - V. Yzocupe (2006)

El orden de una ecuacin diferencial es el de la derivada de mayor orden que figura en la

ecuacin. Por ejemplo: la ecuacin diferencial y + xy = e

es de primer orden; la ecuacin

diferencial y + p( x ) y = 0 , donde p(x) es una funcin dada, es de segundo orden.

Se llama solucin de la ecuacin diferencial a una funcin y = (x ) , determinada en el
intervalo (a, b) junto con sus derivadas sucesivas hasta el orden n inclusive, tal que al hacer
la sustitucin y = (x ) en la ecuacin diferencial, sta, se convierte en una identidad con
respecto a x en el intervalo (a, b).
3.1

Ecuaciones Diferenciales Ordinarias de Valor Inicial

Los problemas de ecuaciones diferenciales ordinarias (EDOs) se clasifican en
problemas con condiciones iniciales y problemas con condiciones de frontera. Muchos
de los problemas con condiciones iniciales dependen del tiempo; en ellos, las
condiciones para la solucin estn dadas en el tiempo inicial. Los mtodos numricos
para los problemas con condiciones iniciales difieren en forma significativa de los que
se utilizan para los problemas con condiciones en la frontera.
El problema de una EDO de primer orden con condiciones iniciales se puede escribir
en la forma:

y (t ) = f ( y , t ), y ( 0 ) = y 0

(25)

donde f (y,t) es una funcin de y y t, en tanto que la segunda ecuacin es una

condicin inicial. En la ecuacin anterior, la primera derivada de y est dada como
una funcin conocida de y y t y queremos calcular la funcin incgnita y integrando
numricamente f (y,t). Si f fuera independiente de y, el clculo sera simplemente una
integracin directa. Sin embargo, el hecho que f sea una funcin de la funcin
desconocida y, hace que la integracin sea distinta.
La condicin inicial siempre es parte de la definicin del problema, debido a que la
solucin de un problema con condiciones iniciales slo se puede determinar de
manera nica si dicha condicin inicial est dada.
Algunos ejemplos de problemas con condiciones iniciales de EDO de primer orden:

a ) y ' ( t ) = 3 y + 5 , y( 0 ) = 1

(26)

b ) y ' (t ) = ty + 1, y( 0 ) = 0
1
c ) y ' (t ) =
, y( 0 ) = 1
1 + y2
Existen tres tipos de mtodos de integracin numrica para problemas con
condiciones iniciales: el de Euler, de Runge-Kutta y el de Predictor-corrector.

Solucin de Ecuaciones Diferenciales Ordinarias 15

FORMULA
RELEVANTE

METODOS

Mtodos de Euler:
1. Hacia delante
2. Modificado
3. Hacia atrs

Runge-Kutta:
1. de segundo orden
2. de tercer orden
3. de cuarto orden

Predictor-corrector:
1. de segundo orden
2. de tercer orden
3. de cuarto orden

3.2

ERROR DE TRUNCAMIENTO
LOCAL

GLOBAL

Diferencia hacia delante

Regla del trapecio
Diferencia hacia atrs

0(h2)
0(h3)
0(h2)

0(h)
0(h2)
0(h)

Regla del trapecio

Regla de 1/3 de Simpson
Regla de 1/3 3/8 de
Simpson

0(h3)
0(h4)
0(h5)

0(h2)
0(h3)
0(h4)

Regla del trapecio

Newton hacia atrs
Newton hacia atrs

0(h3)
0(h4)
0(h5)

0(h2)
0(h3)
0(h4)

Mtodos de Euler
Estos mtodos son adecuados para una programacin rpida debido a su sencillez.
Hay que sealar que cuando el sistema de ecuaciones es cada vez ms complicado, se
utilizan con ms frecuencia los mtodos de Euler. Los mtodos de Euler tienen tres
versiones:
a)
b)
c)

Euler hacia adelante

Euler modificado, y
Euler hacia atrs

3.2.1 Mtodo de Euler Hacia Adelante

'

El mtodo de Euler hacia adelante para la ecuacin y = f ( y, t ) se obtiene de la

aproximacin por diferencias hacia adelante,

yn +1 yn
y n' y n + 1 = y n + t y n'
t
'

reemplazando y n = f ( y n , t n ) , se obtiene:

y n + 1 = y n + t f ( y n , tn )
Mediante la ecuacin anterior se calcula y n en forma recursiva como

(27)

16 Hidrulica Computacional - V. Yzocupe (2006)

(28)

y1 = y 0 + t y 0' = y 0 + t f ( y 0 , t0 )
y 2 = y1 + t f ( y1 , t1 )
y3 = y 2 + t f ( y 2 , t2 )

y n = y n 1 + t f ( y n 1 , tn 1 )

Aunque el mtodo de Euler hacia adelante es muy sencillo, debe utilizarse

cuidadosamente para evitar dos tipos de errores. El primer tipo lo forman los
errores de truncamiento, y el segundo tipo lo constituye la posible inestabilidad,
que aparece cuando la constante del tiempo es negativa (la solucin tiende a
cero sino hay trmino fuente), a menos que el intervalo de tiempo h sea
suficientemente pequeo.
3.2.2 Mtodo de Euler Modificado
El mtodo de Euler modificado tiene dos motivaciones. La primera es que es
ms preciso que el anterior. La segunda es que es ms estable.
'

Este mtodo se obtiene al aplicar la regla del trapecio para integrar y = f ( y , t )

yn +1 = yn +

t
[ f ( y n + 1 , t n + 1 ) + f ( y n , t n )]
2

(29)

3.2.3 Mtodo de Euler Hacia Atrs

Este mtodo se basa en la aproximacin por diferencias hacia atrs, la cual se
escribe como:

y n y n 1
y n' y n = y n 1 + t y n'
t
'

reemplazando y n = f ( y n , t n ) y generalizando para n +1, se obtiene:

y n + 1 = y n + t f (y n + 1 , t n + 1 )

(30)

La precisin de este mtodo es la misma que la de Euler hacia adelante. Adems

si f es una funcin no lineal de y, debe utilizarse un esquema iterativo en cada
paso. Sin embargo, las ventajas son: a) el mtodo es estable para los problemas
rgidos, y b) la solucin es positiva si la solucin exacta es positiva.

Solucin de Ecuaciones Diferenciales Ordinarias 17

3.2.4 Ejemplo de Aplicacin del Mtodo de Euler

Un tanque cilndrico de fondo plano con un dimetro de 1.5 m contiene un
liquido de densidad 1,500 kg/m 3 a una altura h0=3 m. Se desea conocer la altura
del cilindro dentro del tanque 3 minutos despus de que se abre completamente
la vlvula de salida, la cual proporciona un caudal de 0.6 a 2 gh m3/s, donde a
es el rea transversal del tubo de salida; a = 78.5 x 10-4 m2, g = 9.81 m/s2.

SOLUCION:
Aplicando la ley de conservacin de masa dentro del tanque:
Acumulacin = Entrada - Salida

d
( ) = 0 Qs ,
= 4 D 2 h
dt
d
D 2 h = 0.6 a 2 g h
dt 4
dh

D2
= 0.6 a 2 g h
4
dt
2.4 a 2 g
2.4 78.5 x10 - 4 2 x 9.81
dh
=
h
=

h
dt
(1.5 )2
D2

dh
= 0.0118 h
dt

Planteamiento del PVI:

dh
= h = 0.0118 h ,
dt

h (t = 0 s ) = 3 m,

h (t = 180 s ) = ?

Aplicaremos el mtodo de Euler hacia delante para un t=30 s,

hn + 1 = hn + t hn = hn + (30 s ) 0.0118 hn
hn + 1 = hn 0.354 hn

18 Hidrulica Computacional - V. Yzocupe (2006)

Tabulando los resultados se tiene:

3.3

t(s)

30

60

90

120

150

180

h(m)

3.00

2.39

1.84

1.36

0.95

0.60

0.33

Mtodos de Runge-Kutta
Una desventaja fundamental de los mtodos de Euler consiste en que sus rdenes de
precisin son bajos. Esta desventaja tiene dos facetas. Para mantener una alta
precisin se necesita un t pequeo, lo que aumenta el tiempo de clculo y provoca
errores de redondeo.
En los mtodos de Runge-kutta, el orden de precisin aumenta al utilizar puntos
intermedios en cada intervalo. Una mayor precisin implica adems que los errores
decrecen ms rpido al reducir t, en comparacin con los mtodos de precisin baja.
Consideremos una ecuacin diferencial ordinaria

y' =

dy
= f ( y , t ), y ( 0 ) = y 0
dt

(31)

Para calcular y n +1 en tn +1 = tn + t , dado un valor de y n , integramos la ecuacin

anterior en el intervalo [tn , tn + 1 ] :
y n +1

y ]yy n+1 =
n

t n +1

dy =

t n +1

ydt

(32)

f (y , t )dt

yn +1 = yn +

t n +1

tn

f (y , t )dt

Los mtodos de Runge-Kutta se obtienen al aplicar un mtodo de integracin

numrica a la integral del lado derecho de la ecuacin (32).
3.3.1 Mtodo de Runge-Kutta de Segundo Orden
Aplicamos la regla del trapecio al lado derecho de la ecuacin (32):
t n +1

f (y , t ) dt

1
t [ f ( y n , tn ) + f ( y n + 1 , tn + 1 ) ]
2

(33)

En esta ecuacin y n + 1 es una incgnita, por lo que aproximamos el segundo

trmino mediante f ( y n +1 , tn +1 ) , donde y n +1 es la primera estimacin de y n +1

Solucin de Ecuaciones Diferenciales Ordinarias 19

obtenida mediante el mtodo de Euler hacia delante. Este esquema se conoce

como el mtodo de Runge-Kutta de segundo orden y se resume como:

y n + 1 = y n + t f ( y n , tn )
t
[ f ( yn , tn ) + f ( yn + 1 , tn + 1 )]
yn +1 = yn +
2

(34)

En forma cannica lo anterior es:

k1 = t f ( y n , tn )
k 2 = t f ( y n + k1 , tn + 1 )
1
y n + 1 = y n + [ k1 + k 2 ]
2

(35)

3.3.2 Mtodo de Runge-Kutta de Tercer Orden

Un mtodo de Runge-Kutta ms preciso que el anterior es resultado de un
esquema de integracin numrica de orden superior para el segundo trmino
de la ecuacin (32) es:

yn +1 = yn +

t
f ( y n , tn ) + 4 f ( y n + 1 , tn + 1 ) + f ( y n + 1 , tn + 1 )
2
2
6

(36)

Donde y n + 1 y y n + 1 son estimaciones, puesto que no conocemos y n + 1 ni y n + 1 .

2

Se obtiene la estimacin y n + 1 mediante el mtodo de Euler hacia delante:

2

t
f ( y n , tn )
2

yn + 1 = yn +
2

(37)

La estimacin y n + 1 es: y n + 1 = y n + t f ( y n , tn )
o bien

y n + 1 = y n + t f ( y n + 1 , tn + 1 )
2

(38)

o una combinacin lineal de ambas

y n + 1 = y n + t f ( y n , tn ) + (1 ) f ( y n + 1 , tn + 1 )
2

(39)

Donde es un parmetro que hay que determinar de forma que maximice la

precisin del mtodo numrico. En resumen, el mtodo de Runge-Kutta con
una precisin de tercer orden es:

k1 = t f ( y n , t n )
k 2 = t f ( y n + 1 k1 , t n + t )
2

k3 = t f ( y n k1 + 2k 2 , t n + t )
y n + 1 = y n + 1 [ k1 + 4k 2 + k3 ]
6

(40)

20 Hidrulica Computacional - V. Yzocupe (2006)

3.3.3 Mtodo de Runge-Kutta de Cuarto Orden

El mtodo de Runge-Kutta de cuarto orden se obtiene de una manera similar al
de tercer orden, excepto que se utiliza un paso intermedio adicional para
evaluar la derivada. El mtodo de Runge-Kutta de cuarto orden tiene una
precisin hasta el trmino de cuarto orden del desarrollo de Taylor, por lo que
el error local es proporcional a h5.
La siguiente es la versin del mtodo de Runge-Kutta de cuarto orden basada
en la regla de 1/3 de Simpson para integrar numricamente la ecuacin (32):

k1 = t f ( y n , t n )
k 2 = t f ( yn +
k3 = t f ( y n +

k1
, t + t )
2
2 n
k2
, t n + t )
2
2

(41)

k 4 = t f ( y n + k3 , t n + t )

y n + 1 = y n + 1 [ k1 + 2k 2 + 2k3 + k 4 ]
6

La segunda versin se basa en la regla de 3/8 de Simpson:

k1 = t f ( y n , t n )
k 2 = t f ( y n +
k3 = t f ( y n +

k1
3
k1
3

(42)

, t n + t )
3

k2
3

, t n + 2t )
3

k 4 = t f ( y n + k1 k 2 + k3 , t n + t )
y n + 1 = y n + 1 [ k1 + 3k 2 + 3k3 + k 4 ]
8

3.3.4 Ejemplo de Aplicacin del Mtodo Runge-Kutta IV

Calcule el tiempo necesario para que

el nivel del lquido cambie de 4 m a 3
m dentro del tanque esfrico de
radio R0=5 m. La velocidad de salida
por

el

orificio

del

fondo

es

u = 4.895 h (m/s), el dimetro de

dicho orificio, d = 10 cm.

Solucin de Ecuaciones Diferenciales Ordinarias 21

SOLUCION:
Aplicando la ley de conservacin de masa dentro del tanque:
Acumulacin = Entrada - Salida

d
( ) = 0 A u ,
= Ro h 2 h 3 3 , A = d 2
dt
4

d
R o h 2 h 3 3 = d 2 4.895 h
dt
4
dh
dh
4.895
d2 h
h 2
=
2 R o h
dt
dt
4
4.895
2 dh
2 R o h h
d2 h
=
dt
4
2
h
4.895 d
dh
=
2
dt
4 2Ro h h

[ (

)]

dh
0.01224 h
=
2
dt
10 h h
Planteamiento del PVI:
0.01224 h
dh
= h =
,
2
dt
10 h h

h (t = 0 s ) = 4 m,

h (t = ? s ) = 3 m

Aplicaremos el mtodo de Runge-Kutta IV para un t=50 s,

k1 = 0.01224 t hn 10 hn hn2

(hn + k / 2 ) [10(hn + k / 2 ) (hn + k / 2 )2 ]

k 3 = 0.01224 t (hn + k / 2 ) [10(hn + k / 2 ) (hn + k / 2 )2 ]
k 4 = 0.01224 t (hn + k ) [10(hn + k ) (hn + k )2 ]
hn + 1 = hn + 61 (k1 + 2 k 2 + 2 k 3 + k 4 )
k 2 = 0.01224 t

Tabulando los resultados se tiene:

n

t (s)

hn(m)

k1

k2

k3

k4

hn+1

4.000

-0.05100000

-0.05094682

-0.05094687

-0.05089576

3.949

50

3.949

-0.05089576

-0.05084670

-0.05084675

-0.05079974

3.898

100

3.898

-0.05079974

-0.05075478

-0.05075482

-0.05071189

3.847

150

3.847

-0.05071189

-0.05067100

-0.05067103

-0.05063217

3.797

200

3.797

-0.05063217

-0.05059533

-0.05059536

-0.05056055

3.746

250

3.746

-0.05056055

-0.05052776

-0.05052778

-0.05049702

3.696

300

3.696

-0.05049702

-0.05046828

-0.05046830

-0.05044159

3.645

350

3.645

-0.05044159

-0.05041691

-0.05041692

-0.05039427

3.595

400

3.595

-0.05039427

-0.05037366

-0.05037366

-0.05035510

3.544

450

3.544

-0.05035509

-0.05033858

-0.05033858

-0.05032412

3.494

10

500

3.494

-0.05032412

-0.05031172

-0.05031172

-0.05030140

3.444

11

550

3.444

-0.05030140

-0.05029317

-0.05029317

-0.05028703

3.393

22 Hidrulica Computacional - V. Yzocupe (2006)

12

600

3.393

-0.05028703

-0.05028300

-0.05028300

-0.05028109

3.343

13

650

3.343

-0.05028109

-0.05028133

-0.05028133

-0.05028372

3.293

14

700

3.293

-0.05028372

-0.05028828

-0.05028828

-0.05029503

3.243

15

750

3.243

-0.05029503

-0.05030400

-0.05030400

-0.05031520

3.192

16

800

3.192

-0.05031520

-0.05032865

-0.05032866

-0.05034439

3.142

17

850

3.142

-0.05034439

-0.05036242

-0.05036243

-0.05038280

3.092

18

900

3.092

-0.05038280

-0.05040553

-0.05040554

-0.05043066

3.041

19

950

3.041

-0.05043066

-0.05045820

-0.05045821

-0.05048821

2.991

20

1000

2.991

-0.05048821

-0.05052069

-0.05052072

-0.05055573

2.940

300

500

700

900

4.5

h (m)

4.0

3.5

3.0

2.5
0

100

200

400

600

800

1000

Tiempo (seg)

3.4

Solucin de Sistemas de EDOs

La aplicacin del mtodo de Runge-Kutta de cuarto orden a un conjunto de
ecuaciones diferenciales ordinarias es simple. Para simplificar la explicacin,
consideremos un conjunto de ecuaciones:

y ' = f ( y, z, t )

(43)

'

z = g( y , z, t )
El mtodo de Runge-Kutta de cuarto orden para este conjunto es:

k1 = t f ( y n , z n , t n ), l1 = t g ( y n , z n , t n )
k 2 = t f ( y n +
k3 = t f ( y n +

k1
l
, z + 1 , t + t ),
2 n 2 n
2
k2
l2
, z n + , t n + t ),
2
2
2

l2 = t g ( y n +
l3 = t g ( y n +

k1
l
, z + 1 , t + t )
2 n 2 n
2
k2
l2
, z n + , t n + t )
2
2
2

k 4 = t f ( y n + k3 , z n + l3 , t n + t ), l 4 = t g ( y n + k3 , z n + l3 , t n + t )
y n + 1 = y n + 1 [ k1 + 2k 2 + 2k3 + k 4 ]
6

z n + 1 = z n + 1 [l1 + 2l 2 + 2l3 + l 4 ]
6

(44)

Solucin de Ecuaciones Diferenciales Ordinarias 23

3.4.1 Ejemplo de Solucin de un Sistema de EDOs

Considere un sistema ecolgico simple compuesto de coyotes (y) y correcaminos
(x), donde los primeros se alimentan de los segundos. Los tamaos de las
poblaciones cambian de acuerdo con las ecuaciones: ( dx dt ) = c 1 x c 2 xy ,

( dy dt ) = c 3 xy c 4 y , donde: c1 = 0.4 , c 2 = 0.02 , c 3 = 0.001, c 4 = 0.3 y

x (t = 0 ) = 30 , y (t = 0 ) = 3 .
Si no hay coyotes (y) los correcaminos se reproducen con una velocidad de
crecimiento c1x; si no hay correcaminos, la especie de coyotes desaparece con
una velocidad c4y. El trmino xy representa la interaccin de las dos especies y
las constantes c2 y c3 dependen de la habilidad de los depredadores para atrapar
a los correcaminos y de la habilidad de stos para huir. La poblacin de los
coyotes cambia cclicamente. Calcule el ciclo y su periodo.
SOLUCION:
Planteamiento del PVI:

dx
= c1 x c 2 xy
dt
dy
= c 3 xy c 4 y
dt

c1 = 0.4 , c 2 = 0.02 , c 3 = 0.006 , c 4 = 0.3

donde;
x (t = 0 ) = 30 , y (t = 0 ) = 3
Aplicamos el mtodo de Runge-Kutta IV para un sistema de EDOs con un t=2
das.

k1 = t (c1 xn c2 xn y n )
l1 = t (c3 xn y n c4 y n )
k
k
l
k 2 = t c1 ( xn + 1 ) c2 ( xn + 1 )( y n + 1 )
2
2
2

k
l
l
l2 = t c3 ( xn + 1 )( y n + 1 ) c4 ( y n + 1 )
2
2
2

k
k
l
k3 = t c1 ( xn + 2 ) c2 ( xn + 2 )( y n + 2 )
2
2
2

k
l
l
l3 = t c3 ( xn + 2 )( y n + 2 ) c4 ( y n + 2 )
2
2
2

k 4 = t [c1 ( xn + k3 ) c2 ( xn + k3 )( y n + l3 )]

l4 = t [c3 ( xn + k3 )( y n + l3 ) c4 ( y n + l3 )]
xn +1 = xn + 1 [ k1 + 2k 2 + 2k3 + k 4 ]
y n +1 =

6
y n + 1 l1 + 2l 2
6

+ 2l3 + l 4 ]

24 Hidrulica Computacional - V. Yzocupe (2006)

Tabulando los resultados se tiene:

n

t(das)

xn

yn

k1

l1

k2

l2

30

20.40

-1.08

28.20

60

42.51

-0.28

57.98

122

84.85

0.96

234

134.38

7.23

301

30

-125.35

55.03

10

66

73

-139.23

-5.16

12

23

48

-25.63

14

11

29

-4.26

16

10

18

18

13

10

20

11

22

12

k3

l3

k4

l4

-0.68

30.59

0.11

62.80

111.27

2.18

115.90

3.00

129.24

16.84

70.68

25.26

-361.85

75.79

-232.47

24.77

6.98

-43.99

-85.33

-2.18

-19.67

-7.40

-19.47

-13.93

-15.05

-0.71

-11.58

-1.60

0.97

-9.14

2.92

-6.71

2.64

11

4.82

-5.24

7.32

-3.78

20

10.89

-2.87

15.30

36

22.55

-1.39

30.86

24

69

44.92

-0.18

13

26

131

81.83

1.99

14

28

226

10

90.14

15

30

218

38

-155.45

16

32

55

60

17

34

18

18

36

19

38

20

xn+1

yn+1

-0.66

42.80

-0.27

60

0.27

85.58

1.00

122

137.11

7.25

234

-133.43

56.95

301

30

-96.77

-2.74

66

73

38.75

-53.10

23

48

-16.76

-4.02

-16.17

11

29

-12.16

1.03

-9.06

10

18

-7.23

4.85

-5.16

13

11

7.45

-4.05

10.91

-2.84

20

-2.02

16.13

-2.09

22.61

-1.38

36

-0.83

32.98

-0.76

45.11

-0.17

69

59.90

0.47

63.52

0.75

82.38

2.04

131

100.49

4.20

97.90

5.55

91.98

12.29

226

10

12.16

42.08

25.34

-26.84

31.32

-170.50

41.20

218

38

43.28

-221.48

31.07

-143.16

13.78

-94.66

-0.05

55

60

-87.63

-9.69

-15.43

-28.24

-48.70

-10.24

-7.25

-27.42

18

41

41

-14.77

-18.81

-4.73

-16.36

-7.80

-15.68

-2.08

-13.19

11

25

11

25

-2.11

-12.85

0.56

-9.70

-0.05

-10.31

2.26

-7.55

11

15

11

15

2.21

-7.64

4.26

-5.58

4.10

-6.00

6.63

-4.27

15

40

15

6.61

-4.33

9.66

-3.11

9.97

-3.30

14.29

-2.29

25

21

42

25

14.26

-2.31

19.81

-1.59

21.00

-1.60

29.15

-0.97

46

22

44

46

29.06

-0.97

39.43

-0.43

42.09

-0.30

56.73

0.41

87

23

46

87

56.47

0.38

73.82

1.31

77.04

1.78

94.53

4.03

163

24

48

163

93.70

3.95

104.09

8.18

88.47

10.87

35.34

23.26

249

17

25

50

249

17

34.78

22.91

-84.67

42.71

-147.32

39.75

-146.81

11.81

153

50

26

52

153

50

-181.72

30.89

-111.62

-6.93

-101.82

8.06

-76.80

-11.22

38

53

27

54

38

53

-51.26

-15.67

-13.01

-22.72

-28.10

-14.53

-7.74

-20.16

15

35

28

56

15

35

-8.90

-16.88

-2.73

-13.76

-4.39

-13.87

-0.49

-10.94

11

21

29

58

11

21

-0.52

-10.88

1.80

-8.12

1.35

-8.67

3.65

-6.29

12

13

30

60

12

13

3.61

-6.38

5.95

-4.64

5.91

-4.97

8.97

-3.52

18

350
300

Poblacin

250
200
150
100
50
0
0

12

16

20

24

28

32

Tiempo (das)

36

40

44

48

Correcaminos

52

56

60

Coyotes

Solucin de Ecuaciones Diferenciales Ordinarias 25

3.5

Tarea de Aplicacin del Mtodo de Runge-Kutta IV

A un tanque cnico de 6 m de altura y 3 m de radio, que contiene un liquido a
una altura h0=5 m, se le hace llegar un caudal de 7 l/s. Treinta minutos despus
este flujo se interrumpe por falla de la bomba, la cual se repara y arranca una
hora despus. Determine el caudal necesario para que el nivel se recupere y se
mantenga en 5 m; as como el tiempo necesario para alcanzar ese nivel (rgimen
permanente). El caudal de salida es igual a 3.457 h l/s ininterrumpidamente.

26 Hidrulica Computacional - V. Yzocupe (2006)

3.6

Ecuaciones Diferenciales Ordinarias de Valor Frontera

En los problemas de ecuaciones diferenciales ordinarias unidimensionales con valores
en la frontera, se pide que la solucin satisfaga las condiciones de frontera en ambos
extremos del dominio. La definicin de las condiciones de frontera es parte
fundamental de los problemas de este tipo. Por ejemplo, consideremos una varilla
delgada de metal de longitud L, tal que sus extremos estn conectados a distintas
fuentes de calor. Si el calor sale de la superficie de la varilla nicamente mediante la
transferencia de calor por conveccin, se puede escribir la ecuacin gobernante del
proceso fsico como:

d
d
k (x ) T (x ) + h c PT (x ) = h c PT + AS(x )
dx
dx

(45)

donde;
T(x), es la temperatura en cualquier posicin x desde el extremo izquierdo
A, es el rea constante de una seccin transversal de la varilla
k, es la conductividad trmica,
P, es el permetro de la varilla,
hc, es el coeficiente de transferencia de calor por conveccin,
T , es la temperatura neta del aire, y
S, es la fuente de calor.
Las condiciones de frontera son:

T (x = 0 ) = Ti
T (x = L ) = Td

(46)

donde Ti y Td son respectivamente las temperaturas del cuerpo en los extremos

izquierdo y derecho, respectivamente.
3.6.1 Tipos de Condiciones de Frontera
Una ecuacin diferencial de segundo orden puede resolverse analticamente,
obtenindose dos constantes de integracin. Se deben tener dos condiciones de
frontera para poder resolver las dos constantes de integracin. Usualmente,
ellas estn especificadas en los extremos del dominio unidimensional del
problema.

d 2T
+Q=0
dx 2

(47)

Fsicamente, las condiciones de frontera se clasifican como esenciales y

naturales. Matemticamente, una condicin de frontera esencial se conoce como
condicin de Dirichlet, y una condicin de frontera natural se llama condicin
de Neumann.

Solucin de Ecuaciones Diferenciales Ordinarias 27

El tercer tipo de condicin de frontera es el llamado mixto, y corresponde a una

combinacin de los dos tipos anteriores. En realidad existen dos tipos de
condiciones de frontera mixta. La primera es cuando una condicin esencial se
especifica en una frontera y una condicin natural en la otra frontera. El
segundo tipo de condicin de frontera mixta se da cuando se combinan los tipos
esenciales y naturales en una sola condicin de frontera.
El siguiente cuadro muestra ejemplos de condiciones de frontera relacionados
con la ecuacin diferencial de segundo orden (47):

TIPO

DESCRIPCION

Esencial/Dirichlet
Condicin en la frontera
con un valor fijo

Se da un valor de
la funcin solucin

Natural/Neumann
Condicin en la frontera
con una derivada

Se da un valor de la
derivada de la funcin
solucin

Mixta/Cauchy
Condicin en la frontera
de tipo mixto

Se relaciona un valor de
la funcin con la derivada

EJEMPLOS

T( x = 0 ) = T0
T( x = L ) = TL

dT ( x = 0 )
= q0
dx
dT ( x = L )
kL
= qL
dx
k0

dT
+ h (TT ) = 0
dx

3.6.2 Mtodo de Solucin de EDOs con Valor de Frontera

En esta seccin, aplicaremos el mtodo de diferencias finitas para resolver los
problemas de ecuaciones diferenciales ordinarias de segundo orden con
condiciones de frontera.
Para explicar el principio del mtodo, consideremos la ecuacin:

d 2 y (x )
dx 2

+ qy (x ) = S(x ),

0xL

(48)

con las siguientes condiciones de frontera:

dy
(x = 0 ) = 0
dx
y (x = L ) = y D

( frontera izquierda )

(49)

( frontera derecha )

donde, q es un coeficiente constante.

Si dividimos el dominio en N intervalos de igual longitud (segmentos),
obtendremos una retcula donde los intervalos miden h=L/N.

28 Hidrulica Computacional - V. Yzocupe (2006)

Figura 6. Retcula unidimensional

3.6.2.a Algoritmo para los Nodos Internos

Si aplicamos la aproximacin por diferencias centrales (6) al primer

trmino de la ecuacin (48), con lo que se obtiene la ecuacin en
diferencias para la i-sima retcula:

y 2 y i + y i +1
i 1
+ qy i = Si
x 2

(50)

donde yi = y (xi ) , Si =S(xi ) .

yi 1 + 2 + q x 2 yi yi +1 = x 2 S i

(51)

Este algoritmo para nodos internos, se utiliza para todos los puntos de la
retcula, excepto cuando i=1 e i=N+1.
3.6.2.b Algoritmo para la Frontera Izquierda

La condicin de la frontera izquierda dada por la ecuacin (49) es del tipo

Neumann, para su implementacin procedemos de la siguiente manera;
utilizamos una diferencia finita hacia atrs para la derivada especificada
en el punto de la grilla i=1 localizado en x=0.

y1 y0
=0
x

dy
= 0;
dx

y1 = y0

(52)

Luego, la ecuacin (51) para el caso i = 1 ser:

y0 + 2 + q x 2 y1 y 2 = x 2 S1

(53)

Utilizando el resultado de la ecuacin (52), se obtiene el algoritmo para la

frontera izquierda:

(1 + qx 2 ) y1 y2 = x 2 S1

(54)

3.6.2.c Algoritmo para la Frontera Derecha

Escribimos la ecuacin (51) para i=N, conjuntamente con la condicin de

frontera y N +1 = y D para obtener el algoritmo de la frontera derecha:

Solucin de Ecuaciones Diferenciales Ordinarias 29

y N 1 + 2 + q x 2 y N = x 2 S N + y D

(55)

El conjunto de ecuaciones (51), (54) y (55) se escriben en forma conjunta como:

(56)

(1 + qx 2 )y1

y2
+ 2 + q x 2 y2
y2

y1

y3
+ 2 + q x 2 y3

y N 1

= x 2 S1
= x 2 S 2
= x 2 S3

y4

+ 2 + q x y N

= x S N + y D
2

Forman lo que se conoce como un sistema de ecuaciones algebraicas lineales, el

mismo que arreglado en forma matricial como:

1 + q x 2

1
2 + q x 2
1

1
2 + q x 2

2 + q x 2

x 2 S1
(57)

2
x S 2
= x 2 S3

y N x S N + y D
y1
y2
y3

[ A ]{y }= [ D ]
El sistema se compone de una matriz de coeficientes A, el vector de incgnitas y
y el vector de residuos D. La matriz A se conoce como matriz tridiagonal, porque
todos los elementos de la matriz (57) son cero, excepto los elementos de las tres
diagonales principales.
3.6.3 Algoritmo de Solucin para Sistemas Tridiagonales

B1
A
2

C1
B2
A3

C2
B3

C3
An

Bn

y1 D1
y D
2 2
y3 = D3

y N DN

(58)

El algoritmo de solucin para este tipo de sistemas recibe el nombre de solucin

tridiagonal o doble barrido, y se compone de los siguientes pasos:
a) Se inicializan dos nuevas variables:

B1' = B1

D1' = D1

(59)

30 Hidrulica Computacional - V. Yzocupe (2006)

b)

Se calculan en forma recursiva las siguientes ecuaciones, en orden creciente

de i hasta llegar a N:

R = Ai Bi' 1

(60)

Bi' = Bi RC i 1
Di' = Di RDi' 1
c)

para i = 2 ,3,..., N

Se calcula la solucin para la ltima incgnita:

'
'
y N = DN
BN

d)

Se calculan las soluciones de las dems incgnitas en orden decreciente de

i:
(62)
y i = Di' Ci y i +1 Bi'
i = N1, N2 , ..., 1

3.7

(61)

Ejemplo de Solucin de una EDO de Valor Frontera

La temperatura permanente, T, en una aleta disipadora de calor puede ser simulada
mediante el siguiente modelo lineal:

d 2T
P T = P Tb ,
dx 2

0xL

y sus condiciones de frontera:

T ( x = 0) = Ta
dT
h
( x = L ) = [Tb T ( x = L ) ]
dx
k
donde;

L = 9.126 m,
o

Ta = 503 C ,
o

Tb = 91 C ,
x = 0.3042 m,

P=

h(2t + 2 w)
,
12 kwt

h = 3.9
k = 30
w = 1.465
t = 0.228

Encuentre la variacin de la temperatura, T, a lo largo de la distancia x.

SOLUCION:

(63)

Solucin de Ecuaciones Diferenciales Ordinarias 31

La condicin de frontera izquierda establece un valor para la variable temperatura, T,

en el punto x=0, i=0, por lo que se tiene:
(64)

Ti =0 = Ta

Para deducir el algoritmo general de los nodos internos, utilizamos la ecuacin (6)
para reemplazar la derivada de segundo orden en la ecuacin gobernante del proceso
fsico (63), obteniendo de esta manera el algoritmo general a ser aplicado en los nodos
i=1,, 29:
Ti 1 2Ti + Ti + 1

P Ti = P Tb
x 2
Ti 1 2Ti + Ti + 1 x 2 P Ti = x 2 P Tb

Ti 1 ( 2 + x 2 P ) Ti + Ti + 1 = x 2 P Tb

(65)

La condicin de frontera derecha establece una condicin tipo mixta para la variable
temperatura, T, en el punto x=L, i=30. El tratamiento de este tipo de condicin es
diferente:

dT
dT
h
( x = L) =
= [Tb T (i = 30)]

dx
dx i = 30 k

(66)

Primero, en la ecuacin (63) reemplazamos la derivada de segundo orden por dos

derivadas de primer orden con base en i=30 y el intervalo de longitud x 2 :

dT
dT
dx i = 29 1
dx i = 30
2
x
2

P Ti = 30 = P Tb

donde podemos sustituir T )i = 29 1 con una aproximacin en diferencias centrales

2

dT
dx

i = 29 12

Ti = 30 Ti = 29
x

y la ecuacin (66) se utiliza para reemplazar a T )i = 30

Ti = 30 Ti = 29
x

hk (Tb Ti = 30 )
x
2

P Ti = 30 = P Tb

2
2

Ti = 29 + 1 + x h x P Ti = 30 = x h + x P Tb
k
k
2
2

(67)

Luego, el conjunto tridiagonal esta formado por las siguientes ecuaciones algebraicas:

T0 = Ta
Con r = x hk y s = x 2 P , aplicamos la ecuacin (65) para i=1,,29:

i =1 :

T0 ( 2 + s)T1 + T2 = sTb

i =2 :
i =3 :

T1 ( 2 + s)T2 + T3 = sTb
T2 ( 2 + s)T3 + T4 = sTb

i = 28 :
i = 29 :

T27 ( 2 + s)T28 + T29 = sTb

T28 ( 2 + s)T29 + T30 = sTb

32 Hidrulica Computacional - V. Yzocupe (2006)

y la ecuacin (67) para i=30:

( )

T29 + 1 + r s T30 = r + s Tb
2

Las ecuaciones obtenidas pasan a formar parte del sistema matricial AT=b, el cual fue
solucionado mediante el mtodo de eliminacin de Gauss, obtenindose el vector
solucin T. Dicho procedimiento se muestra en el programa de cmputo siguiente:

Program temperatura
!
real coef(31,31), vctb(31)
real A(31,31), x(31), b(31), k
real*8 suma, uc
integer i, j, ii, nn
! Abre archivo de salida
open (2,file='temperatura.out')
! Establece valores de parametros de control
ii=31; nn=31; dx=0.3042
! Asigna datos y calcula coeficientes
Ta=503; Tb=72; h=3.9; k=50; w=1.465; t=0.228
P = h*(2*t+2*w)/(12*k*w*t); r = dx*h/k; s = dx*dx*P
! Limpia la matriz de coeficientes y el vector de residuos
coef = 0; vctb = 0
! Coloca los coeficientes de la condicion de frontera izquierda
coef(1,1) = 1; vctb(1) = Ta
! Coloca los coeficientes de los nodos internos
do i=2,ii-1
coef(i,i-1) = 1
coef(i,i) = -(2+s)
coef(i,i+1) = 1
vctb(i)
= s*Tb
enddo
! Coloca los coeficientes de la condicion de frontera derecha
coef(ii,ii-1) = -1
coef(ii,ii) = 1+r-0.5*s
vctb(ii)
= (r+0.5*s) *Tb
! Transfiere los coeficientes a la matriz A y vector B
do i=1,ii
do j=1,ii
A(i,j) = coef(i,j)
enddo
b(i) = vctb(i)
enddo
! Imprime la matriz de coeficientes y el vector de residuos
write (2,20)
do i=1,nn
write (2,21) ( A(i,j),j=1,nn ), b(i)
enddo
! Realiza la solucin por el mtodo de Gauss
do i=1,nn-1

Solucin de Ecuaciones Diferenciales Ordinarias 33

do k=i+1,nn
if( A(k,i)==0.) cycle
uc = A(k,i)/A(i,i)
do j=1,nn
A(k,j) = A(k,j) - uc*A(i,j)
enddo
b(k) = b(k) - uc*b(i)
enddo
enddo
do i=nn,1,-1
suma = 0.
do j=i+1,nn
suma = suma + A(i,j)*x(j)
enddo
x(i) = ( b(i)-suma )/A(i,i)
enddo
! Imprime la solucin desde la frontera izquierda hacia la derecha
write (2,23)
write (2,26) ( j-1, x(j),j=1,nn )
20
21
23
26

format (/15x, 'Matriz A :',20x,'Vector b :'/)

format (31f6.0, 8x, f10.3)
format (/10x, 'Solucin por el Mtodo de Gauss'/)
format (31(i5,4x,f10.3/))
close (2)
stop
end

0.439
Matriz A :
1. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
1. -2. 1. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
0. 1. -2. 1. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
0. 0. 1. -2. 1. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
0. 0. 0. 1. -2. 1. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 1. -2. 1. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 1. -2. 1. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0. 1. -2. 1. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0. 0. 1. -2. 1. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0. 0. 0. 1. -2. 1. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0. 0. 0. 0. 1. -2. 1. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 1. -2. 1. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 1. -2. 1. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 1. -2. 1. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 1. -2. 1. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 1. -2. 1. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 1. -2. 1. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 1. -2. 1. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 1. -2. 1. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 1. -2. 1. 0. 0. 0. 0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 1. -2. 1. 0. 0. 0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 1. -2. 1. 0. 0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 1. -2. 1. 0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 1. -2. 1. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 1. -2. 1. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 1. -2. 1. 0. 0. 0.
0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 1. -2. 1. 0. 0.
0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 1. -2. 1. 0.
0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 1. -2. 1.
0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 1. -2.
0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. -1.

Vector b:
0.
0.
0.
0.
0.
0.
0.
0.
0.
0.
0.
0.
0.
0.
0.
0.
0.
0.
0.
0.
0.
0.
0.
0.
0.
0.
0.
0.
0.
1.
1.

503.000
.439
.439
.439
.439
.439
.439
.439
.439
.439
.439
.439
.439
.439
.439
.439
.439
.439
.439
.439
.439
.439
.439
.439
.439
.439
.439
.439
.439
.439
1.928

34 Hidrulica Computacional - V. Yzocupe (2006)

Vector Solucin (Mtodo de Gauss):

TEMPERATURA
600
500

T (C)

400
300
200
100
0
0

10

12

14

16

X (m)
i
0
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30

T
503.00
460.874
421.998
386.133
353.062
322.583
294.509
268.671
244.909
223.08
203.05
184.697
167.909
152.584
138.628
125.957
114.493
104.165
94.912
86.677
79.409
73.064
67.604
62.995
59.209
56.224
54.02
52.584
51.908
51.988
52.823

18

20

22

24

26

28

30