D.1: Solvers ODE simples — Derivación
- Page ID
- 119320
El primer orden del día es derivar tres algoritmos simples para generar soluciones numéricas aproximadas al problema de valor inicial
\ begin {align*} y' (t) &=f\ big (t, y (t)\ big)\\ y (t_0) &=y_0\ end {align*}
El primero se llama método de Euler porque fue desarrollado por (¡sorpresa!) Euler 1.
Método de Euler
Nuestro objetivo es aproximar (numéricamente) la función desconocida
\ begin {align*} y (t) &= y (t_0) +\ int_ {t_0} ^t y' (\ tau)\,\, d {\ tau}\\ &= y (t_0) +\ int_ {t_0} ^t f\ grande (\ tau, y (\ tau)\ grande)\,\, d {\ tau}\ end {align*}
para Se\(t\ge t_0\text{.}\) nos dice explícitamente el valor de\(y(t_0)\text{,}\) a saber\(y_0\text{.}\) Así que sabemos\(f\big(\tau,y(\tau)\big)\big|_{\tau=t_0}=f\big(t_0,y_0\big)\text{.}\) Pero no conocemos el\(f\big(\tau,y(\tau)\big)\) integrando para\(\tau>t_0\text{.}\) Por otro lado, si\(\tau\) está cerca\(t_0\text{,}\) entonces\(f\big(\tau,y(\tau)\big)\) permanecerá cerca 2 Este será el caso como siempre y cuando\(f\big(t_0,y_0\big)\text{.}\) Así\(f(t,y)\) is continuous. que elija un número pequeño\(h\) y defina
\ begin {align*} t_1&=t_0+h\\ y_1&=y (t_0) +\ int_ {t_0} ^ {t_1} f (t_0, y_0)\,\, d {\ tau} =y_0+f\ big (t_0, y_0\ big) (t_1-t_0)\\ &=y_0+f\ big (t_0, y_0\ grande) h\ end {alinear*}
Por el argumento anterior
\[ y(t_1)\approx y_1 \nonumber \]
Ahora empezamos de nuevo desde el nuevo punto Ahora\((t_1,y_1)\text{.}\) sabemos un valor aproximado para\(y\) en el momento\(t_1\text{.}\) Si\(y(t_1)\) fuera exactamente\(y_1\text{,}\) entonces la tasa instantánea de cambio de\(y\) en el momento a\(t_1\text{,}\) saber\(y'(t_1)=f\big(t_1,y(t_1)\big)\text{,}\) sería exactamente\(f(t_1,y_1)\) y\(f\big(t,y(t)\big)\) permanecería cerca \(f(t_1,y_1)\)para\(t\) cerca de\(t_1\text{.}\) Definir
\ begin {align*} t_2&=t_1+h=t_0+2h\\ y_2&=y_1 +\ int_ {t_1} ^ {t_2} f (t_1, y_1)\,\, d {t} =y_1+f\ big (t_1, y_1\ big) (t_2-t_1)\ &=y_1+f\ grande (t_1, y_1 grande\) h\ end {alinear*}
tenemos
\[ y(t_2)\approx y_2 \nonumber \]
Simplemente repetimos este argumento ad infinitum. Definir, para\(n=0,1,2,3,\cdots\)
\[ t_n=t_0+nh \nonumber \]
Supongamos que, para algún valor de ya\(n\text{,}\) hemos calculado un valor aproximado\(y_n\) para\(y(t_n)\text{.}\) Entonces la tasa de cambio de\(y(t)\) para\(t\) cerca de\(t_n\) es\(f\big(t,y(t)\big)\approx f\big(t_n,y(t_n)\big)\approx f\big(t_n,y_n\big)\) y
\[ y(t_{n+1})\approx y_{n+1}=y_n+f\big(t_n,y_n\big)h \nonumber \]
Este algoritmo se llama Método de Euler. El parámetro\(h\) se llama el tamaño del paso.
Aquí hay una tabla que aplica algunos pasos del método de Euler al problema del valor inicial
\ begin {align*} y'&=-2t+y\\ y (0) & = 3\ end {alinear*}
con tamaño de paso\(h=0.1\text{.}\) Para este problema de valor inicial
\ begin {align*} f (t, y) &=-2t+y\\ t_0&=0\\ y_0&=3\ end {align*}
Por supuesto, este problema de valor inicial ha sido elegido únicamente con fines ilustrativos. La solución exacta es 3\(y(t)=2+2t+e^t\text{.}\)
| \(n\) | \(t_n\) | \(y_n\) | \(f(t_n,y_n)=-2t_n+y_n\) | \(y_{n+1}=y_n+f(t_n,y_n)*h\) |
| 0 | 0.0 | 3.000 | \(-2*0.0+3.000=3.000\) | \(3.000+3.000*0.1=3.300\) |
| 1 | 0.1 | 3.300 | \(-2*0.1+3.300=3.100\) | \(3.300+3.100*0.1=3.610\) |
| 2 | 0.2 | 3.610 | \(-2*0.2+3.610=3.210\) | \(3.610+3.210*0.1=3.931\) |
| 3 | 0.3 | 3.931 | \(-2*0.3+3.931=3.331\) | \(3.931+3.331*0.1=4.264\) |
| 4 | 0.4 | 4.264 | \(-2*0.4+4.264=3.464\) | \(4.264+3.464*0.1=4.611\) |
| 5 | 0.5 | 4.611 |
La solución exacta en\(t=0.5\) es\(4.6487\text{,}\) a cuatro decimales. Esperamos que el método de Euler sea más preciso a medida que el tamaño del paso se haga más pequeño. Pero, por supuesto, la cantidad de esfuerzo sube también. Si calculamos usando\(h=0.01\text{,}\) obtenemos (después de mucho más trabajo)\(4.6446\text{.}\)
El método mejorado de Euler
El método de Euler es un algoritmo que genera soluciones aproximadas al problema del valor inicial
\ begin {align*} y' (t) &=f\ big (t, y (t)\ big)\\ y (t_0) &=y_0\ end {align*}
En aplicaciones,\(f(t,y)\) es una función dada y\(t_0\) y\(y_0\) se les dan números. Se\(y(t)\) desconoce la función. Denotar por\(\varphi(t)\) la solución exacta 4 Bajo hipótesis razonables sobre\(f\text{,}\) there is exactly one such solution. The interested reader should search engine their way to the Picard-Lindelöf theorem. para este problema de valor inicial. En otras palabras\(\varphi(t)\) es la función que obedece
\ begin {align*}\ varphi' (t) &=f\ big (t,\ varphi (t)\ big)\\ varphi (t_0) &=y_0\ end {align*}
exactamente.
Fijar un tamaño de paso\(h\) y definir\(t_n=t_0+nh\text{.}\) Al convertir el problema en una de aproximar integrales, ahora derivamos otro algoritmo que genera valores aproximados para\(\varphi\) en la secuencia de valores de tiempo igualmente espaciados\(t_0,\ t_1,\ t_2,\ \cdots\text{.}\) Vamos a denotar los valores aproximados\(y_n\) con
\[ y_n\approx\varphi(t_n) \nonumber \]
Por el teorema fundamental del cálculo y la ecuación diferencial, la solución exacta obedece
\ begin {alinear*}\ varphi (t_ {n+1}) &=\ varphi (t_n) +\ int_ {t_n} ^ {t_ {n+1}}\ varphi' (t)\\, d {t}\\ &=\ varphi (t_n) +\ int_ {t_n} ^ {t_ {n+1}} f\ grande (t,\ varphi (t)\ grande)\\, d {t}\ final {alinear*}
Arreglar cualquiera\(n\) y supongamos que ya hemos encontrado\(y_0,\ y_1,\ \cdots,\ y_n\text{.}\) Nuestro algoritmo para la computación\(y_{n+1}\) será de la forma
\[ y_{n+1}=y_n+\text{ approximate value of } \int_{t_n}^{t_{n+1}}f\big(t,\varphi(t)\big)\ \, d{t} \nonumber \]
En el método de Euler, nos aproximamos\(f\big(t,\varphi(t)\big)\)\(t_n\le t\le t_{n+1}\) por la constante\(f\big(t_n,y_n\big)\text{.}\) Así
\ begin {align*} &\ text {valor aproximado de Euler para}\ int_ {t_n} ^ {t_ {n+1}} f\ big (t,\ varphi (t)\ big)\\, d {t}\ text {es}\\ &\ hskip0.6in\ int_ {t_n} ^ {t_ {n+1}} f\ big (t_n, y_n, y_n\ grande)\\, d {t} =f\ grande (t_n, y_n\ grande) h\ final {alinear*}
Entonces Euler aproxima el área de la región complicada\(\ 0\le y\le f\big(t,\varphi(t)\big)\text{,}\)\(t_n\le t\le t_{n+1}\ \) (representada por la región sombreada debajo de la parábola en la mitad izquierda de la figura de abajo) por el área del rectángulo\(\ 0\le y\le f\big(t_n,y_n\big)\text{,}\)\(t_n\le t\le t_{n+1}\ \) (el rectángulo sombreado en la mitad derecha de la figura de abajo).
Nuestro segundo algoritmo, el método mejorado de Euler, obtiene una mejor aproximación usando la regla trapezoidal. Es decir, aproximamos la integral por el área del trapecio a la derecha de abajo, en lugar del rectángulo a la derecha de arriba.
El área exacta de este trapecio es la longitud\(h\) de la base multiplicada por el promedio de las alturas de los dos lados, que es Por\(\frac{1}{2}\big[f\big(t_n,\varphi(t_n)\big)+f\big(t_{n+1},\varphi(t_{n+1})\big)\big]\text{.}\) supuesto que no sabemos\(\varphi(t_n)\) o\(\varphi(t_{n+1})\) exactamente.
Recordemos que ya hemos encontrado\(y_0,\cdots,y_n\) y estamos en proceso de encontrar\(y_{n+1}\text{.}\) Así que ya tenemos una aproximación para\(\varphi(t_n)\text{,}\) saber\(y_n\text{.}\) Pero aún necesitamos aproximarnos ¡\(\varphi(t_{n+1})\text{.}\)Podemos hacerlo usando un paso del método original de Euler! Eso es
\[ \varphi(t_{n+1})\approx \varphi(t_n)+\varphi'(t_n)h \approx y_n+f(t_n,y_n)h \nonumber \]
Entonces nuestra aproximación de\(\frac{1}{2}\big[f\big(t_n,\varphi(t_n)\big)+f\big(t_{n+1},\varphi(t_{n+1})\big)\big]\) es
\[ \frac{1}{2}\Big[f\big(t_n,y_n\big)+f\Big(t_{n+1},y_n+f(t_n,y_n)h\Big)\Big] \nonumber \]
y
\ begin {align*} &\ text {Se mejoró el valor aproximado de Euler para}\ int_ {t_n} ^ {t_ {n+1}} f\ big (t,\ varphi (t)\ big)\\, d {t}\ text {is}\\ &\ hskip0.7in\ frac {1} {2}\ grande [f\ big (t_n, y_n big\) +f\ Grande (t_ {n+1}, y_n+f (t_n, y_n) h\ Grande)\ Grande] h\ final {alinear*}
Armando todo 5, el algoritmo mejorado del método de Euler es
\[ y(t_{n+1})\approx y_{n+1}=y_n+ \frac{1}{2}\Big[f\big(t_n,y_n\big)+f\Big(t_{n+1},y_n+f(t_n,y_n)h\Big)\Big]h \nonumber \]
Estos son los dos primeros pasos del método mejorado de Euler aplicado a
\ begin {align*} y'&=-2t+y\ qquad y (0) = 3\ cr\ end {alinear*}
con\(h=0.1\text{.}\) En cada paso calculamos\(f(t_n,y_n)\text{,}\) seguido por el\(y_n+f(t_n,y_n)h\text{,}\) cual denotamos\(\tilde y_{n+1}\text{,}\) seguido de\(f(t_{n+1},\tilde y_{n+1})\text{,}\) seguido de\(y_{n+1}=y_n+ \frac{1}{2}\big[f\big(t_n,y_n\big)+f\big(t_{n+1},\tilde y_{n+1}\big)\big]h\text{.}\)
\ begin {alignat*} {8} t_0&=0 & y_0&=3 & &\ implica & f (t_0, y_0) &=-2*0+3 =3\\ & & & & & & & & & & & implica &\ tilde y_1&=3+3*0.1 =3.3\\ & & & & & & & implica & f (t_1,\ tilde y_1) &=-2*0.1+3.3 =3.3 3.1\\ & & & & & & &\ implica & y_1& amp; =3+\ frac {1} {2} [3+3.1] *0.1 =3.305\ cr t_1&=0.1\ quad & y_1&=3.305 & & &\ implica & f (t_1, y_1) &=-2*0.1+3.305 =3.105\ & & & & & & & implica &\ tilde y_2&=3.305+3.105*0.1 =3.6155\ & & & & &\ implica & f (t_2,\ tilde y_2) &=-2*0.2+3.6155 = 3.2155\\ & && & & &\ implica & y_2&=3.305+\ frac {1} {2} [3.105+3.2155] *0.1\\ & && & & & & & implica & & = 3.621025\ end {alignat*}
Aquí hay una tabla que da los primeros cinco pasos.
| \(n\) | \(t_n\) | \(y_n\) | \(f(t_n,y_n)\) | \(\tilde y_{n+1}\) | \(f(t_{n+1},\tilde y_{n+1})\) | \(y_{n+1}\) |
| 0 | 0.0 | 3.000 | 3.000 | 3.300 | 3.100 | 3.305 |
| 1 | 0.1 | 3.305 | 3.105 | 3.616 | 3.216 | 3.621 |
| 2 | 0.2 | 3.621 | 3.221 | 3.943 | 3.343 | 3.949 |
| 3 | 0.3 | 3.949 | 3.349 | 4.284 | 3.484 | 4.291 |
| 4 | 0.4 | 4.291 | 3.491 | 4.640 | 3.640 | 4.647 |
| 5 | 0.5 | 4.647 |
Como vimos al final de la Sección D.1.1, el exacto\(y(0.5)\) es 4.6487, a cuatro decimales, y el método de Euler dio 4.611.
El método Runge-Kutta
El algoritmo Runge-Kutta 6 es similar al de Euler y mejorado los métodos de Euler en que también utiliza, en la notación de la última subsección,
\[ y_{n+1}=y_n+{\rm\ approximate\ value\ for \ } \int_{t_n}^{t_{n+1}}f\big(t,\varphi(t)\big)\ dt \nonumber \]
Pero en lugar de aproximarse\(\int_{t_n}^{t_{n+1}}f\big(t,\varphi(t)\big)\ dt\) por el área de un rectángulo, como lo hace Euler, o por el área de un trapecio, como lo hace Euler mejorado, se aproxima por el área bajo una parábola. Es decir, usa la regla de Simpson. De acuerdo con la regla de Simpson (que se deriva en la Sección 1.11.3)
\ begin {align*} &\ int_ {t_n} ^ {t_n+h} f\ grande (t,\ varphi (t)\ grande)\ dt\ &\ hskip0.5in\ approx\ tfrac {h} {6}\ grande [f\ grande (t_n,\ varphi (t_n)\ grande) +4f\ grande (t_n+\ tfrac {h} {2},\ varphi (t_n+\ tfrac {h} {2})\ grande) +f\ grande (t_n+h,\ varphi (t_n+h)\ grande)\ grande]\ final {alinear*}
Análogamente a lo que sucedió en nuestro desarrollo del método mejorado de Euler, no sabemos\(\varphi(t_n)\text{,}\)\(\varphi(t_n+\tfrac{h}{2})\) o\(\varphi(t_n+h)\text{.}\) Entonces tenemos que aproximarlos también. El algoritmo Runge-Kutta, incorporando todas estas aproximaciones, es 7 Está mucho más allá de nuestro alcance derivar este algoritmo, aunque la derivación es similar en sabor a la del método mejorado de Euler. Puedes encontrar más en, por ejemplo, Wikipedia.
\[\begin{align*} k_{1,n}&=f(t_n,y_n) \\ k_{2,n}&=f(t_n+\tfrac{1}{2}h,y_n+\tfrac{h}{2}k_{1,n}) \\ k_{3,n}&=f(t_n+\tfrac{1}{2}h,y_n+\tfrac{h}{2}k_{2,n}) \\ k_{4,n}&=f(t_n+h,y_n+hk_{3,n}) \\ y_{n+1}&=y_n+\tfrac{h}{6}\left[k_{1,n}+2k_{2,n}+2k_{3,n}+k_{4,n}\right] \end{align*}\]
Es decir, Runge-Kutta usa
- \(k_{1,n}\)aproximarse\(f\big(t_n,\varphi(t_n)\big)=\varphi'(t_n)\text{,}\)
- ambos\(k_{2,n}\) y\(k_{3,n}\) aproximarse\(f\big(t_n+\tfrac{h}{2},\varphi(t_n+\tfrac{h}{2})\big) =\varphi'(t_n+\tfrac{h}{2})\text{,}\) y
- \(k_{4,n}\)aproximarse\(f\big(t_n+h,\varphi(t_n+h)\big)\text{.}\)
Estos son los dos primeros pasos del algoritmo Runge-Kutta aplicado a
\ begin {align*} y'&=-2t+y\ qquad y (0) = 3\ end {align*}
con\(h=0.1\text{.}\)
\ begin {alignat*} {2} t_0&=0 & y_0&=3\\ &\ implica & k_ {1,0} &=f (0,3) =-2*0+3 =3\\ &\ implica & &y_0+\ tfrac {h} {2} k_ {1,0} =3+0.05*3=3.15\\ &\ implica & k_ {2,0} &=f (0.05,3.15) =-2*0.05+3.15 =3.05\ &\ implica & &y_0+\ tfrac {h} {2} k_ {2,0} =3+0.05*3.05=3.1525\\ &\ implica & k_ {3,0} &=f (0.05,3.1525) =-2*0.05+3.1525 =3.0525\\ &\ implica & &y_0+hk_ {3,0} =3+0.1*3.0525=3.30525\\ &\ implica & k_ {4,0} &=f (0.1,3.30525) =-2*0.1+3.30525 =3.10525\\ &\ implica & y_1&=3+\ tfrac {0.1} {6} [3+2*3.05+2*3.0525+3.10525] =3.3051708\\ t_1&=0.1 & y_1&=3.3051708\\ &\ implica & k_ {1,1} &=f (0.1,3.3051708) =-2*0.1+3.3051708 =3.1051708\\ &\ implica & y_1+\ tfrac {h} {2} k_ {1,1} =3.3051708+0.05*3.1051708=3.4604293\ &\ implica & k_ {2,1} &=f (0.15,3.4604293) =-2*0.15+3.4604293 =3.1604293\\ &\ implica & &y_1+\ tfrac {h} {2} k_ {2,1} =3.3051708+0.05*3.1604293=3.4631923\\ &\ implica & k_ {3,1} &=f (0.15,3.4631923) =-2*0.15+3.4631923 =3.1631923\\ &\ implica & y_1+hk_ {3,1} =3.3051708+0,11cm *3.4631923=3.62149\\ &\ implica & k_ {4,1} &=f (0.2,3.62149) =-2*0.2+3.62149 =3.22149\\ &\ implica & y_ 2&=3.3051708+\ tfrac {0.1} {6} [3.1051708+2*3.1604293+\\ & & &\ hskip1.0in+2*3.1631923+3.22149] =3.6214025\\ t_2&=0.2 & y_2&=3.6214025\ end {alignat*}
Ahora bien, si bien esto puede parecer intimidante escrito en su totalidad así, hay que tener en cuenta que es bastante fácil escribir un programa para hacer esto. Aquí hay una tabla que da los primeros cinco pasos. Los datos sólo se dan a tres decimales a pesar de que el cálculo se ha hecho a muchos más.
| \(n\) | \(t_n\) | \(y_n\) | \(k_{1,n}\) | \(y_{n,1}\) | \(k_{2,n}\) | \(y_{n,2}\) | \(k_{3,n}\) | \(y_{n,3}\) | \(k_{4,n}\) | \(y_{n+1}\) |
| 0 | 0.0 | 3.000 | 3.000 | 3.150 | 3.050 | 3.153 | 3.053 | 3.305 | 3.105 | 3.305 |
| 1 | 0.1 | 3.305 | 3.105 | 3.460 | 3.160 | 3.463 | 3.163 | 3.621 | 3.221 | 3.621 |
| 2 | 0.2 | 3.621 | 3.221 | 3.782 | 3.282 | 3.786 | 3.286 | 3.950 | 3.350 | 3.949 |
| 3 | 0.3 | 3.950 | 3.350 | 4.117 | 3.417 | 4.121 | 3.421 | 4.292 | 3.492 | 4.291 |
| 4 | 0.4 | 4.292 | 3.492 | 4.466 | 3.566 | 4.470 | 3.570 | 4.649 | 3.649 | 4.648 |
| 5 | 0.5 | 4.6487206 |
Como vimos al final de la Sección D.1.2, el exacto\(y(0.5)\) es 4.6487213, a siete decimales, el método de Euler dio 4.611, y Euler mejorado dio 4.647.
Hasta el momento hemos, ojalá, motivado al Euler, mejorado los algoritmos de Euler y Runge-Kutta. No hemos intentado ver cuán eficientes y precisos son los algoritmos. En la siguiente sección se ofrece un primer vistazo a esas preguntas.