INTRODUCCION

Las leyes que gobiernan los fenómenos de la naturaleza

se expresan habitualmente en forma de ecuaciones
diferenciales. Las ecuaciones del movimiento de los
cuerpos (la segunda ley de Newton) es una ecuación
diferencial de segundo orden, como lo es la ecuación que
describe los sistemas oscilantes, la propagación de las
ondas, la transmisión del calor, la difusión, el movimiento
de partículas subatómicas, etc.
Pocas ecuaciones diferenciales tienen una solución
analítica sencilla, la mayor parte de las veces es necesario
realizar aproximaciones, estudiar el comportamiento del
sistema bajo ciertas condiciones. Así, en un sistema tan
simple como un péndulo, la amplitud de la oscilación ha
de ser pequeña y el rozamiento ha de ser despreciable,
para obtener una solución sencilla que describa
aproximadamente su movimiento periódico.
Para determinar el estado final de un sistema físico, cuyo
comportamiento está descrito por una ecuación
diferencial, se creará un objeto de la clase derivada, y
desde éste llamaremos a la función que describe el
procedimiento numérico, pasándole en uno de sus
parámetros el objeto que describe el estado inicial.
UNIDAD 6: SOLUCION NUMERICA DE ECUACIONES Y
SISTEMAS DE ECUACIONES DIFERENCIALES
Al calcular soluciones numéricas de ecuaciones diferenciales
no se espera hallar una aproximación continua a la solución
del problema de valor inicial, si no aproximaciones en
algunos puntos específicos. Si se requieren valores
intermedios se puede usar interpolación de Lagrange, ajuste
de curvas, etc.

6.1.-Metodos de la serie de Taylor

Una serie de Taylor es una aproximación
de funciones mediante una serie de potencias o suma de
potencias enteras de polinomios como llamados
términos de la serie, dicha suma se calcula a partir de
las derivadas de la función para un determinado valor o
punto suficientemente derivable sobre la función y un
entorno sobre el cual converja la serie. Si esta serie está
centrada sobre el punto cero, , se le denomina serie
de McLaurin.
Esta aproximación tiene tres ventajas importantes:

 la derivación e integración de una de estas series se

puede realizar término a término, que resultan
operaciones triviales;
 se puede utilizar para calcular valores aproximados de
funciones;
 es posible calcular la optimidad de la aproximación.
Algunas funciones no se pueden escribir como serie de
Taylor porque tienen alguna singularidad. En estos casos
normalmente se puede conseguir un desarrollo en serie
utilizando potencias negativas de x (véase Serie de Laurent).

Por ejemplo:
f(x) = exp(−1/x²) se puede desarrollar como serie de
Laurent.
La serie de Taylor de una función f real o compleja
ƒ(x) infinitamente diferenciable en el entorno de un
número real o complejo a es la siguiente serie de potencias:

Que puede ser escrito de una manera más compacta como

la siguiente suma:
Donde:

 n! es el factorial de n
 f (n)(a) denota la n-ésima derivada de f para el valor a de
la variable respecto de la cual se deriva.
La derivada de orden cero de f es definida como la propia f y
tanto (x − a)0 como son ambos definidos como 1 (
= 1). En caso de ser a = 0, como ya se mencionó, la
serie se denomina también de McLaurin.

Cabe destacar que en una serie de Taylor de potencias

centrada en a de la forma

Siempre se puede hacer el cambio de variable (con lo

que en la función a desarrollar original) para

expresarla como

Centrada en 0. Luego hay que deshacer el

cambio de variable. Por ejemplo, si se quiere desarrollar la

función alrededor de a = 1 se puede tomar , de

manera que se desarrollaría centrada en 0.

 6.2 METODO EULER Y EULER MEJORADO
Método de Euler

El método de Euler rara vez se utiliza en la práctica para

obtener la solución aproximada de un problema de valor
inicial, pero se estudia por su simplicidad en la derivación de
la fórmula y de la determinación del error. Los métodos de
orden superior utilizan las mismas técnicas, pero el álgebra
que requieren es mucho más complicada.

Con el método de Euler se obtiene una solución aproximada

de un problema de valor inicial.

Para empezar, se determina la malla {t 0, t1, ... , tN} de paso h,

donde t0 = a y tN = b. En estos puntos es donde se va a
obtener la aproximación de la solución.

Para determinar la fórmula del método, se parte de un

desarrollo de Taylor de la función solución y(t), alrededor de
un punto de la malla, ti, suponiendo que la función y(t)
posee derivadas primera y segunda continuas en (a, b):
Evaluando esta expresión en t = t i+1, para cualquier i, se
tiene:

Pero como ti+1- ti = h, resulta:

Como y(t) satisface la ecuación diferencial, en particular es

y'(ti) = f(ti, yi), entonces reemplazando en la fórmula (4)
resulta:

Si se elimina de la fórmula anterior el término del error, se

puede escribir:

Resultando así la fórmula del método de Euler para

aproximar la solución en un punto de la malla, teniendo una
aproximación en el punto inmediato anterior. Como la
condición en el punto a del problema de valor inicial da el
valor inicial y(t0)= a, se tiene entonces la solución
aproximada en todos los puntos de la malla. Si se llaman y i =
y(ti), se tiene entonces la fórmula de Euler dada en la
fórmula (7):

Ejemplo 1
Consideremos el siguiente problema de valor inicial.

La fórmula de Euler para este problema, tomando N puntos

en el intervalo [1, 2] (sin contar el punto de partida a = 1),
resulta:

Aplicamos el método de Euler para un paso h = 0,2.

Teniendo en cuenta que h = 0,2, la cantidad de puntos en el

intervalo resulta ser N = 5, y entonces la tabla de valores
obtenida con la fórmula dada en (8) resulta:
 t
i y
0 1, 2,0
0 00
0 0
1 1, 2,4
2 00
0 0
2 1, 2,9
4 76
0 0
3 1, 3,8
6 09
0 3
4 1, 5,0
8 28
0 2
5 2, 6,8
0 38
0 4

Ahora, aplicamos la fórmula el método de Euler con N = 20

y N = 50. Representamos gráficamente los puntos
obtenidos, comparándolos con la solución exacta, dada por
la función:
Se ve en los gráficos obtenidos, que a medida que nos
alejamos del valor inicial, la solución aproximada pierde
precisión (se aleja de la solución exacta), para el paso h =
1/20. Cuando se achica el paso, la solución mejora (h =
1/50).

Ejemplo 2.
Use la fórmula de Euler mejorada para aproximar la
solución de y´= 2ty, y(1) = 1, 1 t 1.5, con h = 0.1
Solución :

y*1 = y0 + hf(t0, y0)

y*1 = 1 + 0.1(2(1)(1)) =
1.2
y2 = 1.5479
y3 = 1.9832
y4 = 2.5908
y5 = 3.4509
6.3 MÉTODO DE RUNGE-KUTTA
Uno de los métodos numéricos más utilizado y preciso para
obtener soluciones aproximadas de las ecuaciones
diferenciales es el método de Runge-Kutta de cuarto orden.
Este método consiste en determinar apropiadamente
constantes de manera que una fórmula como:

yi+1 = yi + ak1+ bk2+ ck3+ dk4

Coincida con un desarrollo de Taylor hasta h4, es decir hasta
el quinto termino. Las ki son constantes múltiplos de f(x, y)
evaluadas en ciertos puntos. La deducción del método es
bastante tediosa, por lo que sólo se darán los resultados.

yi+1 = yi + 1/6(k1 + 2k2 +

2k3 + k4)
k1 = h f(ti, yi)
k2 = h f(ti + 1/2h, yi + 1/2k1)
k3 = h f(ti+ 1/2h, yi + 1/2k2)
 k4 = h f(ti + h, yi + k3)

Observe que k2 depende de k1; k3 depende de k2 y así

sucesivamente. Además k2y k3 implican aproximaciones a la
pendiente en el punto medio del intervalo entreti y ti+1=ti+h

Ejemplo
Aproximar la solución de y'= 2ty, y(1)=1, 1 t 1.5, con h =
0.1 usando el método de Runge-Kutta.
Solución:

f(t, y)=2ty, h=0.1, i=0, y0=1, t0=1.0, t1=1.1, t2=1.2, t3=1

1,2,3,4 .3, t4=1.4
y1=y0 + 1/6(k1+2k2+2k3+k4)
k1=h f(ti, yi) = 0.1 (2 (1) (1) ) = 0.2
k2=h f(t0+ 1/2h, y0+ 1/2k1)
k2=(0.1) 2 (1+1/2(0.1) )
(1+1/2(0.2)) = 0.231
k3=h f(t0+1/2h, y0+1/2k2)
k3=(0.1)2(1+1/2(0.1))
(1+1/2(0.231)) = 0.234255
k4=h f(t0+ h, y0+ k3)= 0.2715361

Luego y1=1 + 1/6 (0.2 + 2(0.231) + 2 (0.234255) + 0.2715361)

y1=1.23367435

La siguiente tabla muestra los resultados del método de

Runge-Kutta con h=0.1 redondeados a 4 decimales
 Valor
ti yi
Real
1 1 1
1.1 1.2337 1.2337
1.2 1.5527 1.5527
1.3 1.9937 1.9937
1.4 2.6116 2.6117
1.5 3.4902 3.4904

En la tabla de abajo se comparan los resultados obtenidos a

partir de los métodos de Euler, Taylor de orden 2, Euler
mejorado y Runge-Kutta aplicadas al problemay'=
2ty, y(1)=1, 1 t 1.5 con h = 0.1

Taylor
Euler Runge- Valor
ti Euler de
mejorado Kutta Real
orden 2
1.000
1 1.0000 1.0000 1 1
0
1. 1.200
1.2300 1.2320 1.2337 1.2337
1 0
1. 1.464
1.5427 1.5479 1.5527 1.5527
2 0
1. 1.815
1.9728 1.9832 1.9937 1.9937
3 4
1. 2.287
2.5721 2.5908 2.6116 2.6117
4 4
1. 2.927
3.4188 3.4509 3.4902 3.4904
5 8
 Comparación de métodos numéricos
(Los resultados se redondearon a 4 decimales).

6.4 Solución de sistemas de ecuaciones diferenciales

ordinarias con valores iniciales.
Una ecuación diferencial ordinaria (comúnmente abreviada
"EDO") es la Ecuación Diferencial que relaciona una función
desconocida de una variable independiente con sus
derivadas. Es decir, una sola variable independiente (a
diferencia de las Ecuaciones Diferenciales Parciales que
involucran derivadas parciales de varias variables), y una o
más de sus derivadas respecto de tal variable.

Una ecuación diferencial ordinaria (EDO) puede plantearse,

siendo F una relación o función, como

para representar la EDO en que la función incógnita

(también conocida como variable dependiente), lo es de
una única variable independiente.
En general, una ecuación diferencial lineal de orden n puede
formularse, siendo cada una función dependiente
de t, como:
Una solución de la ecuación (1a) o (1b) será una "familia" de
curvas o funciones del tipo que substituida dentro de
la ecuación la convierte en una igualdad en la que todos los
términos son conocidos.
En la formulación más simple, la función incógnita es una
función para cierto valor real o complejo pero con mayor
generalidad, puede serlo para el valor de un vector o matriz,
lo que lleva a considerar un sistema de ecuaciones
diferenciales ordinarias (EDO) para una única función.
Ecuación diferencial ordinaria

Sea , tal que , la enésima derivada de y,

entonces una ecuación diferencial ordinaria
(EDO)de orden n tiene siguiente forma:

Para funciones vectoriales,

Siendo, tanto ai(x) como r(x) funciones continuas de x. La
función r(x) es llamada el término fuente (traducido del
ingléssource term); si r(x)=0 la ecuación diferencial lineal es
llamada homogénea, de lo contrario es llamada no
homogénea.

Soluciones
Dada una ecuación diferencial
Una solución que no tiene extensión es llamada
una solución general
Una solución general de una ecuación de orden n es una
solución que contiene n variables arbitrarias,
correspondientes a n constantes de integración.
Una solución particular es derivada de la solución general
mediante la fijación de valores particulares para las
constantes, a menudo elegidas para cumplir condiciones
iniciales. Una solución singular es la que no puede derivarse
de la general.
Solución de una EDO de primer orden
Si se considera una ecuación diferencial de la forma:

una ecuación de primer orden resuelta con respecto a la

derivada, se llama su solución general de la anterior
ecuación diferencial, será una función del tipo:

Que depende de una constante arbitraria C. Satisface

ecuación diferencial para cualquier valor de la constante C.
Además cualquiera que sea la condición inicial

Siempre se puede asignar un valor C0 a la constante C, tal

que la función y = φ(x, C0) satisfaga la condición inicial dada.
Se presume que el punto (x0, y0) esta en la región donde se
cumplen las condiciones de existencia y de unicidad de la
solución.

6.5 APLICACIONES
Ejemplo:
Consiste de un resorte o muelle en espiral suspendido de un
soporte rígido con una masa sujeta

Al extremo. Para analizar este fenómeno usamos la ley de

Hooke y la segunda ley de newton.

La ley de Hooke establece que el resorte ejerce una fuerza

de restitución opuesta a la dirección

de alargamiento del resorte peso=kl donde k es la

constante de restitución y l es el

Alargamiento hasta la posición de equilibrio. Para aplicar la

ley de newton, hay q tener en cuenta

q sobre la mas m actúa la fuerza de gravedad (m.g), la

fuerza de restitución (-kx-mg donde xes el alargamiento) , la
fuerza de amortiguación que es proporcional a la velocidad
de la masa )-b(dx/dt) con b la constante de amortiguación )
y las fuerzas externas (F(t)). Así q la ecuacióndiferencial q
resulta es:
M d^2x/dt^2+b dx/dt+kx=f(t)

Aplicaciones a la física:
La Ley de Hooke:
Supongamos que un cuerpo de masa M está sujeto al
extremo de un resorte flexible suspendido de un soporte
rígido (por ejemplo un techo), como se muestra en la figura
5.1b. Cuando M se reemplaza por un cuerpo diferente Mi, el
alargamiento del resorte será, por supuesto, distinto.
Por la Ley de Hooke, el resorte mismo ejerce una fuerza de
restitución F opuesta a la dirección del alargamiento y
proporcional a su magnitud s. Dicho en términos simples, F
= ks, en donde k es una constante de proporcionalidad.
Aunque cuerpos de distinto peso producen distintos
alargamientos del resorte, tal elemento elástico esta
esencialmente caracterizado por él numero k. Por ejemplo,
si un cuerpo que pesa 10lb. Alarga el resorte en 1/2 pie,
entonces,
10 = k (1/2) implica que k = 20 lb. /pie.
Luego, necesariamente una masa que pesa 8 lb. alarga el
mismo resorte en 2/5 pie.
Segunda Ley de Newton:
Después que una masa M se sujeta a un resorte, aquella lo
alargara en una magnitud s y alcanzara la posición de
equilibrio en la cual su peso W es equilibrado por la fuerza
de restitución ks. El peso es definido por:
W=m.g
En donde la masa puede medirse en Kilogramos, gramos o
geolibras (slugs) y g = 9.8 mt/s², p80 cm/s² o 32pie/s²,
respectivamente. Tal como se indica la figura 5.2b,la
condición de equilibrio es m.g = ks o bien m.g - ks = 0. Si
ahora la masa se desplaza de su posición de equilibrio en
una magnitud x y después se suelta, la fuerza
neta F correspondiente a este caso dinámico está dada por
la segunda ley del movimiento de Newton, F = ma, en
donde a es la aceleración d²w/dt². Suponiendo que sobre el
sistema no actúan fuerzas exteriores (movimiento
vibratorio libre), entonces podemos igualar F a la resultante
del peso y la fuerza de restitución:
m d²x/dt² = - k (s + x) + mg
= - kx + mg - ks = - kx
Cero
 Ecuación Diferencial Del Movimiento Libre no
Amortiguado:
Dividiendo la ultima ecuación planteada entre la masa m, se
obtiene la ecuación diferencial de segundo orden:
d²x/dt² + k/m x = 0
o bien d²x/dt² + .²x = 0
En donde .² = k/m. Se dice que la ecuación d²x/dt² + .²x = 0
describe el movimiento armónico simple o movimiento
vibratorio no amortiguado. Hay dos condiciones iniciales
obvias asociadas con dicha ecuación:
x(0) = ., dx/dt% = .
%t = 0
Que representa la magnitud del desplazamiento inicial y la
velocidad inicial, respectivamente. Por ejemplo si . > 0 y . <
0, se trata de una masa que parte de un punto abajo de la
posición de equilibrio y a la cual se ha comunicado una
velocidad dirigida hacia arriba. Si . < 0 y. > 0, se trata de una
masa en reposo que se suelta desde un punto que está %.
%unidades arriba de la posición de equilibrio. Los demás
casos son análogos.
Solución y ecuación de movimiento:
Para resolver la ecuación d²x/dt² + .²x = 0 observemos que
las soluciones de la ecuación auxiliar M² - w² = 0 son los
números complejo M = .i y Mi = - .i. De esta forma se
obtiene una solución general: x (t) = C1 cos .t + C2 sen .t.
El periodo de las vibraciones libres descritas por la ultima
ecuación general planteada es T = 2./. y la frecuencia es . =
1/T = ./2.. Por ejemplo, para x (t) = 2 cos 3t - 4 sen 3t el
periodo es 2./3 y la frecuencia es 3/2.. El primer numero
indica que hay 3 ciclos de la grafica de cada 2. Unidades; en
otras palabras, la masa realiza 3/2. Oscilaciones completas
por unidad de tiempo. Además, se puede demostrar que el
periodo 2. /. Es el intervalo de tiempo entre dos máximos
sucesivos de x (t). Finalmente, una vez que hemos
determinado las constantes C1 y C2 en x (t) = C1 cos .t + C2
sen .t mediante las condiciones iniciales
x(0) = . dx/dt% = .
%t = 0
, Decimos que la solución particular resultante es
la ecuación de movimiento.
Ejemplo:
Resolver e interpretar el problema de valor inicial:
D²x/dt² + 16 x = 0
x(0) = 10, dx/dt% = 0
%t = 0
Solución:
Una formulación equivalente del problema es: se estira
hacia abajo de un cuerpo que pende de un resorte hasta
que esté 10 unidades bajo la posición de equilibrio y luego
se le retiene hasta t = 0; se le suelta a continuación de
manera que parta de un estado de reposo. Aplicando las
condiciones iniciales a la solución:
x (t) = C1 cos 4t + C2 sen 4t.
Resulta x (0) = 10 = C1 . 1 + C2 . 0
de modo que C1 = 10 y por lo tanto
x (t) = 10 cos 4t + C2 sen 4t.
dx/dt = 40 sen 4t + 4C2 cos 4t
dx/dt% = 0 = 4C2 . 1
%t = 0
La ultima ecuación implica que C = 0 y por lo tanto la
ecuación de movimiento es x (t) = 10 cos 4t.
La solución muestra claramente que una vez que el sistema
se pone en movimiento, permanece en tal estado, con la
masa deslazándose alternadamente 10 unidades hacia cada
lado de la posición de equilibrio x = 0. El periodo de
oscilación es 2./4 = ./2 segundos.
Ejemplo:
Un cuerpo que pesa 2lb. se estira un resorte 6plg. Dicho cuerpo se
suelta en t = 0 desde un punto que está 8plg bajo la posición de
equilibrio, con una velocidad dirigida hacia arriba de 4/3 pie/seg.
Determine la función x (t) que describe el movimiento libre
resultante.
Solución:
Puesto que estamos usando el sistema de unidades inglesas
gravitatorias, las magnitudes dadas en pulgadas deben
expresarse en pies: 6plg = 6/12 = 1/2 pie, 8plg = 8/12 = 2/3
pie. Además, debemos convertir las unidades de peso en
unidades de masa. M = W/g
Tenemos M = 2/32 = 1/16slug; Además, por la Ley de Hooke
se tiene: 2 - k (1/2) lo que implica que k = 4lb/pie.
Por consiguiente, se tiene: 1/16 d²x/dt² = -4x y d²x/dt² + 64x
= 0.
El desplazamiento y la velocidad iniciales están dados por:
x(0) = 2/3, dx/dt% = - 4/3
%t = 0
En donde el signo negativo que aparece en la última
condición en consecuencia de que a la masa se le da una
velocidad inicial con dirección negativa, esto es, dirigida
hacia arriba.
Ahora bien, .² = 64, ósea. = 8, de modo que la solución
general de la ecuación diferencial es:
x (t) = C1 cos 8t + C2 sen 8t.
Aplicando las condiciones iniciales a esta ecuación tenemos
que:
x (0) = 2/3 = C1 . 1 + C2 . 0 (C1 = 2/3)
x (t) = 2/3 cos 8t + 8C2 cos 8t
x´(t) = - 16/3 sen 8t + 8C2 cos 8t
x´(0) = - 4/3 = - 16/3 . 0 + 8C2 .1, (C = -1/6)
Luego C2 = - 1/6. Por consiguiente, la ecuación de
movimiento es:
x (t) = 2/3 cos 8t - 1/6 sen 8t.
 CONCLUSIONES:
Finalmente y para concluir se determinó que, la resolución de problemas de
ingeniería está asociada, por lo general, a resultados numéricos puesto que se
requieren respuestas prácticas.
La mayor parte de las leyes científicas de expresan en términos de rapidez de
variación de una variable con respecto otra.
Además proporcionan una herramienta esencial para modelar muchos problemas
en Ingeniería, Física, Economía y Biología, puesto que estos, por lo general,
requieren la determinación de una función que satisface a una ecuación
diferencial.
El Método de Runge Kutta es mejor que el método de Euler, pero aún así es
posible aumentar la precisión achicando los pasos entre los puntos o
implementando el método de orden superior.
Es el método más utilizado para resolver numéricamente problemas de
ecuaciones diferenciales ordinarias con condiciones iniciales es el método de
Runge-Kutte, el cual proporciona un pequeño margen de error con respecto a la
solución real del problema y es fácilmente programable en un software para
realizar iteraciones necesarias.
El dominio de los métodos numéricos, en combinación con las capacidades y
potencialidades de la programación de computadoras resuelve problemas de
ingeniería de manera más fácil y eficientemente.