8.03: Runge-Kutta 2nd-Order Method for Solving Ordinary Differential Equations
- Page ID
- 126430
\( \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}}} \)
\(\newcommand{\avec}{\mathbf a}\) \(\newcommand{\bvec}{\mathbf b}\) \(\newcommand{\cvec}{\mathbf c}\) \(\newcommand{\dvec}{\mathbf d}\) \(\newcommand{\dtil}{\widetilde{\mathbf d}}\) \(\newcommand{\evec}{\mathbf e}\) \(\newcommand{\fvec}{\mathbf f}\) \(\newcommand{\nvec}{\mathbf n}\) \(\newcommand{\pvec}{\mathbf p}\) \(\newcommand{\qvec}{\mathbf q}\) \(\newcommand{\svec}{\mathbf s}\) \(\newcommand{\tvec}{\mathbf t}\) \(\newcommand{\uvec}{\mathbf u}\) \(\newcommand{\vvec}{\mathbf v}\) \(\newcommand{\wvec}{\mathbf w}\) \(\newcommand{\xvec}{\mathbf x}\) \(\newcommand{\yvec}{\mathbf y}\) \(\newcommand{\zvec}{\mathbf z}\) \(\newcommand{\rvec}{\mathbf r}\) \(\newcommand{\mvec}{\mathbf m}\) \(\newcommand{\zerovec}{\mathbf 0}\) \(\newcommand{\onevec}{\mathbf 1}\) \(\newcommand{\real}{\mathbb R}\) \(\newcommand{\twovec}[2]{\left[\begin{array}{r}#1 \\ #2 \end{array}\right]}\) \(\newcommand{\ctwovec}[2]{\left[\begin{array}{c}#1 \\ #2 \end{array}\right]}\) \(\newcommand{\threevec}[3]{\left[\begin{array}{r}#1 \\ #2 \\ #3 \end{array}\right]}\) \(\newcommand{\cthreevec}[3]{\left[\begin{array}{c}#1 \\ #2 \\ #3 \end{array}\right]}\) \(\newcommand{\fourvec}[4]{\left[\begin{array}{r}#1 \\ #2 \\ #3 \\ #4 \end{array}\right]}\) \(\newcommand{\cfourvec}[4]{\left[\begin{array}{c}#1 \\ #2 \\ #3 \\ #4 \end{array}\right]}\) \(\newcommand{\fivevec}[5]{\left[\begin{array}{r}#1 \\ #2 \\ #3 \\ #4 \\ #5 \\ \end{array}\right]}\) \(\newcommand{\cfivevec}[5]{\left[\begin{array}{c}#1 \\ #2 \\ #3 \\ #4 \\ #5 \\ \end{array}\right]}\) \(\newcommand{\mattwo}[4]{\left[\begin{array}{rr}#1 \amp #2 \\ #3 \amp #4 \\ \end{array}\right]}\) \(\newcommand{\laspan}[1]{\text{Span}\{#1\}}\) \(\newcommand{\bcal}{\cal B}\) \(\newcommand{\ccal}{\cal C}\) \(\newcommand{\scal}{\cal S}\) \(\newcommand{\wcal}{\cal W}\) \(\newcommand{\ecal}{\cal E}\) \(\newcommand{\coords}[2]{\left\{#1\right\}_{#2}}\) \(\newcommand{\gray}[1]{\color{gray}{#1}}\) \(\newcommand{\lgray}[1]{\color{lightgray}{#1}}\) \(\newcommand{\rank}{\operatorname{rank}}\) \(\newcommand{\row}{\text{Row}}\) \(\newcommand{\col}{\text{Col}}\) \(\renewcommand{\row}{\text{Row}}\) \(\newcommand{\nul}{\text{Nul}}\) \(\newcommand{\var}{\text{Var}}\) \(\newcommand{\corr}{\text{corr}}\) \(\newcommand{\len}[1]{\left|#1\right|}\) \(\newcommand{\bbar}{\overline{\bvec}}\) \(\newcommand{\bhat}{\widehat{\bvec}}\) \(\newcommand{\bperp}{\bvec^\perp}\) \(\newcommand{\xhat}{\widehat{\xvec}}\) \(\newcommand{\vhat}{\widehat{\vvec}}\) \(\newcommand{\uhat}{\widehat{\uvec}}\) \(\newcommand{\what}{\widehat{\wvec}}\) \(\newcommand{\Sighat}{\widehat{\Sigma}}\) \(\newcommand{\lt}{<}\) \(\newcommand{\gt}{>}\) \(\newcommand{\amp}{&}\) \(\definecolor{fillinmathshade}{gray}{0.9}\)Lesson 1: Theory of Runge-Kutta 2nd-Order Method
Learning Objectives
After successful completion of this lesson, you should be able to:
1) list the formulas of the Runge-Kutta 2nd order method for ordinary differential equations and know how to use them.
What is the Runge-Kutta 2nd order method?
The Runge-Kutta 2nd order method is a numerical technique used to solve an ordinary differential equation of the form
\[\frac{dy}{dx} = f\left( x,y \right),\ y\left( x_0 \right) = y_{0}\;\;\;\;\;\;\;\;\;\;\;\; (\PageIndex{1.1}) \nonumber\]
Only first-order ordinary differential equations of the form of Equation \((\PageIndex{1.1})\) can be solved by using the Runge-Kutta 2nd order method. In other sections, we discuss how the Euler and Runge-Kutta methods are used to solve higher-order ordinary or coupled (simultaneous) ordinary differential equations.
How does one write a first-order ordinary differential equation in the above form?
Runge-Kutta 2nd order method
Euler’s method is given by
\[y_{i + 1} = y_{i} + f\left( x_{i},y_{i} \right)h\;\;\;\;\;\;\;\;\;\;\;\; (\PageIndex{1.2}) \nonumber\]
where
\[y(x_0)=y_0 \nonumber\]
\[h = x_{i + 1} - x_{i} \nonumber\]
To understand the Runge-Kutta 2nd order method, we need to first derive Euler’s method from the Taylor series.
\[y_{i + 1} = y_{i} \ + \ \left. \ \frac{dy}{dx} \right|_{x_{i},y_{i}}\left( x_{i + 1} - x_{i} \right) \ + \ \frac{1}{2!}\left. \ \frac{d^{2}y}{dx^{2}} \right|_{x_{i},y_{i}}\left( x_{i + 1} - x_{i} \right)^{2} \ + \ \frac{1}{3!}\left. \ \frac{d^{3}y}{dx^{3}} \right|_{x_{i},y_{i}}\left( x_{i + 1} - x_{i} \right)^{3} \ + ... \nonumber\]
Since \(\displaystyle \frac{dy}{dx}=f(x,y)\),
\[y_{i+1} = y_{i} \ + \ f(x_{i},y_{i})\left( x_{i + 1} - x_{i} \right) \ + \ \frac{1}{2!}f'(x_{i},y_{i})\left( x_{i + 1} - x_{i} \right)^{2} \ + \ \frac{1}{3!}f''(x_{i},y_{i})\left( x_{i + 1} - x_{i} \right)^{3} \ + ...\;\;\;\;\;\;\;\;\;\;\;\; (\PageIndex{1.3}) \nonumber\]
As you can see, the first two terms of the Taylor series in Equation \((\PageIndex{1.3})\)
\[y_{i + 1} = y_{i} + f\left( x_{i},y_{i} \right)h \nonumber\]
represent Euler’s method and hence can be considered to be a Runge-Kutta 1st order method.
The true error in the approximation is given by
\[E_{t} = \frac{f^{\prime}\left( x_{i},y_{i} \right)}{2!}h^{2} + \frac{f^{\prime\prime}\left( x_{i},y_{i} \right)}{3!}h^{3} + ...\;\;\;\;\;\;\;\;\;\;\;\; (\PageIndex{1.4}) \nonumber\]
So what would a 2nd order method formula look like? It would include one more term of the Taylor series as follows.
\[y_{i + 1} = y_{i} + f\left( x_{i},y_{i} \right)h + \frac{1}{2!}f^{\prime}\left( x_{i},y_{i} \right)h^{2}\;\;\;\;\;\;\;\;\;\;\;\; (\PageIndex{1.5}) \nonumber\]
Let us take a generic example of a first-order ordinary differential equation
\[\frac{dy}{dx} = e^{- 2x} - 3y, \ y\left( 0 \right) = 5 \nonumber\]
\[f\left( x,y \right) = e^{- 2x} - 3y \nonumber\]
Now since \(y\) is a function of \(x\),
\[\begin{split} f^{\prime}\left( x,y \right) &= \frac{\partial f\left( x,y \right)}{\partial x} + \frac{\partial f\left( x,y \right)}{\partial y}\frac{dy}{dx}\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; (\PageIndex{1.6})\\ &= \frac{\partial}{\partial x}\left( e^{- 2x} - 3y \right) + \frac{\partial}{\partial y}\left\lbrack \left( e^{- 2x} - 3y \right) \right\rbrack\left( e^{- 2x} - 3y \right)\\ &= - 2e^{- 2x} + ( - 3)\left( e^{- 2x} - 3y \right)\\ &= - 5e^{- 2x} + 9y \end{split}\]
The 2nd order formula for the above example would be
\[\begin{split} y_{i + 1} &= y_{i} + f\left( x_{i},y_{i} \right)h + \frac{1}{2!}f^{\prime}\left( x_{i},y_{i} \right)h^{2}\\ &= y_{i} + \left( e^{- 2x_{i}} - 3y_{i} \right)h + \frac{1}{2!}\left( - 5e^{- 2x_{i}} + 9y_{i} \right)h^{2} \end{split}\]
However, we already recognize the pitfall of having to find \(f^{\prime}\left( x,y \right)\) in the above method. What Runge and Kutta did was write the 2nd order method as
\[y_{i + 1} = y_{i} + \left( a_{1}k_{1} + a_{2}k_{2} \right)h\;\;\;\;\;\;\;\;\;\;\;\; (\PageIndex{1.7}) \nonumber\]
where
\[k_{1} = f\left( x_{i},y_{i} \right)\;\;\;\;\;\;\;\;\;\;\;\; (\PageIndex{1.8}\text{a}) \nonumber\]
\[k_{2} = f\left( x_{i} + p_{1}h,y_{i} + q_{11}k_{1}h \right)\;\;\;\;\;\;\;\;\;\;\;\; (\PageIndex{1.8}\text{b}) \nonumber\]
This form allows one to take advantage of the 2nd order method without having to calculate \(f^{\prime}\left( x,y \right)\).
So how do we find the unknowns \(a_{1}\), \(a_{2}\), \(p_{1},\) and \(q_{11}\)? Without proof, equating Equation \((\PageIndex{1.5})\) and \((\PageIndex{1.7})\), gives three equations.
\[a_{1} + a_{2} = 1\;\;\;\;\;\;\;\;\;\;\;\; (\PageIndex{1.9}\text{a}) \nonumber\]
\[a_{2}p_{1} = \frac{1}{2}\;\;\;\;\;\;\;\;\;\;\;\; (\PageIndex{1.9}\text{b}) \nonumber\]
\[a_{2}q_{11} = \frac{1}{2}\;\;\;\;\;\;\;\;\;\;\;\; (\PageIndex{1.9}\text{c}) \nonumber\]
Since we have 3 equations and 4 unknowns, we can assume the value of one of the unknowns. The other three are then determined from the three equations. Generally, the value of \(a_{2}\) is chosen to evaluate the other three constants. The three values used for \(a_{2}\) are \(\displaystyle \frac{1}{2}\), \(1\) and \(\displaystyle \frac{2}{3}\), and are known as Heun’s Method, the midpoint method, and Ralston’s method, respectively.
Heun’s Method
Here \(a_{2} = \displaystyle \frac{1}{2}\) is chosen, and from Equations \((\PageIndex{1.9} \text{a-c})\),
\[a_{1}\ = \frac{1}{2} \nonumber\]
\[p_{1}\ = 1 \nonumber\]
\[q_{11} = 1 \nonumber\]
resulting in
\[y_{i + 1} = y_{i} + \left( \frac{1}{2}k_{1} + \frac{1}{2}k_{2} \right)h\;\;\;\;\;\;\;\;\;\;\;\; (\PageIndex{1.10}) \nonumber\]
where
\[k_{1} = f\left( x_{i},y_{i} \right)\;\;\;\;\;\;\;\;\;\;\;\; (\PageIndex{1.11}\text{a}) \nonumber\]
\[k_{2} = f\left( x_{i} + h,y_{i} + k_{1}h \right)\;\;\;\;\;\;\;\;\;\;\;\; (\PageIndex{1.11}\text{b}) \nonumber\]
This method is graphically explained in Figure \(\PageIndex{1.1}\).
Midpoint Method
Here \(a_{2} = 1\) is chosen, and from Equations \((\PageIndex{1.9} \text{a-c})\),
\[a_{1} = 0 \nonumber\]
\[p_{1} = \frac{1}{2} \nonumber\]
\[q_{11} = \frac{1}{2} \nonumber\]
resulting in
\[y_{i + 1} = y_{i} + k_{2}h\;\;\;\;\;\;\;\;\;\;\;\; (\PageIndex{1.12}) \nonumber\]
where
\[k_{1} = f\left( x_{i},y_{i} \right)\;\;\;\;\;\;\;\;\;\;\;\; (\PageIndex{1.13}\text{a}) \nonumber\]
\[k_{2} = f\left( x_{i} + \frac{1}{2}h,y_{i} + \frac{1}{2}k_{1}h \right)\;\;\;\;\;\;\;\;\;\;\;\; (\PageIndex{1.13}\text{b}) \nonumber\]
Ralston’s Method
Here \(a_{2} = \displaystyle\frac{2}{3}\) is chosen, and from Equations \((\PageIndex{1.9} \text{a-c})\),
\[a_{1} = \frac{1}{3} \nonumber\]
\[p_{1} = \frac{3}{4} \nonumber\]
\[q_{11} = \frac{3}{4} \nonumber\]
resulting in
\[y_{i + 1} = y_{i} + \left( \frac{1}{3}k_{1} + \frac{2}{3}k_{2} \right)h\;\;\;\;\;\;\;\;\;\;\;\;(\PageIndex{1.14}) \nonumber\]
where
\[k_{1} = f\left( x_{i},y_{i} \right)\;\;\;\;\;\;\;\;\;\;\;\; (\PageIndex{1.15}\text{a}) \nonumber\]
\[k_{2} = f\left( x_{i} + \frac{3}{4}h,y_{i} + \frac{3}{4}k_{1}h \right)\;\;\;\;\;\;\;\;\;\;\;\; (\PageIndex{1.15}\text{b}) \nonumber\]
Lesson 2: Application of Runge-Kutta 2nd-Order Method
Learning Objectives
After successful completion of this lesson, you should be able to:
1) Use Runge-Kutta 2nd order method to solve first-order ordinary differential equations.
Recap
In the previous lesson, we discussed the theory behind the Runge-Kutta 2nd order method of solving an ordinary differential equation of the form \(\dfrac{dy}{dx}=f(x,y)\) where \(y(x_0)=y_0\). In this lesson, we take an example of how to apply the algorithm of the Runge-Kutta 2nd order method.
A ball at \(1200\ \text{K}\) is allowed to cool down in air at an ambient temperature of \(300\ \text{K}\). Assuming heat is lost only due to radiation, the differential equation for the temperature of the ball is given by
\[\frac{d \theta}{dt} = - 2.2067 \times 10^{-12}\ \left(\theta^{4} - 81 \times 10^{8}\right) \nonumber\]
where \(\theta\) is in \(\text{K}\) and \(t\) in seconds. Find the temperature at \(t = 480\) seconds using Runge-Kutta 2nd order method. Assume a step size of \(h = 240\) seconds.
Solution
\[\frac{d \theta}{dt} = - 2.2067 \times 10^{- 12}\left( \theta^{4} - 81 \times 10^{8} \right) \nonumber\]
\[f\left( t,\theta \right) = - 2.2067 \times 10^{- 12}\left( \theta^{4} - 81 \times 10^{8} \right) \nonumber\]
As per Heun’s method given in the previous lesson for an ordinary differential equation,
\[\frac{d \theta}{dt} = f(t,\theta) \nonumber\]
Heun’s method formula is given by
\[\theta_{i + 1} = \theta_{i} + \left( \frac{1}{2}k_{1} + \frac{1}{2}k_{2} \right)h \nonumber\]
\[k_{1} = f\left( t_{i},\theta_{i} \right) \nonumber\]
\[k_{2} = f\left( t_{i} + h,\theta_{i} + k_{1}h \right) \nonumber\]
For Step 1,
\[i = 0,\ t_{0} = 0,\ \theta_{0} = \theta(0) = 1200\ \text{K} \nonumber\]
\[\begin{split} t_1&=t_{0}+h\\ &=0+240\\ &=240\ \text{s}\end{split}\]
\[\begin{split} k_{1} &= f\left( t_{0},\theta_{o} \right)\\ &= f\left( 0,1200 \right)\\ &= - 2.2067 \times 10^{- 12}\left( 1200^{4} - 81 \times 10^{8} \right)\\ &= - 4.5579 \end{split}\]
\[\begin{split} k_{2} &= f\left( t_{0} + h,\theta_{0} + k_{1}h \right)\\ &= f\left( 0 + 240,1200 + \left( - 4.5579 \right)240 \right)\\ &= f\left( 240,106.09 \right)\\ &= - 2.2067 \times 10^{- 12}\left( 106.09^{4} - 81 \times 10^{8} \right)\\ &= 0.017595 \end{split}\]
\[\begin{split} \theta_{1} &= \theta_{0} + \left( \frac{1}{2}k_{1} + \frac{1}{2}k_{2} \right)h\\ &= 1200 + \left( \frac{1}{2}\left( - 4.5579 \right) + \frac{1}{2}\left( 0.017595 \right) \right)240\\ &= 1200 + \left( - 2.2702 \right)240\\ &= 655.16\ \text{K}\\ &\approx \theta(240) \end{split}\]
For Step 2
\[i = 1,t_{1} = 240 \ \text{s},\theta_{1} = 655.16\ \text{K} \nonumber\]
\[\begin{split} &=t_{1}+h\\ &=240+240\\ &=480\ \text{s}\end{split}\]
\[\begin{split} k_{1} &= f\left( t_{1},\theta_{1} \right)\\ &= f\left( 240,655.16 \right)\\ &= - 2.2067 \times 10^{- 12}\left( 655.16^{4} - 81 \times 10^{8} \right)\\ &= - 0.38869 \end{split}\]
\[\begin{split} k_{2} &= f\left( t_{1} + h,\theta_{1} + k_{1}h \right)\\ &= f\left( 240 + 240,655.16 + \left( - 0.38869 \right)240 \right)\\ &= f\left( 480,561.87 \right)\\ &= - 2.2067 \times 10^{- 12}\left( 561.87^{4} - 81 \times 10^{8} \right)\\ &= - 0.20206 \end{split}\]
\[\begin{split} \theta_{2} &= \theta_{1} + \left( \frac{1}{2}k_{1} + \frac{1}{2}k_{2} \right)h\\ &= 655.16 + \left( \frac{1}{2}\left( - 0.38869 \right) + \frac{1}{2}\left( - 0.20206 \right) \right)240\\ &= 655.16 + \left( - 0.29538 \right)240\\ &= 584.27\ \text{K}\\ &\approx \theta(480) \end{split}\]
The results from Heun’s method are compared with the exact results in Figure \(\PageIndex{2.1}\).
The exact solution of the ordinary differential equation is given by the solution of a nonlinear equation as
\[0.92593\ln\frac{\theta - 300}{\theta + 300} - 1.8519\tan^{- 1}\left( 0.0033333\theta \right) = - 0.22067 \times 10^{- 3}t - 2.9282 \nonumber\]
The solution to this nonlinear equation at \(t = 480\ \text{s}\) is
\[\theta(480) = 647.57\ \text{K} \nonumber\]
Using a smaller step size would increase the accuracy of the result, as given in Table \(\PageIndex{2.1}\) and Figure \(\PageIndex{2.2}\) below.
\(\text{Step size}, \ h\) | \(\theta\left ( 480 \right)\) | \(E_{t}\) | \(\left| \epsilon_ {t} \right|\%\) |
---|---|---|---|
\(480\) | \(-393.87\) | \(1041.4\) | \(160.82\) |
\(240\) | \(584.27\) | \(63.304\) | \(9.7756\) |
\(120\) | \(651.35\) | \(-3.7762\) | \(0.58313\) |
\(60\) | \(649.91\) | \(-2.3406\) | \(0.36145\) |
\(30\) | \(648.21\) | \(-0.63219\) | \(0.097625\) |
In Table \(\PageIndex{2.2}\), the Euler’s method and the Runge-Kutta 2nd order method results are shown as a function of step size,
\(\text{Step size}, \ h\) | \(\theta(480)\) | |||
---|---|---|---|---|
\(\text{Euler}\) | \(\text{Heun}\) | \(\text{Midpoint}\) | \(\text{Ralston}\) | |
\(480\) | \(-987.84\) | \(-393.87\) | \(1208.4\) | \(449.78\) |
\(240\) | \(110.32\) | \(584.27\) | \(976.87\) | \(690.01\) |
\(120\) | \(546.77\) | \(651.35\) | \(690.20\) | \(667.71\) |
\(60\) | \(614.97\) | \(649.91\) | \(654.85\) | \(652.25\) |
\(30\) | \(632.77\) | \(648.21\) | \(649.02\) | \(648.61\) |
while in Figure \(\PageIndex{2.3}\), the comparison is shown over time.
How do these three methods compare with results obtained if we found \(\mathbf{f^\prime (x,y)}\) directly?
Of course, we know that since we are including the first three terms in the series, if the solution is a polynomial of order two or less (that is, quadratic, linear, or constant), any of the three methods are exact. But for any other case, the results are going to be different.
Let us take the example of showing the difference in results.
\[\frac{dy}{dx} = e^{- 2x} - 3y,y\left( 0 \right) = 5. \nonumber\]
If we directly find \(f^{\prime}\left( x,y \right)\), the first three terms of the Taylor series give
\[y_{i + 1} = y_{i} + f\left( x_{i},y_{i} \right)h + \frac{1}{2!}f^{\prime}\left( x_{i},y_{i} \right)h^{2} \nonumber\]
where
\[f\left( x,y \right) = e^{- 2x} - 3y \nonumber\]
\[f^{\prime}\left( x,y \right) = - 5e^{- 2x} + 9y \nonumber\]
For a step size of \(h = 0.2\), using Heun’s method, we find
\[y\left( 0.6 \right) = 1.0930 \nonumber\]
The exact solution
\[y\left( x \right) = e^{- 2x} + 4e^{- 3x} \nonumber\]
gives
\[\begin{split} y\left( 0.6 \right) &= e^{- 2\left( 0.6 \right)} + 4e^{- 3\left( 0.6 \right)}\\ &= 0.96239 \end{split}\]
Then the absolute relative true error is
\[\begin{split} \left| \epsilon_{t} \right| &= \left| \frac{0.96239 - 1.0930}{0.96239} \right| \times 100\\ &= 13.571\% \end{split}\]
For the same problem, the results from Euler’s method and the three Runge-Kutta methods are given in Table \(\PageIndex{2.3}\).
\(Exact\) | \(Euler\) | \(Direct\ 2nd\) | \(Heun\) | \(Midpoint\) | \(Ralston\) | |
---|---|---|---|---|---|---|
\(y(0.6)\) | \(0.96239\) | \(0.4955\) | \(1.0930\) | \(1.1012\) | \(1.0974\) | \(1.0994\) |
\(\left| \epsilon_{t} \right|\) % | \(48.514\) | \(13.571\) | \(14.423\) | \(14.029\) | \(14.236\) |
Lesson 3: Derivation of Runge-Kutta 2nd-Order Method
Learning Objectives
After successful completion of this lesson, you should be able to:
1) derive the Runge-Kutta 2nd order method for the first-order ordinary differential equation.
How do we get the 2nd order Runge-Kutta method equations?
We wrote the 2nd order Runge-Kutta equations without proof to solve
\[\frac{dy}{dx} = f\left( x,y \right),\ y\left( 0 \right) = y_{0}\;\;\;\;\;\;\;\;\;\;\;\; (\PageIndex{3.1}) \nonumber\]
as
\[y_{i + 1} = y_{i} + \left( a_{1}k_{1} + a_{2}k_{2} \right)h \;\;\;\;\;\;\;\;\;\;\;\; (\PageIndex{3.2}) \nonumber\]
where
\[k_{1} = f\left( x_{i},y_{i} \right)\;\;\;\;\;\;\;\;\;\;\;\; (\PageIndex{3.3}\text{a}) \nonumber\]
\[k_{2} = f\left( x_{i} + p_{1}h,y_{i} + q_{11}k_{1}h \right)\;\;\;\;\;\;\;\;\;\;\;\; (\PageIndex{3.3}\text{b}) \nonumber\]
and
\[a_{1} + a_{2} = 1\;\;\;\;\;\;\;\;\;\;\;\; (\PageIndex{3.4}\text{a}) \nonumber\]
\[a_{2}p_{2} = \frac{1}{2}\;\;\;\;\;\;\;\;\;\;\;\; (\PageIndex{3.4}\text{b}) \nonumber\]
\[a_{2}q_{11} = \frac{1}{2}\;\;\;\;\;\;\;\;\;\;\;\; (\PageIndex{3.4}\text{c}) \nonumber\]
The advantage of using 2nd order Runge-Kutta method equations is based on not having to find the derivative of \(f\left( x,y \right)\) symbolically in the ordinary differential equation
So how do we get the above three Equations \((\PageIndex{3.4}\text{a-c})\)? The question is answered below.
Writing out the first three terms of the Taylor series are
\[y_{i + 1} = y_{i} + \left. \ \frac{dy}{dx} \right|_{x_{i}y_{i}}h + \left. \ \frac{1}{2!}\frac{d^{2}y}{dx^{2}} \right|_{x_{i}y_{i}}h^{2} + O\left( h^{3} \right)\;\;\;\;\;\;\;\;\;\;\;\; (\PageIndex{3.5}) \nonumber\]
where
\[h = x_{i + 1} - x_{i} \nonumber\]
Since
\[\frac{dy}{dx} = f\left( x,y \right) \nonumber\]
we can rewrite the Taylor series as
\[y_{i + 1} = y_{i} + f\left( x_{i},y_{i} \right)h + \frac{1}{2!}f^{\prime}\left( x_{i},y_{i} \right)h^{2} + O\left( h^{3} \right)\;\;\;\;\;\;\;\;\;\;\;\; (\PageIndex{3.6}) \nonumber\]
Now
\[f^{\prime}\left( x,y \right) = \frac{\partial f\left( x,y \right)}{\partial x} + \frac{\partial f\left( x,y \right)}{\partial y}\frac{dy}{dx}.\;\;\;\;\;\;\;\;\;\;\;\; (\PageIndex{3.7}) \nonumber\]
Hence
\[\begin{split} y_{i + 1} &= y_{i} + f\left( x_{i},y_{i} \right)h + \frac{1}{2!}\left( \left. \ \frac{\partial f}{\partial x} \right|_{x_{i},y_{i}} + \left. \ \frac{\partial f}{\partial y} \right|_{x_{i},y_{i}} \times \left. \ \frac{dy}{dx} \right|_{x_{i},y_{i}} \right)h^{2} + O\left( h^{3} \right)\\ &= y_{i} + f\left( x_{i},y_{i} \right)h + \frac{1}{2}\left. \ \frac{\partial f}{\partial x} \right|_{x_{i},y_{i}}h^{2} + \frac{1}{2}\left. \ \frac{\partial f}{\partial y} \right|_{x_{i},y_{i}}f\left( x_{i},y_{i} \right)h^{2} + O\left( h^{3} \right)\;\;\;\;\;\;\;\;\;\;\;\; (\PageIndex{3.8}) \end{split}\]
Now the term used in the Runge-Kutta 2nd order method for \(k_{2}\) can be written as a Taylor series of two variables with the first three terms as
\[\begin{split} k_{2} &= f\left( x_{i} + p_{1}h,y_{i} + q_{11}k_{1}h \right)\\ &= f\left( x_{i},y_{i} \right) + p_{1}h\left. \ \frac{\partial f}{\partial x} \right|_{x_{i},y_{i}} + q_{11}k_{1}h\left. \ \frac{\partial f}{\partial y} \right|_{x_{i},y_{i}} + O\left( h^{2} \right)\;\;\;\;\;\;\;\;\;\;\;\; (\PageIndex{3.9}) \end{split}\]
Hence
\[\begin{split} y_{i + 1} &= y_{i} + \left( a_{1}k_{1} + a_{2}k_{2} \right)h\\ &= y_{i} + \left( a_{1}f\left( x_{i},y_{i} \right) + a_{2}\left\{ f\left( x_{i},y_{i} \right) + p_{1}h\left. \ \frac{\partial f}{\partial x} \right|_{x_{i},y_{i}} + q_{11}k_{1}h\left. \ \frac{\partial f}{\partial y} \right|_{x_{i},y_{i}} + O\left( h^{2} \right) \right\} \right)h\\ &= y_{i} + \left( a_{1} + a_{2} \right)hf\left( x_{i},y_{i} \right) + a_{2}p_{1}h^{2}\left. \ \frac{\partial f}{\partial x} \right|_{x_{i},y_{i}} + a_{2}q_{11}f\left( x_{i},y_{i} \right)h^{2}\left. \ \frac{\partial f}{\partial y} \right|_{x_{i},y_{i}} + O\left( h^{3} \right)\;\;\;\;\;\;\;\;\;\;\;\; (\PageIndex{3.10}) \end{split}\]
Equating the terms in Equation \((\PageIndex{3.8})\) and Equation \((\PageIndex{3.10})\), we get
\[a_{1} + a_{2} = 1 \nonumber\]
\[a_{2}p_{1} = \frac{1}{2} \nonumber\]
\[a_{2}q_{11} = \frac{1}{2} \nonumber\]
Multiple Choice Test
(1). To solve the ordinary differential equation
\[3\frac{dy}{dx} + xy^{2} = \sin x,\ y\left( 0 \right) = 5 \nonumber\]
by the Runge-Kutta 2nd order method, you need to rewrite the equation as
(A) \(\displaystyle \frac{dy}{dx} = \sin x - xy^{2},\ y\left( 0 \right) = 5\)
(B) \(\displaystyle \frac{dy}{dx} = \frac{1}{3}\left( \sin x - xy^{2} \right),\ y\left( 0 \right) = 5\)
(C) \(\displaystyle \frac{dy}{dx} = \frac{1}{3}\left( - \cos x - \frac{xy^{3}}{3} \right),\ y\left( 0 \right) = 5\)
(D) \(\displaystyle \frac{dy}{dx} = \frac{1}{3}\sin x,\ y\left( 0 \right) = 5\)
(2). Given
\[3\frac{dy}{dx} + 5y^{2} = \sin x,\ y\left( 0.3 \right) = 5 \nonumber\]
and using a step size of \(h = 0.3\), the value of \(y\left( 0.9 \right)\) using the Runge-Kutta 2nd order Heun’s method is most nearly
(A) \(-4297.4\)
(B) \(-4936.7\)
(C) \(-0.21336\times 10^{14}\)
(D) \(-0.24489\times 10^{14}\)
(3). Given
\[3\frac{dy}{dx} + 5\sqrt{y} = e^{0.1x},\ y\left( 0.3 \right) = 5 \nonumber\]
and using a step size of \(h = 0.3\), the best estimate of \(\displaystyle \frac{dy}{dx}\left( 0.9 \right)\) using the Runge-Kutta 2nd order midpoint method most nearly is
(A) \(-2.2473\)
(B) \(-2.2543\)
(C) \(-2.6188\)
(D) \(-3.2045\)
(4). The velocity \((\text{m}/\text{s})\) of a body is given as a function of time (seconds) by
\[v\left( t \right) = 200\ln\left( 1 + t \right) - t,\ t \geq 0 \nonumber\]
Using the Runge-Kutta 2nd order Ralston method with a step size of \(5\) seconds, the distance in meters traveled by the body from \(t = 2\) to \(t = 12\) seconds is estimated most nearly as
(A) \(3904.9\)
(B) \(3939.7\)
(C) \(6556.3\)
(D) \(39397\)
(5). The Runge-Kutta 2nd order method can be derived by using the first three terms of the Taylor series of writing the value of \(y_{i + 1}\) (that is, the value of \(y\) at \(x_{i + 1}\)) in terms of \(y_{i}\) (that is, the value of \(y\) at \(x_{i}\)) and all the derivatives of \(y\) at \(x_{i}\). If \(h = x_{i + 1} - x_{i}\), the explicit expression for \(y_{i + 1}\) if the first three terms of the Taylor series are chosen for solving the ordinary differential equation
\[\frac{dy}{dx} + 5y = 3e^{- 2x}, \ y\left( 0 \right) = 7 \nonumber\]
would be
(A) \(\displaystyle y_{i + 1} = y_{i} + \left( 3e^{- 2x_{i}} - 5y_{i} \right)h + 5\frac{h^{2}}{2}\)
(B) \(\displaystyle y_{i + 1} = y_{i} + \left( 3e^{- 2x_{i}} - 5y_{i} \right)h + \left( - 21e^{- 2x_{i}} + 25y_{i} \right)\frac{h^{2}}{2}\)
(C) \(\displaystyle y_{i + 1} = y_{i} + \left( 3e^{- 2x_{i}} - 5y_{i} \right)h + \left( - 6e^{- 2x_{i}} \right)\frac{h^{2}}{2}\)
(D) \(\displaystyle y_{i + 1} = y_{i} + \left( 3e^{- 2x_{i}} - 5y_{i} \right)h + \left( - 6e^{- 2x_{i}} + 5 \right)\frac{h^{2}}{2}\)
(6). A spherical ball is taken out of a furnace at \(1200\ \text{K}\) and is allowed to cool in air. You are given the following
\[\text{radius of ball} = 2\ \text{cm} \nonumber\]
\[\text{specific heat of ball} = 420\ \text{J/kg} \cdot \text{K} \nonumber\]
\[\text{density of ball} = 7800\ \text{kg/m}^{3} \nonumber\]
\[\text{convection coefficient} = 350\ \text{J/s} \cdot \text{m}^{2} \cdot \text{K} \nonumber\]
\[\text{ambient temperature} = 300 \ \text{K} \nonumber\]
The ordinary differential equation that is given for the temperature \(\theta\) of the ball is
\[\frac{d \theta}{dt} = - 2.20673 \times 10^{- 13}\left( \theta^{4} - 81 \times 10^{8} \right) \nonumber\]
if only radiation is accounted for. If convection is accounted for in addition to radiation, the ordinary differential equation is
(A) \(\displaystyle \frac{d \theta}{dt} = - 2.20673 \times 10^{- 13}\left( \theta^{4} - 81 \times 10^{8} \right) - 1.6026 \times 10^{- 2}\left( \theta - 300 \right), \theta(0)=1200\ \text{K}\)
(B) \(\displaystyle \frac{d \theta}{dt} = - 2.20673 \times 10^{- 13}\left( \theta^{4} - 81 \times 10^{8} \right) - 4.3982 \times 10^{- 2}\left( \theta - 300 \right), \theta(0)=1200\ \text{K}\)
(C) \(\displaystyle \frac{d \theta}{dt} = - 1.6026 \times 10^{- 2}\left( \theta - 300 \right), \theta(0)=1200\ \text{K}\)
(D) \(\displaystyle \frac{d \theta}{dt} = - 4.3982 \times 10^{- 2}\left( \theta - 300 \right), \theta(0)=1200\ \text{K}\)
For complete solution, go to
http://nm.mathforcollege.com/mcquizzes/08ode/quiz_08ode_runge2nd_solution.pdf
Problem Set
(1). For the ordinary differential equation
\[2\frac{dy}{dx} + 3 xy = e^{- 1.5x},y(0) = 5 \nonumber\]
use a step size of \(h = 2\) and Runge-Kutta 2nd order Ralston’s method to find \(y(6)\).
- Answer
-
\(y(6) \approx -23684\)
(2). For the ordinary differential equation
\[2\frac{dy}{dx} + 3 y = e^{- 1.5x},y(0) = 5 \nonumber\]
a) use Runge-Kutta 2nd order midpoint method to find \(y(2.5)\) using \(h = 1.25\),
b) use Runge-Kutta 2nd order midpoint method to find \(\displaystyle \frac{dy}{dx}\ (2.5)\) using \(h = 1.25\),
c) true value of \(y(2.5)\),
d) absolute relative true error for part (a),
e) true error of \(\displaystyle \frac{dy}{dx}\ (2.5),\)
f) absolute relative true error for part (b).
- Answer
-
\(a)\ y\left(2.5\right) \approx 3.5433\)
\(b)\ -5.3031\)
\(c)\ 0.14699\)
\(d)\ 2310.6\%\)
\(e)\ -0.20872\)
\(f)\ 2440.8\%\)
(3). Find the approximate value of the integral \(\displaystyle \int_{3}^{8}{2e^{0.6x} \ dx}\) using Runge-Kutta 2nd order Heun’s method.
a) Use \(h = 2.5\) and compare the value with the exact value.
b) Use \(h = 1.25\) and compare the value with the exact value.
- Answer
-
\(a)\ 454.46,\ 18.083\%\)
\(b)\ 402.74,\ 4.6441\%\)
(4). From problem (3a), is the value obtained using Heun’s Runge-Kutta 2nd order method same as 2-segment LRAM (Left Endpoint Rectangular Approximation), or 2-segment MRAM (Midpoint Rectangular Approximation), or 2-segment RRAM (Right Endpoint Rectangular Approximation) method of integration or composite trapezoidal rule with 2 segments? Explain.
- Answer
-
\(454.46\), MRAM
(5). A water tank with a hole at the top and a circular hole at the bottom is shown in the figure.
Given
\[\text{Diameter of the tank,}\ D = 6\ \text{m} \nonumber\]
\[\text{Initial height of the water,}\ h = 5\ \text{m} \nonumber\]
\[\text{Diameter of the hole,}\ d = 4.95\ \text{cm} \nonumber\]
The flow rate \(\dot{Q}\) of water through the bottom hole is given by
\[\dot{Q} = -Av, \nonumber\]
where
\[A =\text{cross-sectional area of the hole,} \nonumber\]
\[v =\text{velocity of the water flowing through the hole.} \nonumber\]
Since
\[v = \sqrt{2gh} \nonumber\]
where
\[h =\text{height of the water from the bottom of the tank,} \nonumber\]
\[g = \text{acceleration due to gravity.} \nonumber\]
From this, we get
\[\dot{Q} = - A\sqrt{2gh} \nonumber\]
Also, since the volume \(V\) of the water at height \(h\) is given by
\[V = \frac{1}{3}\pi h^{2}\left( \frac{3D}{2} - h \right) \nonumber\]
where
\[\frac{dV}{dt} = \left( \pi hD - \pi h^{2} \right)\frac{dh}{dt} \nonumber\]
Since
\[\dot{Q} = \frac{dV}{dt} \nonumber\]
\[\left( \pi hD - \pi h^{2} \right)\frac{dh}{dt} = - A\sqrt{2gh} \nonumber\]
\[\left( \pi hD - \pi h^{2} \right)\frac{dh}{dt} = - \frac{\pi d^{2}}{4}\sqrt{2gh} \nonumber\]
\[\frac{dh}{dt} = - \frac{d^{2}\sqrt{2gh}}{4\left( hD - h^{2} \right)} \nonumber\]
a) Use Runge-Kutta 2nd order Heun’s method with a step size of 10 minutes to find the height of water at \(t = 50\) minutes.
b) Based on the results from part (a), estimate the time when the height of the water is \(4\) meters.
c) Compare the results of part (a) and part (b) with the exact values.
- Answer
-
\(a)\ 2.9282\ \text{m}\)
\(b)\ 19.248\ \text{min}\)
\(c)\) part (a): \(h=2.9371\ \text{m},\ 0.303\%\); part (b): \(t=19.41\ \text{min},\ 0.835\%\)
(6). For a spherical ball taken out of furnace at 2500 K that loses heat due to radiation and convection, the following is given.
\[\text{Radius of ball,}\ r = 1.0\ \text{cm}, \nonumber\]
\[\text{Density of the ball,}\ \rho = 3000\ \text{kg/m}^{3}, \nonumber\]
\[\text{Specific heat,}\ C = 1000\ \text{J}/(\text{kg} \cdot \text{K}), \nonumber\]
\[\text{Emittance,}\ \epsilon = 0.5, \nonumber\]
\[\text{Stefan-Boltzmann Constant,}\ \sigma = 5.67 \times 10^{- 8}\ \text{J}/(\text{s} \cdot \text{m}^{2} \cdot \text{K}^{4}), \nonumber\]
\[\text{Convective cooling coefficient,}\ h = 500\ \text{J/}(\text{s} \cdot \text{m}^{2} \cdot \text{K}), \nonumber\]
\[\text{Initial temperature of the ball,}\ \theta(0) = 2500\ \text{K}, \nonumber\]
\[\text{Ambient temperature,}\ \theta_{a} = 300\ \text{K}, \nonumber\]
the ordinary differential equation that governs the temperature of the ball is given by
\[\frac{d\theta}{dt} = - 2.8349 \times 10^{- 12}(\theta^{4} - 81 \times 10^{8}) - 0.05(\theta - 300). \nonumber\]
Hint: The above formula was calculated by using the knowledge that
\[\text{Rate of heat lost due to radiation} = A \epsilon \sigma(\theta^{4} - \theta_{a}^{4}) \nonumber\]
\[\text{Rate of heat lost due to convection} = hA(\theta - \theta_{a}) \nonumber\]
\[\text{Heat stored} = mc\theta \nonumber\]
where
\[A =\text{surface area of the ball,} \nonumber\]
\[m =\text{mass of the ball.} \nonumber\]
a) What are following at \(t = 0\ \text{s}\)?
- temperature,
- rate of change of temperature,
- rate of heat loss due to radiation,
- rate of heat loss due to convection,
- rate of heat stored,
b) Find the following at \(t = 10\ \text{s}\) using Runge Kutta 2nd order Ralston’s method and a step size of \(h = 2.5\).
- temperature,
- rate of change of temperature,
- rate of heat loss due to radiation,
- rate of heat loss due to convection,
- rate of heat stored,
c) Can either convection or radiation be neglected if we are interested in finding the temperature profile between \(0\) and \(10\ \text{s}\)? Base your answer on quantitative reasoning.
d) Find the time when the temperature will be \(2000\ \text{K}\) using Ralston’s method and a step size of \(h = 2.5\ \text{s}\).
- Answer
-
\(a)\ 2500\ \text{K},\ -220.72\ \text{K/s},\ 1391.3\ \text{W},\ 1382.3\ \text{W},\ -2773.6\ \text{W}\)
\(b)\ 1381.9\ \text{K},\ -64.41\ \text{K/s},\ 129.63\ \text{W},\ 679.78\ \text{W},\ -809.41\ \text{W}\)
\(c)\) Neither can be neglected, as rate of heat lost due to convection and radiation is of same order.
\(d)\ 3.1591\ \text{s}\)