3.3: El método Runge-Kutta
- Page ID
- 114808
\( \newcommand{\vecs}[1]{\overset { \scriptstyle \rightharpoonup} {\mathbf{#1}} } \)
\( \newcommand{\vecd}[1]{\overset{-\!-\!\rightharpoonup}{\vphantom{a}\smash {#1}}} \)
\( \newcommand{\id}{\mathrm{id}}\) \( \newcommand{\Span}{\mathrm{span}}\)
( \newcommand{\kernel}{\mathrm{null}\,}\) \( \newcommand{\range}{\mathrm{range}\,}\)
\( \newcommand{\RealPart}{\mathrm{Re}}\) \( \newcommand{\ImaginaryPart}{\mathrm{Im}}\)
\( \newcommand{\Argument}{\mathrm{Arg}}\) \( \newcommand{\norm}[1]{\| #1 \|}\)
\( \newcommand{\inner}[2]{\langle #1, #2 \rangle}\)
\( \newcommand{\Span}{\mathrm{span}}\)
\( \newcommand{\id}{\mathrm{id}}\)
\( \newcommand{\Span}{\mathrm{span}}\)
\( \newcommand{\kernel}{\mathrm{null}\,}\)
\( \newcommand{\range}{\mathrm{range}\,}\)
\( \newcommand{\RealPart}{\mathrm{Re}}\)
\( \newcommand{\ImaginaryPart}{\mathrm{Im}}\)
\( \newcommand{\Argument}{\mathrm{Arg}}\)
\( \newcommand{\norm}[1]{\| #1 \|}\)
\( \newcommand{\inner}[2]{\langle #1, #2 \rangle}\)
\( \newcommand{\Span}{\mathrm{span}}\) \( \newcommand{\AA}{\unicode[.8,0]{x212B}}\)
\( \newcommand{\vectorA}[1]{\vec{#1}} % arrow\)
\( \newcommand{\vectorAt}[1]{\vec{\text{#1}}} % arrow\)
\( \newcommand{\vectorB}[1]{\overset { \scriptstyle \rightharpoonup} {\mathbf{#1}} } \)
\( \newcommand{\vectorC}[1]{\textbf{#1}} \)
\( \newcommand{\vectorD}[1]{\overrightarrow{#1}} \)
\( \newcommand{\vectorDt}[1]{\overrightarrow{\text{#1}}} \)
\( \newcommand{\vectE}[1]{\overset{-\!-\!\rightharpoonup}{\vphantom{a}\smash{\mathbf {#1}}}} \)
\( \newcommand{\vecs}[1]{\overset { \scriptstyle \rightharpoonup} {\mathbf{#1}} } \)
\( \newcommand{\vecd}[1]{\overset{-\!-\!\rightharpoonup}{\vphantom{a}\smash {#1}}} \)
En general, si\(k\) hay algún número entero positivo y\(f\) satisface supuestos apropiados, existen métodos numéricos con error de truncamiento local\(O(h^{k+1})\) para resolver un problema de valor inicial
\[\label{eq:3.3.1} y'=f(x,y),\quad y(x_0)=y_0.\]
Además, se puede demostrar que un método con error de truncamiento local\(O(h^{k+1})\) tiene un error de truncamiento global\(O(h^k)\). En las Secciones 3.1 y 3.2 estudiamos métodos numéricos donde\(k=1\) y\(k=2\). Saltaremos los métodos para los cuales\(k=3\) y procederemos al método Runge - Kutta, el método más utilizado, para el cual\(k=4\). La magnitud del error de truncamiento local está determinada por la quinta derivada\(y^{(5)}\) de la solución del problema del valor inicial. Por lo tanto, el error de truncamiento local será mayor donde\(|y^{(5)}|\) sea grande, o menor donde\(|y^{(5)}|\) sea pequeño. El método Runge-Kutta calcula valores aproximados\(y_1\),,...\(y_2\),\(y_n\) de la solución de la Ecuación\ ref {eq:3.3.1} at\(x_0\),,...\(x_0+h\),\(x_0+nh\) de la siguiente manera: Dado\(y_i\), computar
\[\begin{align*} k_{1i}&=f(x_i,y_i),\\ k_{2i}&=f \left(x_i+{h\over2},y_i+{h\over2}k_{1i}\right),\\ k_{3i}&=f\left(x_i+{h\over2},y_i+{h\over2}k_{2i}\right),\\ k_{4i}&=f(x_i+h,y_i+hk_{3i}),\end{align*}\]
y
\[y_{i+1}=y_i+{h\over6}(k_{1i}+2k_{2i}+2k_{3i}+k_{4i}).\nonumber \]
El siguiente ejemplo, que trata del problema de valor inicial considerado en Ejemplos y Ejemplo 3.3.1 , ilustra el procedimiento computacional indicado en el método Runge-Kutta.
Utilice el método Runge-Kutta con\(h=0.1\) para encontrar valores aproximados para la solución del problema del valor inicial
\[\label{eq:3.3.2} y'+2y=x^3e^{-2x},\quad y(0)=1,\]
en\(x=0.1,0.2\).
Solución
Nuevamente reescribimos la ecuación\ ref {eq:3.3.2} como
\[y'=-2y+x^3e^{-2x},\quad y(0)=1, \nonumber\]
que es de la forma Ecuación\ ref {eq:3.3.1}, con
\[f(x,y)=-2y+x^3e^{-2x},\ x_0=0,\mbox{ and}\ y_0=1. \nonumber\]
El método Runge-Kutta rinde
\[\begin{aligned} k_{10} & = f(x_0,y_0) = f(0,1)=-2,\\ k_{20} & = f(x_0+h/2,y_0+hk_{10}/2)=f(.05,1+(.05)(-2))\\ &= f(.05,.9)=-2(.9)+(.05)^3e^{-.1}=-1.799886895,\\ k_{30} & = f(x_0+h/2,y_0+hk_{20}/2)=f(.05,1+(.05)(-1.799886895))\\ &= f(.05,.910005655)=-2(.910005655)+(.05)^3e^{-.1}=-1.819898206,\\ k_{40} & = f(x_0+h,y_0+hk_{30})=f(.1,1+(.1)(-1.819898206))\\ &=f(.1,.818010179)=-2(.818010179)+(.1)^3e^{-.2}=-1.635201628,\\ y_1&=y_0+{h\over6}(k_{10}+2k_{20}+2k_{30}+k_{40}),\\ &=1+{.1\over6}(-2+2(-1.799886895)+2(-1.819898206) -1.635201628)=.818753803,\\[4pt] k_{11} & = f(x_1,y_1) = f(.1,.818753803)=-2(.818753803))+(.1)^3e^{-.2}=-1.636688875,\\ k_{21} & = f(x_1+h/2,y_1+hk_{11}/2)=f(.15,.818753803+(.05)(-1.636688875))\\ &= f(.15,.736919359)=-2(.736919359)+(.15)^3e^{-.3}=-1.471338457,\\ k_{31} & = f(x_1+h/2,y_1+hk_{21}/2)=f(.15,.818753803+(.05)(-1.471338457))\\ &= f(.15,.745186880)=-2(.745186880)+(.15)^3e^{-.3}=-1.487873498,\\ k_{41} & = f(x_1+h,y_1+hk_{31})=f(.2,.818753803+(.1)(-1.487873498))\\ &=f(.2,.669966453)=-2(.669966453)+(.2)^3e^{-.4}=-1.334570346,\\ y_2&=y_1+{h\over6}(k_{11}+2k_{21}+2k_{31}+k_{41}),\\ &=.818753803+{.1\over6}(-1.636688875+2(-1.471338457)+2(-1.487873498)-1.334570346) \\&=.670592417.\end{aligned}\]
El método Runge-Kutta es lo suficientemente preciso para la mayoría de las aplicaciones.
La Tabla 3.3.1 muestra los resultados del uso del método Runge-Kutta con tamaños de paso\(h=0.1\) y\(h=0.05\) para encontrar valores aproximados de la solución del problema del valor inicial
\[y'+2y=x^3e^{-2x},\quad y(0)=1 \nonumber\]
al\(x=0\),\(0.1\),\(0.2\),\(0.3\),...,\(1.0\). A modo de comparación, también se muestran los valores aproximados correspondientes obtenidos con el método mejorado de Euler en el Ejemplo 3.2.2, y los valores de la solución exacta
\[y={e^{-2x}\over4}(x^4+4).\nonumber \]
Los resultados obtenidos por el método Runge-Kutta son claramente mejores que los obtenidos por el método mejorado de Euler de hecho; los resultados obtenidos por el método Runge-Kutta con\(h=0.1\) son mejores que los obtenidos por el método mejorado de Euler con\(h=0.05\).
| Euler mejorado | Runge-Kutta | ||||
|---|---|---|---|---|---|
| x | h=0.1 | h=0.05 | h=0.1 | h-0.05 | “exacta” |
| 0.0 | 1.000000000 | 1.000000000 | 1.000000000 | 1.000000000 | 1.000000000 |
| 0.1 | 0.820040937 | 0.819050572 | 0.818753803 | 0.818751370 | 0.818751221 |
| 0.2 | 0.672734445 | 0.671086455 | 0.670592417 | 0.670588418 | 0.670588174 |
| 0.3 | 0.552597643 | 0.550543878 | 0.549928221 | 0.549923281 | 0.549922980 |
| 0.4 | 0.455160637 | 0.452890616 | 0.452210430 | 0.452205001 | 0.452204669 |
| 0.5 | 0.376681251 | 0.374335747 | 0.373633492 | 0.373627899 | 0.373627557 |
| 0.6 | 0.313970920 | 0.311652239 | 0.310958768 | 0.310953242 | 0.310952904 |
| 0.7 | 0.264287611 | 0.262067624 | 0.261404568 | 0.261399270 | 0.261398947 |
| 0.8 | 0.225267702 | 0.223194281 | 0.222575989 | 0.222571024 | 0.222570721 |
| 0.9 | 0.194879501 | 0.192981757 | 0.192416882 | 0.192412317 | 0.192412038 |
| 1.0 | 0.171388070 | 0.169680673 | 0.169173489 | 0.169169356 | 0.169169104 |
La Tabla 3.3.2 muestra resultados análogos para el problema del valor inicial no lineal
\[y'=-2y^2+xy+x^2,\ y(0)=1. \nonumber\]
Aplicamos el método mejorado de Euler a este problema en el Ejemplo 3.2.3.
| Euler mejorado | Runge-Kutta | ||||
|---|---|---|---|---|---|
| x | h=0.1 | h=0.05 | h=0.1 | h-0.05 | “exacta” |
| 0.0 | 1.000000000 | 1.000000000 | 1.000000000 | 1.000000000 | 1.000000000 |
| 0.1 | 0.840500000 | 0.838288371 | 0.837587192 | 0.837584759 | 0.837584494 |
| 0.2 | 0.733430846 | 0.730556677 | 0.729644487 | 0.729642155 | 0.729641890 |
| 0.3 | 0.661600806 | 0.658552190 | 0.657582449 | 0.657580598 | 0.657580377 |
| 0.4 | 0.615961841 | 0.612884493 | 0.611903380 | 0.611901969 | 0.611901791 |
| 0.5 | 0.591634742 | 0.588558952 | 0.587576716 | 0.587575635 | 0.587575491 |
| 0.6 | 0.586006935 | 0.582927224 | 0.581943210 | 0.581942342 | 0.581942225 |
| 0.7 | 0.597712120 | 0.594618012 | 0.593630403 | 0.593629627 | 0.593629526 |
| 0.8 | 0.626008824 | 0.622898279 | 0.621908378 | 0.621907553 | 0.621907458 |
| 0.9 | 0.670351225 | 0.667237617 | 0.666251988 | 0.666250942 | 0.666250842 |
| 1.0 | 0.730069610 | 0.726985837 | 0.726017378 | 0.726015908 | 0.726015790 |
Las tablas 3.3.3 y 3.3.4 muestran los resultados obtenidos aplicando los métodos semilineales Runge-Kutta y Runge-Kutta al problema de valor inicial
\[y'-2xy=1,\ y(0)=3, \nonumber\]
que consideramos en Ejemplos 3.3.3 y 3.3.4
| \(x\) | \(h=0.2\) | \(h=0.1\) | \(h=0.05\) | “Exacto” |
|---|---|---|---|---|
| \ (x\)” style="text-align:center; ">0.0 | \ (h=0.2\)” style="text-align:right; ">3.000000000 | \ (h=0.1\)” style="text-align:right; ">3.000000000 | \ (h=0.05\)” style="text-align:right; ">3.000000000 | 3.000000000 |
| \ (x\)” style="text-align:center; ">0.2 | \ (h=0.2\)” style="text-align:right; ">3.327846400 | \ (h=0.1\)” style="text-align:right; ">3.327851633 | \ (h=0.05\)” style="text-align:right; ">3.327851952 | 3.327851973 |
| \ (x\)” style="text-align:center; ">0.4 | \ (h=0.2\)” style="text-align:right; ">3.966044973 | \ (h=0.1\)” style="text-align:right; ">3.966058535 | \ (h=0.05\)” style="text-align:right; ">3.966059300 | 3.966059348 |
| \ (x\)” style="text-align:center; ">0.6 | \ (h=0.2\)” style="text-align:right; ">5.066996754 | \ (h=0.1\)” style="text-align:right; ">5.067037123 | \ (h=0.05\)” style="text-align:right; ">5.067039396 | 5.067039535 |
| \ (x\)” style="text-align:center; ">0.8 | \ (h=0.2\)” style="text-align:right; ">6.936534178 | \ (h=0.1\)” style="text-align:right; ">6.936690679 | \ (h=0.05\)” style="text-align:right; ">6.936700320 | 6.936700945 |
| \ (x\)” style="text-align:center; ">1.0 | \ (h=0.2\)” style="text-align:right; ">10.184232252 | \ (h=0.1\)” style="text-align:right; ">10.184877733 | \ (h=0.05\)” style="text-align:right; ">10.184920997 | 10.184923955 |
| \ (x\)” style="text-align:center; ">1.2 | \ (h=0.2\)” style="text-align:right; ">16.064344805 | \ (h=0.1\)” style="text-align:right; ">16.066915583 | \ (h=0.05\)” style="text-align:right; ">16.067098699 | 16.067111677 |
| \ (x\)” style="text-align:center; ">1.4 | \ (h=0.2\)” style="text-align:right; ">27.278771833 | \ (h=0.1\)” style="text-align:right; ">27.288605217 | \ (h=0.05\)” style="text-align:right; ">27.289338955 | 27.289392347 |
| \ (x\)” style="text-align:center; ">1.6 | \ (h=0.2\)” style="text-align:right; ">49.960553660 | \ (h=0.1\)” style="text-align:right; ">49.997313966 | \ (h=0.05\)” style="text-align:right; ">50.000165744 | 50.000377775 |
| \ (x\)” style="text-align:center; ">1.8 | \ (h=0.2\)” style="text-align:right; ">98.834337815 | \ (h=0.1\)” style="text-align:right; ">98.971146146 | \ (h=0.05\)” style="text-align:right; ">98.982136702 | 98.982969504 |
| \ (x\)” style="text-align:center; ">2.0 | \ (h=0.2\)” style="text-align:right; ">211.393800152 | \ (h=0.1\)” style="text-align:right; ">211.908445283 | \ (h=0.05\)” style="text-align:right; ">211.951167637 | 211.954462214 |
| \(x\) | \(h=0.2\) | \(h=0.1\) | \(h=0.05\) | “Exacto” |
|---|---|---|---|---|
| \ (x\) ">0.0 | \ (h=0.2\) ">3.000000000 | \ (h=0.1\) ">3.000000000 | \ (h=0.05\) ">3.000000000 | 3.000000000 |
| \ (x\) ">0.2 | \ (h=0.2\) ">3.327853286 | \ (h=0.1\) ">3.327852055 | \ (h=0.05\) ">3.327851978 | 3.327851973 |
| \ (x\) ">0.4 | \ (h=0.2\) ">3.966061755 | \ (h=0.1\) ">3.966059497 | \ (h=0.05\) ">3.966059357 | 3.966059348 |
| \ (x\) ">0.6 | \ (h=0.2\) ">5.067042602 | \ (h=0.1\) ">5.067039725 | \ (h=0.05\) ">5.067039547 | 5.067039535 |
| \ (x\) ">0.8 | \ (h=0.2\) ">6.936704019 | \ (h=0.1\) ">6.936701137 | \ (h=0.05\) ">6.936700957 | 6.936700945 |
| \ (x\) ">1.0 | \ (h=0.2\) ">10.184926171 | \ (h=0.1\) ">10.184924093 | \ (h=0.05\) ">10.184923963 | 10.184923955 |
| \ (x\) ">1.2 | \ (h=0.2\) ">16.067111961 | \ (h=0.1\) ">16.067111696 | \ (h=0.05\) ">16.067111678 | 16.067111677 |
| \ (x\) ">1.4 | \ (h=0.2\) ">27.289389418 | \ (h=0.1\) ">27.289392167 | \ (h=0.05\) ">27.289392335 | 27.289392347 |
| \ (x\) ">1.6 | \ (h=0.2\) ">50.000370152 | \ (h=0.1\) ">50.000377302 | \ (h=0.05\) ">50.000377745 | 50.000377775 |
| \ (x\) ">1.8 | \ (h=0.2\) ">98.982955511 | \ (h=0.1\) ">98.982968633 | \ (h=0.05\) ">98.982969450 | 98.982969504 |
| \ (x\) ">2.0 | \ (h=0.2\) ">211.954439983 | \ (h=0.1\) ">211.954460825 | \ (h=0.05\) ">211.954462127 | 211.954462214 |
El caso donde xno es el punto final izquierdo
En lo que va de este capítulo hemos considerado métodos numéricos para resolver un problema de valor inicial
\[\label{eq:3.3.3} y'=f(x,y),\quad y(x_0)=y_0\]
en un intervalo\([x_0,b]\), para el cual\(x_0\) es el punto final izquierdo. No hemos discutido métodos numéricos para resolver la Ecuación\ ref {eq:3.3.3} en un intervalo\([a,x_0]\), para lo cual\(x_0\) es el punto final correcto. Para ser específicos, ¿cómo podemos obtener valores aproximados\(y_{-1}\),\(y_{-2}\),...,\(y_{-n}\) de la solución de la Ecuación\ ref {eq:3.3.3} at\(x_0-h, \dots,x_0-nh\), ¿dónde\(h=(x_0-a)/n\)? Aquí está la respuesta a esta pregunta:
Considerar el problema de valor inicial
\[\label{eq:3.3.4} z'=-f(-x,z),\quad z(-x_0)=y_0,\]en el intervalo\([-x_0,-a]\), para el cual\(-x_0\) es el punto final izquierdo. Utilizar un método numérico para obtener valores aproximados\(z_1\)\(z_2\),,...,\(z_n\) de la solución de\(\eqref{eq:3.3.4}\) at\(-x_0+h\),\(-x_0+2h\),...,\(-x_0+nh=-a\). Entonces\(y_{-1}=z_1\),\(y_{-2}=z_2\),\(\dots\),\(y_{-n}=z_n\) son valores aproximados de la solución de\(\eqref{eq:3.3.3}\) at\(x_0-h\),\(x_0-2h\),...,\(x_0-nh=a\).
La justificación de esta respuesta se esboza en el Ejercicio 3.3.23. Observe lo fácil que es hacer el cambio el problema dado Ecuación\ ref {eq:3.3.3} al problema modificado Ecuación\ ref {eq:3.3.4}: primero reemplazar\(f\) por\(-f\) y luego reemplazar\(x\),\(x_0\), y\(y\) por\(-x\), y\(-x_0\), y\(z\), respectivamente.
Utilice el método Runge-Kutta con tamaño de paso\(h=0.1\) para encontrar valores aproximados de la solución de
\[\label{eq:3.3.5} (y-1)^2y'=2x+3,\quad y(1)=4\]
al\(x=0\),\(0.1\),\(0.2\),...,\(1\).
Solución
Primero reescribimos la Ecuación\ ref {eq:3.3.5} en la forma Ecuación\ ref {eq:3.3.3} como
\[\label{eq:3.3.6} y'={2x+3\over(y-1)^2},\quad y(1)=4.\]
Dado que la condición inicial\(y(1)=4\) se impone en el punto final correcto del intervalo\([0,1]\), aplicamos el método Runge-Kutta al problema del valor inicial
\[\label{eq:3.3.7} z'={2x-3\over(z-1)^2},\quad z(-1)=4\]
en el intervalo\([-1,0]\). (Debe verificar que la Ecuación\ ref {eq:3.3.7} está relacionada con la Ecuación\ ref {eq:3.3.6} ya que la Ecuación\ ref {eq:3.3.4} está relacionada con la Ecuación\ ref {eq:3.3.3}.) En el cuadro [tabla:3.3.5} se muestran los resultados. Al invertir el orden de las filas en la Tabla [tabla:3.3.5} y al cambiar los signos de los valores de se\(x\) obtienen las dos primeras columnas de la Tabla [tabla:3.3.6}. La última columna de la Tabla [tabla:3.3.6} muestra los valores exactos de\(y\), los cuales están dados por
\[y=1+(3x^2+9x+15)^{1/3}.\nonumber \]
(Dado que la ecuación diferencial en la Ecuación\ ref {eq:3.3.6} es separable, esta fórmula se puede obtener por el método de la Sección 2.2.)
| \(x\) | \(z\) |
|---|---|
| \ (x\) ">-1.0 | \ (z\) ">4.000000000 |
| \ (x\) ">-0.9 | \ (z\) ">3.944536474 |
| \ (x\) ">-0.8 | \ (z\) ">3.889298649 |
| \ (x\) ">-0.7 | \ (z\) ">3.834355648 |
| \ (x\) ">-0.6 | \ (z\) ">3.779786399 |
| \ (x\) ">-0.5 | \ (z\) ">3.725680888 |
| \ (x\) ">-0.4 | \ (z\) ">3.672141529 |
| \ (x\) ">-0.3 | \ (z\) ">3.619284615 |
| \ (x\) ">-0.2 | \ (z\) ">3.567241862 |
| \ (x\) ">-0.1 | \ (z\) ">3.516161955 |
| \ (x\) ">0.0 | \ (z\) ">3.466212070 |
| \(x\) | \(y\) | Exacto |
|---|---|---|
| \ (x\) ">0.00 | \ (y\) ">3.466212070 | 3.466212074 |
| \ (x\) ">0.10 | \ (y\) ">3.516161955 | 3.516161958 |
| \ (x\) ">0.20 | \ (y\) ">3.567241862 | 3.567241864 |
| \ (x\) ">0.30 | \ (y\) ">3.619284615 | 3.619284617 |
| \ (x\) ">0.40 | \ (y\) ">3.672141529 | 3.672141530 |
| \ (x\) ">0.50 | \ (y\) ">3.725680888 | 3.725680889 |
| \ (x\) ">0.60 | \ (y\) ">3.779786399 | 3.779786399 |
| \ (x\) ">0.70 | \ (y\) ">3.834355648 | 3.834355648 |
| \ (x\) ">0.80 | \ (y\) ">3.889298649 | 3.889298649 |
| \ (x\) ">0.90 | \ (y\) ">3.944536474 | 3.944536474 |
| \ (x\) ">1.00 | \ (y\) ">4.000000000 | 4.000000000 |
Te dejamos desarrollar un procedimiento para manejar la solución numérica de la Ecuación\ ref {eq:3.3.3} en un intervalo\([a,b]\) tal que\(a<x_0<b\) (Ejercicios 3.3.26 y 3.3.27).