| 9. Ecuaciones diferenciales |

9.3. Métodos numéricos: Objetivos de aprendizaje

       En este capítulo estudiamos métodos numéricos para resolver una ecuación diferencial de primer orden.

La SECCIÓN 9.3.1 trata del método de Euler, que en realidad es demasiado tosco para ser de mucha utilidad en aplicaciones prácticas. Sin embargo, su simplicidad permite una introducción a las ideas necesarias para comprender los mejores métodos discutidos en las otras dos secciones.
La SECCIÓN 9.3.2 analiza las mejoras en el método de Euler.
La SECCIÓN 9.3.3 trata del método de Runge-Kutta, quizás el método más utilizado para la solución numérica de ecuaciones diferenciales.

9.3.1 MÉTODO DE EULER

Si un problema de valor inicial

y′ = f (x, y),  y(x0) = y0        (9.3.1.1)

no se puede resolver analíticamente, es necesario recurrir a métodos numéricos para obtener aproximaciones útiles a una solución de (9.3.1.1). Consideraremos estos métodos en este capítulo.

        Estamos interesados en calcular valores aproximados de la solución de (9.3.1.1) en puntos igualmente espaciados x0, x1,. . . , xn = b en un intervalo [x0, b]. Por lo tanto,

xi = x0 + ih, i = 0, 1, . . ., n,

donde

Denotaremos los valores aproximados de la solución en estos puntos por y0, y1,. . . , yn; por tanto, yi es una aproximación de y(xi). Llamaremos

ei = y(xi) − yi

el error en el i-ésimo paso. Debido a la condición inicial y(x0) = y0, siempre tendremos e0 = 0. Sin embargo, en general ei ≠ 0 si i > 0.

       Encontramos dos fuentes de error al aplicar un método numérico para resolver un problema de valor inicial:
• Las fórmulas que definen el método se basan en algún tipo de aproximación. Los errores debidos a la inexactitud de la aproximación se denominan errores de truncamiento.
• Las computadoras hacen aritmética con un número fijo de dígitos y, por lo tanto, cometen errores al evaluar las fórmulas que definen los métodos numéricos. Los errores debidos a la incapacidad de la computadora para realizar operaciones aritméticas exactas se denominan errores de redondeo.
       Dado que un análisis cuidadoso del error de redondeo está más allá del alcance de este libro, consideraremos solo los errores de truncamiento.

Método de Euler

El método numérico más simple para resolver (9.3.1.1) es el método de Euler. Este método es tan tosco que rara vez se utiliza en la práctica; sin embargo, su simplicidad lo hace útil con fines ilustrativos.

      El método de Euler se basa en el supuesto de que la recta tangente a la curva integral de (9.3.1.1) en (xi, y(xi)) se aproxima a la curva integral en el intervalo [xi, xi + 1]. Dado que la pendiente de la curva integral de (9.3.1.1) en (xi, y(xi)) es y′(xi) = f (xi, y(xi)), la ecuación de la recta tangente a la curva integral en (xi, y(xi)) es

y = y(xi) + f (xi, y(xi)) (xxi).        (9.3.1.2)

Al establecer x = xi + 1 = xi + h en (9.3.1.2) se obtiene

yi + 1 = y(xi) + h f (xi, y(xi))        (9.3.1.3)

como una aproximación a y(xi + 1). Como se conoce y(x0) = y0, podemos usar (9.3.1.3) con i = 0 para calcular

y1 = y0 + h f (x0, y0).

Sin embargo, al establecer i = 1 en (9.3.1.3) se obtiene

y2 = y(x1) + h f (x1, y(x1)),

lo cual no es útil, ya que no sabemos el valor de y(x1). Por lo tanto, reemplazamos y(x1) por su valor aproximado y1 y redefinimos

y2 = y1 + h f (x1, y1).

Habiendo calculado y2, podemos calcular

y3 = y2 + h f (x2, y2).

En general, el método de Euler comienza con el valor conocido y(x0) = y0 y calcula y1, y2,. . . , yn sucesivamente con la fórmula

yi + 1 = yi + h f (xi, yi),  0 ≤ in − 1.      (9.3.1.4)

    El siguiente ejemplo ilustra el procedimiento de cálculo indicado en el método de Euler.

Ejemplo ilustrativo 9.3.1.1

Use el método de Euler con h = 0.1 para encontrar valores aproximados para la solución del problema de valor inicial

       (9.3.1.5)

en x = 0.1, 0.2, 0.3.

Solución:

Reescribimos (9.3.1.5) como

Esta imagen tiene un atributo ALT vacío; su nombre de archivo es image-177.png

que es de la forma (9.3.1.1), con

Esta imagen tiene un atributo ALT vacío; su nombre de archivo es image-178.png

aplicando el método de Euler, se obtiene

Esta imagen tiene un atributo ALT vacío; su nombre de archivo es image-179.png

      Hemos escrito los detalles de estos cálculos para asegurarnos de que comprenda el procedimiento. Sin embargo, en el resto de los ejemplos, así como en los ejercicios de este capítulo, asumiremos que puede utilizar una calculadora programable o una computadora para realizar los cálculos necesarios.

Ejemplos que ilustran el error en el método de Euler

Ejemplo ilustrativo 9.3.1.2

Utilice el método de Euler con tamaños de paso h = 0.1, h = 0.05 y h = 0.025 para encontrar valores aproximados de la solución del problema de valor inicial.en x = 0, 0,1, 0,2, 0,3,. . . , 1.0. Compare estos valores aproximados con los valores de la solución exacta

Esta imagen tiene un atributo ALT vacío; su nombre de archivo es image-173.png        (9.3.1.6

que se puede obtener mediante el método de la Sección 9.2.1. (Verificar)

Solución:

La tabla 9.3.1.1 muestra los valores de la solución exacta (9.3.1.6) en los puntos especificados y los valores aproximados de la solución en estos puntos obtenidos por el método de Euler con tamaños de paso h = 0.1, h = 0.05 y h = 0.025 . Al examinar esta tabla, tenga en cuenta que los valores aproximados en la columna correspondiente a h = .05 son en realidad los resultados de 20 pasos con el método de Euler. No hemos enumerado las estimaciones de la solución obtenida para x = 0.05, 0.15,. . . , ya que no hay nada con qué compararlos en la columna correspondiente a h = 0.1. De manera similar, los valores aproximados en la columna correspondiente a h = 0.025 son en realidad los resultados de 40 pasos con el método de Euler.

Esta imagen tiene un atributo ALT vacío; su nombre de archivo es image-174.png

(Tabla 9.3.1.1. Solución numérica de y′ + 2y = x3e − 2x, y(0) = 1, por el método de Euler.)

Puede ver en la tabla 9.3.1.1 que disminuir el tamaño del paso mejora la precisión del método de Euler.
Por ejemplo,

Con base en esta escasa evidencia, podría adivinar que el error al aproximar la solución exacta a un valor fijo de x por el método de Euler se reduce aproximadamente a la mitad cuando el tamaño del paso se reduce a la mitad. Puede encontrar más evidencia para apoyar esta conjetura examinando la tabla 9.3.1.2, que enumera los valores aproximados de yexactoyaproximado en x = 0.1, 0.2,. . . , 1.0.

(Tabla 9.3.1.2. Errores en soluciones aproximadas de y′ + 2y = x3e − 2x, y(0) = 1, obtenido por el método de Euler.)

Ejemplo ilustrativo 9.3.1.3

Las tablas 9.3.1.3 y 9.3.1.4 muestran resultados análogos para el problema de valor inicial no lineal

y′ = −2y2 + xy + x2y(0) = 1,        (9.3.1.7)

excepto que en este caso no podemos resolver (9.3.1.7) exactamente. Los resultados en la columna “Exacto” se obtuvieron utilizando un método numérico más preciso conocido como el método de Runge-Kutta con un tamaño de paso pequeño. Son exactos a ocho lugares decimales. 

(Tabla 9.3.1.3. Solución numérica de y′ = −2y2 + xy + x2y(0) = 1, por el método de Euler.)
(Tabla 9.3.1.4. Errores en soluciones aproximadas de y′ = −2y2 + xy + x2y(0) = 1, obtenidas por el método de Euler.)

      Dado que creemos que es importante para evaluar la precisión de los métodos numéricos que estudiaremos en este capítulo, a menudo incluimos una columna que enumera los valores de la solución exacta del problema de valor inicial, incluso si las instrucciones del ejemplo o del ejercicio no lo piden específicamente. Si se incluyen comillas en el encabezado, los valores se obtuvieron aplicando el método de Runge-Kutta de la manera que se explica en la Sección 9.3.3. Si no se incluyen las comillas, los valores se obtuvieron de una fórmula conocida para la solución. En cualquier caso, los valores son exactos hasta ocho lugares a la derecha del punto decimal.

Error de truncamiento en el método de Euler

De acuerdo con los resultados indicados en las tablas 9.3.1.1 a 9.3.1.4, ahora mostraremos que, bajo supuestos razonables sobre f, existe una constante K tal que el error al aproximar la solución del problema de valor inicial

y′ = f(x, y), y(x0) = y0,

en un punto dado b > x0 por el método de Euler con tamaño de paso h = (bx0)/n satisface la desigualdad

|y(b) − yn| ≤ Kh,

donde K es una constante independiente de n.

 

       Hay dos fuentes de error (sin contar el redondeo) en el método de Euler:

1. El error cometido al aproximar la curva integral por la recta tangente y = y(xi) + f (xi, y(xi)) (xxi) (9.3.1.2) sobre el intervalo [xi, xi+1].

2. El error cometido al reemplazar y(xi) por yi en (9.3.1.2) y usar yi + 1 = yi + h f (xi, yi),  0 ≤ in − 1 (9.3.1.4) en lugar de (9.3.1.2) para calcular yi+1.

 

      El método de Euler asume que yi+1 definido en (9.3.1.2) es una aproximación a y(xi+1). Llamamos al error en esta aproximación el error de truncamiento local en el i-ésimo paso, y lo denotamos por Ti; por lo tanto,

Ti = y(xi+1) − y(xi) − hf(xi, y(xi)).      (9.3.1.8)

Ahora usaremos el teorema de Taylor para estimar Ti, asumiendo por simplicidad que f, fx y fy son continuas y acotadas para todo (x, y). Entonces y″ existe y está acotado en [x0, b]. Para ver esto, diferenciamos

y′(x) = f(x, y(x))

para obtener

Esta imagen tiene un atributo ALT vacío; su nombre de archivo es image-124.png

Como asumimos que f, fx y fy están acotadas, existe una constante M tal que

Esta imagen tiene un atributo ALT vacío; su nombre de archivo es image-125.pnglo que implica que

y″(x)| ≤ M, x0 < x < b.        (9.3.1.9)

Dado que xi+1 = xi + h, el teorema de Taylor implica que

Esta imagen tiene un atributo ALT vacío; su nombre de archivo es image-126.pngdonde Esta imagen tiene un atributo ALT vacío; su nombre de archivo es image-127.pnges un número entre xi y xi+1. Dado que y′(xi) = f(xi, y(xi)) esto se puede escribir como

Esta imagen tiene un atributo ALT vacío; su nombre de archivo es image-129.png

o equivalente,

Esta imagen tiene un atributo ALT vacío; su nombre de archivo es image-130.png

Comparando esto con (9.3.1.8) se observa que

Esta imagen tiene un atributo ALT vacío; su nombre de archivo es image-131.png

Recordando (9.3.1.9), podemos establecer el límite

Esta imagen tiene un atributo ALT vacío; su nombre de archivo es image-132.png       (9.3.1.10)

Aunque puede ser difícil determinar la constante M, lo importante es que existe una M tal que (9.3.1.10) se cumple. Decimos que el error de truncamiento local del método de Euler es de orden h2, que escribimos como O(h2).

      Tenga en cuenta que la magnitud del error de truncamiento local en el método de Euler está determinada por la segunda derivada y″ de la solución del problema de valor inicial. Por lo tanto, el error de truncamiento local será mayor donde |y″| es grande, o menor donde |y″| es pequeño.

      Dado que el error de truncamiento local para el método de Euler es O(h2), es razonable esperar que reducir a la mitad h reduzca el error de truncamiento local en un factor de 4. Esto es cierto, pero reducir a la mitad el tamaño del paso también requiere el doble de pasos para aproximar el solución en un punto dado. Para analizar el efecto general del error de truncamiento en el método de Euler, es útil derivar una ecuación que relacione los errores

ei+1 = y(xi+1) − yi+1  y  ei = y(xi) − yi.

Para ello, recuerda que

Esta imagen tiene un atributo ALT vacío; su nombre de archivo es image-133.png      (9.3.1.11)

y

Esta imagen tiene un atributo ALT vacío; su nombre de archivo es image-134.png       (9.3.1.12)

Restando (9.3.1.12) de (9.3.1.11) se obtiene

Esta imagen tiene un atributo ALT vacío; su nombre de archivo es image-135.png       (9.3.1.13)

El último término de la derecha es el error de truncamiento local en el i-ésimo paso. Los otros términos reflejan la forma en que los errores cometidos en los pasos anteriores afectan a ei+1. Desde |Ti| ≤ Mh2/2, vemos de (9.3.1.13) que

Esta imagen tiene un atributo ALT vacío; su nombre de archivo es image-136.png      (9.3.1.14)

Como asumimos que fy es continua y acotada, el teorema del valor medio implica que

Esta imagen tiene un atributo ALT vacío; su nombre de archivo es image-137.png

dondeEsta imagen tiene un atributo ALT vacío; su nombre de archivo es image-138.pngestá entre yi y y(xi). Por lo tanto

Esta imagen tiene un atributo ALT vacío; su nombre de archivo es image-139.png

para alguna R constante. De esto y (9.3.1.14),

Esta imagen tiene un atributo ALT vacío; su nombre de archivo es image-140.png      (9.3.1.15)

Por conveniencia, sea C = 1 + Rh. Como e0 = y(x0) − y0 = 0, aplicando (9.3.1.15) repetidamente se obtiene

Esta imagen tiene un atributo ALT vacío; su nombre de archivo es image-142.png

Recordando la fórmula para la suma de una serie geométrica, vemos que

Esta imagen tiene un atributo ALT vacío; su nombre de archivo es image-143.png

(ya que C = 1 + Rh). De esto y (9.3.1.16),

Esta imagen tiene un atributo ALT vacío; su nombre de archivo es image-144.png          (9.3.1.17)

Dado que el teorema de Taylor implica que

1 + Rh < eRh

(verificar),

Esta imagen tiene un atributo ALT vacío; su nombre de archivo es image-145.png

Esto y (9.3.1.17) implican que

y(b) − yn| ≤ Kh,         (9.3.1.18)

con

Esta imagen tiene un atributo ALT vacío; su nombre de archivo es image-146.png

Por (9.3.1.18) decimos que el error de truncamiento global del método de Euler es de orden h, que escribimos como O(h).

Ecuaciones Semilineales y Variación de Parámetros

Una ecuación que se puede escribir en la forma

y′ + p(x)y = h(x, y)         (9.3.1.19)

con pEsta imagen tiene un atributo ALT vacío; su nombre de archivo es image-147.png0 se dice que es semilineal. (Por supuesto, (9.3.1.19) es lineal si h es independiente de y). Una forma de aplicar el método de Euler a un problema de valor inicial

y′ + p(x)y = h(x, y), y(x0) = y0        (9.3.1.20)

para (9.3.1.19) es pensar en ello como

y′ = f(x, y), y(x0) = y0,

donde

f(x, y) = −p(x)y + h(x, y).

Sin embargo, también podemos comenzar aplicando variación de parámetros a (9.3.1.20), como en las Secciones 9.2.1 y 9.2.4; por lo tanto, escribimos la solución de (9.3.1.20) como y = uy1, donde y1 es una solución no trivial de la ecuación complementaria y′ + p(x)y = 0. Entonces y = y = uy1 es una solución de (9.3.1.20) si y solo si u es una solución del problema de valor inicial

u′ = h(x, uy1(x))/y1(x), u(x0) = y(x0)/y1(x0).         (9.3.1.21)

Podemos aplicar el método de Euler para obtener valores aproximados u0, u1, . . . , un de este problema de valor inicial, y luego tomar

yi = uiy1(xi)

como valores aproximados de la solución de (9.3.1.20). Llamaremos a este procedimiento el método semilineal de Euler.

 

      Los siguientes dos ejemplos muestran que los métodos semilineales de Euler y Euler pueden arrojar resultados drásticamente diferentes.

Ejemplo ilustrativo 9.3.1.4

En el Ejemplo 9.2.1.7 tuvimos que dejar la solución del problema de valor inicial

y′ − 2xy = 1, y(0) = 3        (9.3.1.22)

en la forma

Esta imagen tiene un atributo ALT vacío; su nombre de archivo es image-148.png        (9.3.1.23)

porque era imposible evaluar esta integral exactamente en términos de funciones elementales. Use tamaños de paso h = 0.2, h = 0.1 y h = 0.05 para encontrar valores aproximados de la solución de (9.3.1.22) en x = 0, 0.2, 0.4, 0.6, . . . , 2.0 por (a) el método de Euler; (b) el método semilineal de Euler.

Solución:

(a) Reescribiendo (9.3.1.22) como

y′ = 1 + 2xy, y(0) = 3        (9.3.1.24)

y aplicando el método de Euler con f(x, y) = 1 + 2xy se obtienen los resultados que se muestran en la tabla 9.3.1.5. Debido a las grandes diferencias entre las estimaciones obtenidas para los tres valores de h, sería evidente que estos resultados son inútiles incluso si los valores “exactos” no se incluyeran en la tabla.

Tabla 9.3.1.5. Solución numérica de y′ − 2xy = 1, y(0) = 3, con el método de Euler

Es fácil ver por qué el método de Euler produce resultados tan pobres. Recuerde que la constante M en Esta imagen tiene un atributo ALT vacío; su nombre de archivo es image-132.png(9.3.1.10) – que juega un papel importante en la determinación del error de truncamiento local en el método de Euler –  debe ser un límite superior para los valores de la segunda derivada y″ de la solución del problema de valor inicial (9.3.1.22) en (0, 2). El problema es que y″ asume valores muy grandes en este intervalo. Para ver esto derivamos (9.3.1.24) para obtener

y″(x) = 2y(x) + 2xy′(x) = 2y(x) + 2x(1 + 2xy(x)) = 2(1 + 2x2)y(x) + 2x

donde la segunda igualdad sigue de nuevo de (9.3.1.24). Dado que (9.3.1.23) implica que Esta imagen tiene un atributo ALT vacío; su nombre de archivo es image-150.png si x > 0,

Esta imagen tiene un atributo ALT vacío; su nombre de archivo es image-151.png

Por ejemplo, haciendo que x = 2 muestra que y″(2) > 2952.

(b) Como Esta imagen tiene un atributo ALT vacío; su nombre de archivo es image-152.png es una solución de la ecuación complementaria y′ − 2xy = 0, podemos aplicar el método semilineal de Euler a (9.3.1.22), con

Esta imagen tiene un atributo ALT vacío; su nombre de archivo es image-153.png

Los resultados enumerados en la Tabla 9.3.1.6 son claramente mejores que los obtenidos por el método de Euler.

 

Tabla 9.3.1.6. Solución numérica de y′ − 2xy = 1, y(0) = 3, por el método semilineal de Euler

      No podemos dar un procedimiento general para determinar de antemano si el método de Euler o el método semilineal de Euler producirán mejores resultados para un problema de valor inicial semilineal dado y′ + p(x)y = h(x, y), y(x0) = y0 (9.3.1.20). Como regla general, el método semilineal de Euler dará mejores resultados que el método de Euler si |u″| es pequeño en [x0, b], mientras que el método de Euler da mejores resultados si |u″| es grande en [x0, b]. En muchos casos los resultados obtenidos por los dos métodos no difieren apreciablemente. Sin embargo, proponemos una forma intuitiva de decidir cuál es el mejor método: probar ambos métodos con varios tamaños de paso, como hicimos en el Ejemplo 9.3.1.4, y aceptar los resultados obtenidos por el método para el que las aproximaciones cambian menos a medida que disminuye el tamaño del paso.

Ejemplo ilustrativo 9.3.1.5

Aplicar el método de Euler con tamaños de paso h = 0.1, h = 0.05 y h = 0.025 al problema de valor inicial

Esta imagen tiene un atributo ALT vacío; su nombre de archivo es image-155.png      (9.3.1.25)

en [1, 2] arroja los resultados de la Tabla 9.3.1.7. Aplicando el método semilineal de Euler con

Esta imagen tiene un atributo ALT vacío; su nombre de archivo es image-156.png

da los resultados de la Tabla 9.3.1.8. Dado que estos últimos son claramente menos dependientes del tamaño del paso que los primeros, concluimos que el método semilineal de Euler es mejor que el método de Euler para (9.3.1.25). Esta conclusión se sustenta comparando los resultados aproximados obtenidos por los dos métodos con los valores “exactos” de la solución.

 

Tabla 9.3.1.7. Solución numérica de y′ − 2y = x/(1 + y2), y(1) = 7, por el método de Euler.

 

Tabla 9.3.1.8. Solución numérica de y′ − 2y = x/(1 + y2), y(1) = 7, por el método semilineal de Euler.

Ejemplo ilustrativo 9.3.1.6

Aplicar el método de Euler con tamaños de paso h = 0.1, h = 0.05 y h = 0.025 al problema de valor inicial

y′ + 3x2y = 1 + y2y(2) = 2        (9.3.1.26)

en [2, 3] arroja los resultados de la Tabla 9.3.1.9. Aplicando el método semilineal de Euler con

Esta imagen tiene un atributo ALT vacío; su nombre de archivo es image-159.png

da los resultados de la Tabla 9.3.1.10. Observando la estrecha concordancia entre las tres columnas de la Tabla 9.3.1.9 (al menos para valores mayores de x) y la falta de dicha concordancia entre las columnas de la Tabla 9.3.1.10, concluimos que el método de Euler es mejor que el método semilineal de Euler. para (9.3.1.26). La comparación de los resultados con los valores exactos apoya esta conclusión.

 

Tabla 9.3.1.9. Solución numérica de y′ + 3x2y = 1 + y2, y(2) = 2, por el método de Euler.

 

Tabla 9.3.1.10. Solución numérica de y′ + 3x2y = 1 + y2, y(2) = 2, por el método semilineal de Euler.

      En las próximas dos secciones estudiaremos otros métodos numéricos para resolver problemas de valor inicial, llamados el método mejorado de Euler, el método del punto medio, el método de Heun y el método de Runge-Kutta. Si el problema de valor inicial es semilineal como en (9.3.1.19), también tenemos la opción de usar la variación de parámetros y luego aplicar el método numérico dado al problema de valor inicial (9.3.1.21) para u. Por analogía con la terminología utilizada aquí, llamaremos al procedimiento resultante método semilineal de Euler mejorado, método semilineal del punto medio, método semilineal de Heun o método semilineal de Runge-Kutta, según sea el caso.


Deja un comentario

Tu dirección de correo electrónico no será publicada.