Alex W.

Pilco Nuñez, Facultad de Ingeniería Química y Textil, Universidad Nacional de Ingeniería

CAPÍTULO 9

ECUACIONES DIFERENCIALES ORDINARIAS

Objetivos

➢ Estudiar los métodos numéricos de Euler, Taylor y Runge – Kutta de cuarto

orden para resolver una ecuación diferencial ordinaria.
➢ Comparar los resultados que se obtienen usando e los métodos numéricos
con los resultados analíticos en el caso que se pueda resolver.

Cuando estudiamos cálculo diferencial e integral y ecuaciones diferenciales,

resolvemos modelos matemáticos que se utilizan en muchos problemas de las
ciencias, como ingeniería, economía y otras disciplinas.

En ese capítulo se busca exponer una solución aproximada por métodos

numéricos, si es que existen condiciones iniciales del problema.

Hay muchos programas en diferentes paquetes de software que están diseñados

para resolver cualquier tipo de ecuación diferencial que se presente.

Formulación del problema de valor inicial (PVI)

La ecuación diferencial ordinaria (EDO) general de primer orden es

= f (x, y )
dy
(8.1)
dx

En la teoría de las EDO se establece que su solución general debe contener una
constante arbitraria c , de tal modo que la solución general de la ecuación 8.1 es

F ( x, y , c ) = 0 (8.2)

1
Alex W.Pilco Nuñez, Facultad de Ingeniería Química y Textil, Universidad Nacional de Ingeniería

La ecuación 8.2 representa una familia de curvas en el plano x − y , obtenida cada

una de ellas para un valor particular de c , como se muestra en la figura 8.1.

F3 = 0

y0
F2 = 0, con y ( x0 ) = y0

F1 = 0

x
x0
Figura 8.1. Representación gráfica de la solución general, ecuación (8.2).

Cada una de estas curvas corresponde a una solución particular de la EDO,

ecuación 8.1, y analíticamente dichas constantes se obtienen exigiendo que la
solución de esa ecuación pase por algún punto ( x0 , y 0 ) ; esto es, que

y (x0 ) = y 0 (8.3)

En la práctica la gran mayoría de las ecuaciones no pueden resolverse utilizando

las técnicas analíticas, y se deberá recurrir a los métodos numéricos.

Cuando se usan métodos numéricos no se encuentran soluciones de la forma

F ( x, y, c ) = 0 , ya que trabajan con números y dan por resultado números. Sin
embargo, el propósito usual de encontrar una solución es determinar valores de y
(números) correspondientes a valores específicos de x , lo cual es factible con los
mencionados métodos numéricos sin tener que encontrar F ( x, y, c ) = 0 .

El problema de valor inicial (PVI) por resolver numéricamente queda formulado

como sigue

2
Alex W.Pilco Nuñez, Facultad de Ingeniería Química y Textil, Universidad Nacional de Ingeniería

a) Una ecuación diferencial de primer orden, ecuación 8.1

b) El valor de y en un punto conocido x 0 (condición inicial)
c) El valor de x f donde se quiere conocer el valor de y o ( y f )
que en lenguaje matemático quedará así


= f ( x, y )
dy
 dx

PVI  y (x0 ) = y 0
y (x f ) = ?
 (8.4)



Formulado el problema de valor inicial, a continuación se describe las técnicas

numéricas para resolverlo.

Método de Euler

El método de Euler es el método más simple de los métodos numéricos para

resolver un problema de valor inicial, ecuación 8.4. Consiste en dividir el intervalo
que va de x 0 a x f en n subintervalos de ancho h (tamaño de paso); o sea

x f − x0
h= (8.5)
n

de manera que se obtiene un conjunto discreto de n puntos: x 0 , x1 , x 2 , . . . , x n

 
( x f se convierte en x n ) del intervalo de interés x0 , x f . Para cualquiera de estos
puntos se cumple que

xi = x0 + i h i = 0, 1, 2, . . . , n − 1 (8.6)

La condición inicial y ( x0 ) = y 0 representa el P0 = (x0 , y 0 ) por donde pasa la curva

solución de la ecuación 8.4, la cual por simplicidad se denotará como F (x ) = y en
lugar de F ( x, y, c1 ) = 0 .

Con el punto P0 se puede evaluar la primera derivada de F (x) en ese punto

dy
F ´( x ) = = f ( x0 , y0 ) (8.7)
dx P0

3
Alex W.Pilco Nuñez, Facultad de Ingeniería Química y Textil, Universidad Nacional de Ingeniería

Con esta información se traza una recta, aquella que pasa por P0 y de pendiente
f (x0 , y 0 ) . Esta recta aproxima F (x) en una vecindad de x 0 . Tómese la recta como
reemplazo de F (x) localice en ella (la recta) el valor de y correspondiente a x1 .

F (x f )
y

F (x1 ) f (x0 , y 0 )

y1

y0
h

x
x0 x1 x f = xn

Figura 8.2. Deducción gráfica del método de Euler.

Entonces de la figura 8.2

y1 − y0
= f ( x0 , y 0 ) (8.8)
x1 − x0

Se resuelve para y1

y1 = y 0 + (x1 − x0 ) f (x0 , y 0 ) = y 0 + h f (x0 , y 0 ) (8.9)

Es evidente que la ordenada y1 calculada de esta manera no es igual a F (x1 ) ,

pues existe un pequeño error. No obstante, el valor de y1 sirve para aproximar
F´ (x) en el punto P = (x1 , y1 ) y repetir el procedimiento anterior a fin de generar la
sucesión de aproximaciones siguiente

4
Alex W.Pilco Nuñez, Facultad de Ingeniería Química y Textil, Universidad Nacional de Ingeniería

x0 y0
x1 = x0 + h y1 = y 0 + h f ( x0 , y 0 )
x 2 = x1 + h y 2 = y1 + h f ( x1 , y1 )
.
. (8.10)
.
xi +1 = xi + h y i +1 = y i + h f (xi , y i )
.
.
.
x n = x n −1 + h y n = y n −1 + h f ( x n −1 , y n −1 )

Como se muestra en la figura 8.3, en esencia se trata de aproximar la curva

y = F (x ) por medio de una serie de segmentos de línea recta.

Como la aproximación a una curva mediante una línea recta no es exacta, se

comete un error propio del método mismo. Dicho error puede disminuirse tanto
como se quiera (al menos teóricamente) reduciendo el valor del tamaño de paso
h , pero a cambio de un mayor número de cálculos y tiempo de máquina.

y = F (x )
y
Error final

y Euler

x0 x1 x3 x4 x f = xn

Figura 8.3. Aplicación repetida del método de Euler.

5
Alex W.Pilco Nuñez, Facultad de Ingeniería Química y Textil, Universidad Nacional de Ingeniería

Ejemplo 8.1: Resuelva el siguiente

 dy
 = x− y
dx

PVI  y( 0 ) = 2
 y( 1 ) = ?



mediante el método de Euler y divida en cinco subintervalos.

Solución

Como el intervalo de interés es 0,1, el valor del tamaño de paso h es

1− 0
h= = 0.2
5

con lo cual se generan los argumentos x y los valores de y , mediante el método

de Euler

x0 = 0 y0 = 2
x1 = x0 + h = 0 + 0.2 = 0.2 y1 = y0 + h f ( x0 , y0 ) = 2 + 0.2  0 − 2 = 1.6
x2 = x1 + h = 0.2 + 0.2 = 0.4 y2 = y1 + h f ( x1 , y1 ) = 1.6 + 0.2  0.2 − 1.6 = 1.32
x3 = x2 + h = 0.4 + 0.2 = 0.6 y3 = y2 + h f ( x2 , y2 ) = 1.32 + 0.2  0.4 − 1.32 = 1.136
x4 = x3 + h = 0.6 + 0.2 = 0.8 y4 = y3 + h f ( x3 , y3 ) = 1.136 + 0.2  0.6 − 1.136 = 1.0288
x5 = x4 + h = 0.8 + 0.2 = 1.0 y5 = y4 + h f ( x4 , y4 ) = 1.0288 + 0.2  0.8 − 1.0288 = 0.98304

i xi yi
0 0.0 2.00000
1 0.2 1.60000
2 0.4 1.32000
3 0.6 1.13600
4 0.8 1.02880
5 1.0 0.98304

El valor de y5 = y(x5 ) = y(1.0) = 0.98304 obtenido por este método puede

compararse con la solución analítica que es 1.10364; el error relativo porcentual
es 11%.

6
Alex W.Pilco Nuñez, Facultad de Ingeniería Química y Textil, Universidad Nacional de Ingeniería

Ejemplo 8.2: Resuelva el siguiente

 dy
 =x
dx

PVI  y (2 ) = 4
 y( 4 ) = ?



mediante el método de Euler. Hacer cuatro subintervalos. Además resuelva

analíticamente para encontrar el error relativo porcentual.

Solución

Como el intervalo de interés es 2, 4 , el valor del tamaño de paso h es

4−2
h= = 0.5
4

con lo cual se generan los argumentos x y los valores de y , mediante el método

de Euler
x0 = 2 y0 = 4
x1 = 2.5 y1 = y 0 + h f (x0 , y 0 ) = 4 + 0.5 2 = 5
x2 = 3 y 2 = y1 + h f (x1 , y1 ) = 5 + 0.5 2.5 = 6.25
x3 = 3.5 y3 = y 2 + h f (x 2 , y 2 ) = 6.25 + 0.5 3 = 7.75
x4 = 4 y 4 = y 3 + h f (x3 , y 3 ) = 7.75 + 0.5 3.5 = 9.5

El valor obtenido por este método es y4 = y(x4 ) = y(4) = 9.5 . Resolviendo

analíticamente y considerando la condición inicial, se tiene
x2
y= +2
2

A partir de la solución analítica se puede presentar una tabla con los siguientes
resultados

yi , aprox. − y i , analítico
i xi yi , aproximado yi , analítico ERP = x100
y i , analítico
0 2.000 4.000 4.000 0.000
1 2.500 5.000 5.125 2.439
2 3.000 6.250 6.500 3.846
3 3.500 7.750 8.125 4.615
4 4.000 9.500 10.000 5.000

7
Alex W.Pilco Nuñez, Facultad de Ingeniería Química y Textil, Universidad Nacional de Ingeniería

Ejemplo 8.3: Un tanque cilíndrico de fondo plano con un diámetro de 1.5 m ,

contiene un líquido hasta una altura de 3 m respecto del fondo plano. Se desea
saber la altura del líquido en metros dentro del tanque 90 segundos después de
que se abre completamente la válvula de salida, la cual da una velocidad de
0.6 2 g a m / s , el diámetro de tubo de salida es 0.1 m y la aceleración de la
gravedad es 9.81 m / s 2 . Resolver con el método de Euler y considerar sobre el
intervalo de la variable independiente tres subintervalos.

Solución

El vaciado de un tanque cilíndrico se modela haciendo un balance de masa con la

siguiente expresión
dm • •
= m entrada − m salida
dt
d (V  )
= 0−v A
dt
•
Donde m es el masa del líquido en el tanque en un tiempo t , y m es el flujo
másico del líquido.

Se considera que el descenso del nivel del líquido se realiza a temperatura

constante

dV
= −v A
dt

8
Alex W.Pilco Nuñez, Facultad de Ingeniería Química y Textil, Universidad Nacional de Ingeniería

 D2 d a  d2 
= −0.6 2 g a  
4 dt  4 
2
da d
= −0.6 2 g   a
dt D

Reemplazando los datos, queda

2
da d
= −0.6 2 g   a
dt D
da
= −0.01181 a
dt

Al considerar como tiempo cero el momento de abrir la válvula y además la altura

buscada a un tiempo de 180 segundos, se llega a

 da
 = −0.01181 a
dt

PVI  a( 0 ) = 3

 a ( 90 ) = ?


Se usa el método de Euler para resolver este PVI . Sobre el intervalo de interés
 0, 90  , el valor del tamaño de paso h es
90 − 0
h= = 30
3

con lo cual se generan los argumentos x y los valores de y mediante el método

de Euler
xi +1 = xi + h yi +1 = yi + h f ( xi , yi ) (8.10)

En la siguiente tabla se muestran los resultados de los cálculos para la altura del
líquido a para un tiempo dado t mediante el método de Euler

i ti ai
0 0 3.00
1 30
2 60
3 90

9
Alex W.Pilco Nuñez, Facultad de Ingeniería Química y Textil, Universidad Nacional de Ingeniería

Ejemplo 8.4: Un tanque con agitación contiene 20 m3 de una solución acuosa de

sal A cuya concentración es 0.5 kg / m3 . Si en un instante determinado se deja
ingresar un flujo volumétrico de 3 m3 / h , con una concentración de la sal de
2 kg / m3 y al mismo tiempo se hace salir un flujo volumétrico de 2 m3 / h . Asuma en
todo instante la densidad de la solución acuosa de sal A de los flujos que entra y
sale, como también del tanque, iguales y constantes.

a) Calcular la concentración de la sal en kg / m3 , en el instante en que inicia el

rebalse, si la capacidad máxima del tanque con agitación es 40 m3 . Usar el
método de Euler. Considerar sobre el intervalo de la variable independiente
cuatro subintervalos.
b) Resolver analíticamente y calcular el error relativo porcentual.

Solución

a) Aplicación de la ecuación de balance de masa a la sal

d ( mA ) • •
= V entrada C A0 − V salida C A
dt
d (VC A )
= ( 3)( 2 ) − 2C A
dt

Donde V es el volumen de la solución en un tiempo t

10
Alex W.Pilco Nuñez, Facultad de Ingeniería Química y Textil, Universidad Nacional de Ingeniería

dCA dV
V + CA = 6 − 2CA
dt dt

Aplicación de la ecuación de balance de masa a todo el sistema

d ( V )
= ( 3)(  ) − ( 2 )(  )
dt
dV
=1
dt


V
dV =  dt
t

20 0

V = 20 + t

Reemplazando en la ecuación de balance de masa para la sal

dC A
( 20 + t ) + C A (1) = 6 − 2C A
dt
dC A 6 − 3C A
=
dt 20 + t

Se inicia el rebalse cuando el volumen de la solución acuosa de sal alcanza los

40 m3
V = 20 + t
40 = 20 + t
t = 20 horas

Entonces nuestro problema de valor inicial será

dC A 6 − 3C A
=
dt 20 + t
C A ( 0 ) = 0.5
C A ( 20 ) = ?

Sobre el intervalo de interés es  0, 20 , el valor del tamaño de paso h es

20 − 0
h= =5
4

con lo cual se generan los argumentos x y los valores de y , mediante el método

de Euler

11
Alex W.Pilco Nuñez, Facultad de Ingeniería Química y Textil, Universidad Nacional de Ingeniería

xi +1 = xi + h yi +1 = yi + h f ( xi , yi ) (8.10)

En la siguiente tabla se muestran los resultados de los cálculos para C A mediante

el método de Euler

i ti C A ,i
0 0 0.5000
1 5 1.6250
2 10 1.8500
3 15 1.9250
4 20 1.9571

b) Resolviendo analíticamente

De la ecuación de balance de masa para la sal

dC A
6 − 2C A = ( 20 + t ) + C A (1)
dt
dC A 6 − 3C A
=
dt 20 + t
dC dt
0.5 6 − 3CA A = 0 20 + t
CA 20

1
− ln ( C A − 2 ) C0.5A = ln ( 20 + t ) 20
0
3
C A = 1.8125 kg / m3

El valor de C A,4 ( 20 ) = 1.9571 obtenido por el método de Euler puede compararse

con la solución analítica que es C A = 1.8125 ; el error relativo porcentual es 8%.

Método de Taylor

Antes de explicar este método, conviene hacer una acotación al método de Euler.

Puede decirse que el método de Euler utiliza los primeros dos términos de la serie
de Taylor para su primera iteración; o sea

F ( x1 )  y1 = F ( x0 ) + F ' ( x0 )( x1 − x0 ) (8.11)

donde se señala que y1 no es igual a F ( x1 ) .

12
Alex W.Pilco Nuñez, Facultad de Ingeniería Química y Textil, Universidad Nacional de Ingeniería

Esto pudo hacer pensar que para encontrar y2 , se expandió de nuevo F ( x ) en

serie de Taylor, como sigue

F ( x2 )  y2 = F ( x1 ) + F ' ( x1 )( x2 − x1 ) (8.12)

sin embargo, no se dispone de los valores exactos de F ( x1 ) y F ' ( x1 ) y,

rigurosamente hablando, son los que deben usarse en una expansión de Taylor de
F ( x ) - en este caso alrededor de x1 - ; por tanto, el lado derecho de la ecuación
8.12 no es evaluable. Por ello, sólo en la primera iteración, para encontrar y1 , se
usa realmente una expansión en serie de Taylor de F ( x ) , aceptando desde luego
que se tienen valores exactos en la condición inicial y0 = F ( x0 ) . Después de eso,
se emplea la ecuación

yi +1 = yi + f ( xi , yi )( xi +1 − xi )
(8.13)
= F ( xi ) + F ' ( xi )( xi +1 − xi )

que guarda similitud con una expansión en serie de Taylor.

Aclarando este punto, a continuación se aplicará la información acerca de las

series de Taylor para mejorar la exactitud del método de Euler y obtener
extensiones que constituyen la familia de métodos llamados algoritmos de Taylor.

Si se usan tres términos en lugar de dos en la expansión de F ( x1 ) , entonces

( x1 − x0 )
2

F ( x1 )  y1 = F ( x0 ) + F ' ( x0 )( x1 − x0 ) + F '' ( x0 ) (8.14)

2!
Como
dF ' ( x ) df ( x, y )
F '' ( x ) = =
dx dx

y
h = x1 − x0

La primera iteración de la ecuación 8.14 tomaría la forma

h2 df ( x, y )
y1 = y0 + h f ( x0 , y0 ) + x0 , y0 (8.15)
2! dx

13
Alex W.Pilco Nuñez, Facultad de Ingeniería Química y Textil, Universidad Nacional de Ingeniería

Ahora cabe pensar que usando una fórmula de iteración basada en la ecuación
8.15 para obtener y2 , y3 , . . . yn mejoraría la exactitud obtenida con la ecuación
8.11. Se propone entonces la fórmula

h2 df ( x, y )
yi +1 = yi + h f ( xi , yi ) + xi , yi (8.16)
2! dx

La utilidad de esta ecuación depende de cuán fácil sea la diferenciación de

f ( x, y ) . Si f ( x, y ) es una función sólo de x , la diferenciación con respecto a x
es relativamente fácil y la fórmula propuesta es muy práctica.

Si, como es el caso general, f ( x, y ) es una función de x y y , habrá que usar

derivadas totales. La derivada total de f ( x, y ) con respecto a x está dada por

df ( x, y ) f ( x, y ) f ( x, y ) dy
= + (8.17)
dx x y dx

Si se aplican las ideas vistas en el método de Euler pero empleando como fórmula
la ecuación 8.16, se obtiene el método de Taylor de segundo de orden. Este último
es indicativo de la derivada de mayor orden que se emplea y de cierta exactitud.
Con esta terminología, al método de Euler le correspondería el nombre de método
de Taylor de primer orden.

Ejemplo 8.5: Resuelva el siguiente

 dy
 = x− y
dx

PVI  y( 0 ) = 2
 y( 1 ) = ?



mediante el método de Taylor de segundo orden y divida en cinco subintervalos.

Solución

Como el intervalo de interés es 0,1, el valor del tamaño de paso h es

1− 0
h= = 0.2
5
Se aplica la ecuación 8.16

14
Alex W.Pilco Nuñez, Facultad de Ingeniería Química y Textil, Universidad Nacional de Ingeniería

h2 df ( x, y )
yi +1 = yi + h f ( xi , yi ) + xi , yi (8.16)
2! dx

se evalúa primero la ecuación 8.17

df ( x, y ) f ( x, y ) f ( x, y ) dy
= +
dx x y dx
df ( x, y )
= 1 + ( −1)( x − y ) = 1 − x + y
dx

con lo cual se generan los valores de y

x0 = 0 y0 = 2
h 2 df ( x, y )
x1 = x0 + h = 0 + 0.2 = 0.2 y1 = y0 + h f ( x0 , y0 ) + x0 , y0
2! dx
0.22
= 2 + 0.2  0 − 2 + (1 − 0 + 2 ) = 1.66
2
h 2 df ( x, y )
x2 = x1 + h = 0.2 + 0.2 = 0.4 y2 = y1 + h f ( x1 , y1 ) + x1 , y1
2! dx
0.22
= 1.66 + 0.2  0.2 − 1.66 + (1 − 0.2 + 1.66 ) = 1.4172
2
h 2 df ( x, y )
x3 = x2 + h = 0.4 + 0.2 = 0.6 y3 = y2 + h f ( x2 , y2 ) + x2 , y2
2! dx
= 1.254104
h 2 df ( x, y )
x4 = x3 + h = 0.6 + 0.2 = 0.8 y4 = y3 + h f ( x3 , y3 ) + x3 , y3
2! dx
= 1.15636528
h 2 df ( x, y )
x5 = x4 + h = 0.8 + 0.2 = 1.0 y5 = y4 + h f ( x4 , y4 ) + x4 , y4
2! dx
= 1.11221953

i xi yi
0 0.0 2.00000000
1 0.2 1.66000000
2 0.4 1.41720000
3 0.6 1.25410400
4 0.8 1.15636528
5 1.0 1.11221953

15
Alex W.Pilco Nuñez, Facultad de Ingeniería Química y Textil, Universidad Nacional de Ingeniería

El valor de y5 = y ( x5 ) = y (1.0 ) = 1.11222 obtenido por este método puede

compararse con la solución analítica que es 1.10364; el error relativo porcentual
es 0.8%.

Ejemplo 8.6: Un tanque con agitación contiene 20 m3 de una solución acuosa de

sal A cuya concentración es 0.5 kg / m3 . Si en un instante determinado se deja
ingresar un flujo volumétrico de 3 m3 / h , con una concentración de la sal de
2 kg / m3 y al mismo tiempo se hace salir un flujo volumétrico de 2 m3 / h . Asuma en
todo instante la densidad de la solución acuosa de sal A de los flujos que entra y
sale, como también del tanque, iguales y constantes.

a) Calcular la concentración de la sal en kg / m3 , en el instante en que inicia el

rebalse, si la capacidad máxima del tanque con agitación es 40 m3 . Usar el
método de Taylor de segundo orden. Considerar sobre el intervalo de la
variable independiente cuatro subintervalos.
b) Resolver analíticamente y calcular el error relativo porcentual.

Solución

a) Aplicación de la ecuación de balance de masa a la sal

d (VCA )
( 3)( 2) − 2CA =
dt

16
Alex W.Pilco Nuñez, Facultad de Ingeniería Química y Textil, Universidad Nacional de Ingeniería

Donde V es el volumen de la solución en un tiempo t

dCA dV
6 − 2CA = V + CA
dt dt

Aplicación de la ecuación de balance de masa a todo el sistema

d ( V )
( 3)(  ) − ( 2 )(  ) =
dt
dV
=1
dt


V
dV =  dt
t

20 0

V = 20 + t

Reemplazando en la ecuación de balance de masa para la sal

dC A
6 − 2C A = ( 20 + t ) + C A (1)
dt
dC A 6 − 3C A
=
dt 20 + t

Se inicia el rebalse cuando el volumen de la solución acuosa de sal alcanza los

40 m3
V = 20 + t
40 = 20 + t
t = 20 horas

Entonces nuestro problema de valor inicial será

dC A 6 − 3C A
=
dt 20 + t
C A ( 0 ) = 0.5
C A ( 20 ) = ?

Por el método de Taylor de segundo orden, sobre el intervalo de interés es  0, 20 ,

el valor del tamaño de paso h es
20 − 0
h= =5
4

Se aplica la ecuación 8.16

17
Alex W.Pilco Nuñez, Facultad de Ingeniería Química y Textil, Universidad Nacional de Ingeniería

h2 df ( x, y )
yi +1 = yi + h f ( xi , yi ) + xi , yi (8.16)
2! dx

se evalúa primero la ecuación 8.17

df ( x, y ) f ( x, y ) f ( x, y ) dy
= +
dx x y dx
df ( x, y ) ( 6 − 3C A ) − 3 ( 6 − 3C A ) = − 4 ( 6 − 3CA )
=−
( t + 20 ) ( t + 20 ) ( t + 20 ) ( t + 20 )
2 2
dx

En la siguiente tabla se muestran los resultados de los cálculos para C A

i ti C A ,i
0 0 0.5000
1 5 1.0625
2 10 1.4000
3 15 1.5000
4 20 1.7224

b) Resolviendo analíticamente

De la ecuación de balance de masa para la sal

dC A
6 − 2C A = ( 20 + t ) + C A (1)
dt
dC A 6 − 3C A
=
dt 20 + t
dC dt
0.5 6 − 3CA A = 0 20 + t
CA 20

1
− ln ( C A − 2 ) C0.5A = ln ( 20 + t ) 20
0
3
C A = 1.8125 kg / m3

El valor de C A,4 ( 20 ) = 1.7224 obtenido por método de Taylor de segundo orden

puede compararse con la solución analítica que es C A = 1.8125 ; el error relativo
porcentual es 5%.

18
Alex W.Pilco Nuñez, Facultad de Ingeniería Química y Textil, Universidad Nacional de Ingeniería

Método de Runge – Kutta de cuarto orden

Este método da resultados sorprendentemente casi exactos sin necesidad de

utilizar valores extremadamente pequeños de h . No daremos justificación alguna
para el método sino meramente listaremos las diversas fórmulas que intervienen y
explicaremos cómo se utilizan.

Vamos a considerar ahora el método llamado método de Runge – Kutta de cuarto

orden para aproximar los valores de la solución del problema de valor inicial


= f ( x, y )
dy
 dx

PVI  y (x0 ) = y 0
y (x f ) = ?
 (8.4)



Formulado el problema de valor inicial, a continuación se describe el método de

Runge – Kutta de cuarto orden

x0 y0
x1 = x0 + h y1 = y0 + h  0
x2 = x1 + h y2 = y1 + h 1
.
.
.
xi +1 = xi + h yi +1 = yi + h  i
.
. (8.18)
.
xn = xn −1 + h yn = yn −1 + h  n −1

donde
1
i = ( k1 + 2k2 + 2k3 + k4 )
6
y además
k1 = f ( xi , yi )
k2 = f ( xi + h / 2, yi + h k1 / 2 )
k3 = f ( xi + h / 2, yi + h k2 / 2 )
k4 = f ( xi + h, yi + h k3 )
19
Alex W.Pilco Nuñez, Facultad de Ingeniería Química y Textil, Universidad Nacional de Ingeniería

Ejemplo 8.7: Se hacen reaccionar isotérmicamente 260 g de acetato de etilo

( CH 3COOC2 H 5 ) con 175 g de hidróxido de sodio ( NaOH ) en solución acuosa
(ajustando el volumen total de la mezcla reaccionante a 5 litros de solución) para
producir acetato de sodio ( CH 3COONa ) y etanol ( C2 H 5OH ), de acuerdo con la
ecuación estequiométrica

CH3COOC2 H5 + NaOH ⎯⎯
k
→ CH3COONa + C2 H5OH

L
Si la constante de velocidad de reacción k está dada por k = 1.44 x10−2
mol . min

Determine la cantidad de acetato de sodio y etanol en moles, presentes 30

minutos después de iniciada la reacción. Use el método de Runge-Kutta de cuarto
orden y considere como tamaño de paso con respecto al tiempo 3 min.

Solución

Denotamos para las especies químicas de la reacción

C1 : concentración del acetato de etilo

C2 : concentración del hidróxido de sodio
C3 : concentración del acetato de sodio
C4 : concentración del etanol

Se considera que la velocidad de reacción respecto al acetato de etilo está dada

por la siguiente expresión
dC
r = − 1 = k C1 C2
dt

Como la mezcla reaccionante se encuentra en fase líquida, se puede asumir que

el volumen de la mezcla se mantiene constante durante la reacción, por tanto se
tiene
1 dn1 k
− = 2 n1 n2
V dt V
dn k
− 1 = n1 n2
dt V

donde n1 y n2 son moles de las especies químicas reactantes en un tiempo

determinado. Si x es la cantidad de moles que reaccionó el acetato de etilo en un
tiempo determinado, entonces se tiene

20
Alex W.Pilco Nuñez, Facultad de Ingeniería Química y Textil, Universidad Nacional de Ingeniería

d ( n1,0 − x )
−
dt
=
k
V
( n1,0 − x )( n2,0 − x )
donde n1,0 y n2,0 son moles iniciales de la especies químicas reactantes, por
consiguiente, la última ecuación se puede expresar de la siguiente forma

= ( n1,0 − x )( n2,0 − x )
dx k
dt V
dx 1.44 x10−2  260  175 
=  − x  − x
dt 5  88  40 
dx
= 2.88 x10−3 ( 2.954 − x )( 4.375 − x )
dt

Entonces nuestro problema de valor inicial será

dx
= 2.88 x10−3 ( 2.954 − x )( 4.375 − x )
dt
x ( 0) = 0
x ( 30 ) = ?

Por el método de Runge-Kutta de cuarto orden, sobre el intervalo de interés es

0,30 , el valor del tamaño de paso h es
30 − 0
h= =3
10

En la siguiente tabla se muestran los resultados de los cálculos para las

concentraciones del soluto en los tanques 1, 2 y 3

Iteración t , min x, mol

0 0 0
1 3
2 6
3 9
4 12
5 15
6 18
7 21
8 24
9 27
10 30

21
Alex W.Pilco Nuñez, Facultad de Ingeniería Química y Textil, Universidad Nacional de Ingeniería

Ecuaciones diferenciales ordinarias de orden superior y sistemas de

ecuaciones diferenciales ordinarias

Cuando en el problema de valor inicial aparecen una ecuación diferencial de orden

n , n condiciones especificadas en un punto x 0 y un punto x f donde hay que
encontrar el valor de y (x f ) , se tiene el problema de valor inicial general (PVIG):



dny
dx n
(
= f x, y, y ' , y ' ' , . . . , y ( n −1) )



y (x0 ) = y 0 , y ' (x0 ) = y 0 ' , y ' ' ( x0 ) = y 0 ' ' , . . . , y ( n −1) (x0 ) = y 0
( n −1)
PVIG  (8.19)


 y (x f ) = ?



Para resolver la ecuación anterior no se desarrollan nuevos métodos, sino que se

emplea una extensión de los métodos estudiados en este capítulo. Para ello
necesitaremos primero pasar la ecuación diferencial ordinaria o EDO de la
ecuación 8.19 a un sistema de ecuaciones diferenciales simultáneas de primer
orden cada una.

Ejemplo 8.8: Una de las ecuaciones diferenciales ordinarias más empleadas en la

matemática física es la ecuación de Bessel

x 2 y' ' + x y' + ( x 2 − n 2 ) y = 0

donde n puede tener cualquier valor, pero generalmente toma un valor entero.
Escriba esta ecuación como un sistema de ecuaciones diferenciales ordinarias de
primer orden.

Solución

La ecuación se pone en la forma normal

1  n2 
y ' ' = − y ' +  2 − 1 y
x x 

Algunas veces es más conveniente para los cálculos computacionales emplear

y =y
y' = z

22
Alex W.Pilco Nuñez, Facultad de Ingeniería Química y Textil, Universidad Nacional de Ingeniería

como nuevas variables. Se deriva la segunda y se tiene:

y '' = z '
El sistema queda
y' = z
1  n2 
z' = − z +  2 − 1 y
x x 

sistema que solo podrá resolverse para valores de x distintos de cero.

Ejemplo 8.9: Pase la ecuación diferencial ordinaria

d 2 y dy
2
+ = x2 + y 2
dx dx

a un sistema de dos ecuaciones diferenciales ordinarias simultáneas de primer

orden.

Solución

Con el despeje de la derivada de segundo orden se tiene

d2y dy
2
=− + x2 + y 2
dx dx

El cambio de variable para los cálculos computacionales son

y =y
dy
=z
dx
Se deriva la segunda
d2y dz
2
=
dx dx

y las nuevas variables se sustituyen en la ecuación diferencial, con lo cual resulta

dy
=z
dx
dz
= − z + x2 + y 2
dx

23
Alex W.Pilco Nuñez, Facultad de Ingeniería Química y Textil, Universidad Nacional de Ingeniería

Si ha a las dos ecuaciones simultáneas que forman un sistema, se aplica el

método de Runge – Kutta de cuarto orden, se tiene

y i +1 = y i +
h
(k1 + 2k 2 + 2k 3 + k 4 )
6
= z i + (c1 + 2 c 2 + 2 c3 + c 4 )
h
z i +1
6

las cuales se calculan alternadamente, y las k y c se obtienen de

k1 = f 1 ( x i , y i , z i )
c1 = f 2 (xi , y i , z i )
k 2 = f1 (xi + h / 2, y i + h k1 / 2, z i + h c1 / 2)
c 2 = f 2 (xi + h / 2, y i + h k1 / 2, z i + h c1 / 2)
k 3 = f1 (xi + h / 2, y i + h k 2 / 2, z i + h c 2 / 2)
c3 = f 2 (xi + h / 2, y i + h k 2 / 2, z i + h c 2 / 2)
k 4 = f1 (xi + h, y i + h k 3 , z i + h c3 )
c 4 = f 2 (xi + h, y i + h k 3 , z i + h c3 )

calculadas en ese orden.

Ejemplo 8.10: Resuelva el siguiente problema de valor inicial general por el

método de Runge – Kutta de cuarto orden. (Use un programa)

 1  1 
 y '' = − y' +  2 − 1 y

x x 
 y ( 1 ) =1
PVIG 
 y '( 1 ) = 2

 y (3)=?

Nótese que la EDO es la ecuación de Bessel con n = 1 (véase ejemplo 8.5). Divida
el intervalo de interés 1, 3 en ocho subintervalos.

Solución

Con el tamaño de paso de h = 0.25 , se tiene el nuevo PVIG

24
Alex W.Pilco Nuñez, Facultad de Ingeniería Química y Textil, Universidad Nacional de Ingeniería

d2 y 1 dy  1 
 dx 2 = − x dx +  x 2 − 1 y
  
 y ( x0 =1 ) = 1
PVIG 
 dy =2
 dx ( x0 =1)

 y ( 3 ) = ?

El cambio de variable para los cálculos computacionales son

y =y
dy
=z
dx
Se deriva la segunda
d2y dz
2
=
dx dx

y las nuevas variables se sustituyen en la ecuación diferencial, con lo cual resulta

 dy
 dx = z

 dz = − 1 z +  1 − 1 y
 dx  2 
x x 

PVIG  y (1) = 1

 z (1) = 2

 y ( 3) = ?
 z ( 3) = ?


y (1.00 ) = 1 z (1.00 ) = 2
y (1.25 ) = 1.441151281 z (1.25 ) = 1.539492187
y (1.50 ) = 1.771863249 z (1.50 ) = 1.107293533
y (1.75 ) = 1.994766280 z (1.75 ) = 0.675599895
y (2.00 ) = 2.109754328 z (2.00 ) = 0.245291635
y (2.25 ) = 2.118486566 z (2.25 ) = − 0.172076357
y (2.50 ) = 2.026089844 z (2.50 ) = − 0.561053191
y (2.75 ) = 1.841680320 z (2.75 ) = − 0.905578495
y (3.00 ) = 1.578253875 z (3.00 ) = − 1.190934201

25
Alex W.Pilco Nuñez, Facultad de Ingeniería Química y Textil, Universidad Nacional de Ingeniería

Ejemplo 8.11: En el diagrama de la figura se muestra una corriente de 300 L/min

y concentración de soluto de 50 g/L que ingresa al tanque 1. Del tanque 1 sale una
corriente de 300 L/min e ingresa al tanque 2. Del tanque 2 sale una corriente de
450 L/min e ingresa al tanque 3. Del tanque 3 sale una corriente de 450 L/min,
parte de ella, se toma una corriente de recirculación de 150 L/min y se lleva al
tanque 2 y el resto sale del sistema.

El volumen se conserva constante en cada tanque e igual a 1000 litros, determine

la concentración del soluto (g/L) en cada tanque 10 minutos después de iniciado el
proceso. Use el método de Runge-Kutta de cuarto orden y considere como
tamaño de paso con respecto al tiempo 2 min.

Tanque 1 Tanque 2 Tanque 3

C1 (0) =30 g/L C2 (0)=30 g/L C3(0)=30 g/L
V1=1000 L V2=1000 L V3=1000 L

Solución

Aplicación de la ecuación de balance de masa al soluto en cada tanque

d (V1 C1 )
Tanque 1 = ( 300 )( 50 ) − 300 C1
dt
d (V2 C2 )
Tanque 2 = 300 C1 + 150 C3 − 450 C2
dt
d (V3 C3 )
Tanque 3 = 450 C2 − 450 C3
dt

Donde V1 , V2 y V3 , son los volúmenes de las soluciones en los tanques 1, 2 y 3

que de acuerdo al problema son constantes en el proceso

26
Alex W.Pilco Nuñez, Facultad de Ingeniería Química y Textil, Universidad Nacional de Ingeniería

dC1
Tanque 1 = 15 − 0.3 C1
dt
dC2
Tanque 2 = 0.3 C1 + 0.15 C3 − 0.45 C2
dt
dC3
Tanque 3 = 0.45 C2 − 0.45 C3
dt

El sistema de ecuaciones diferenciales ordinarias es

dC1
= 15 − 0.3 C1
dt
dC2
= 0.3 C1 + 0.15 C3 − 0.45 C2
dt
dC3
= 0.45 C2 − 0.45 C3
dt
C1 ( 0 ) = 30
C2 ( 0 ) = 30
C3 ( 0 ) = 30
C1 (10 ) = ?
C2 (10 ) = ?
C3 (10 ) = ?

Por el método de Runge-Kutta de cuarto orden, sobre el intervalo de interés es

 0,10 , el valor del tamaño de paso es h = 2 es
En la siguiente tabla se muestran los resultados de los cálculos para las
concentraciones del soluto en los tanques 1, 2 y 3

Iteración t , min C1 , g / L C2 , g / L C3 , g / L
0 0 30 30 30
1 2
2 4
3 6
4 8
5 10

27