Differential equations (contd)

In the first part we have dealt with differential equations of the form $$ y'(x) = f(x,y), $$ where $y(x_0) = y_0.$ We were given the function $f(x,y)$ and the numbers $x_0,y_0.$ And we started to solve it.

We had side-stepped the important question: does such an equation always possess a solution? And even if it does, is the solution unique?

The following theorem makes this problem meaningful.

Theorem Let $(x_0, y_0)\in{\mathbb R}^2.$ Let $\alpha,\beta>0$ be any two numbers. Let $f(x,y)$ be a continuous function defined over $[x_0-\alpha,x_0+\alpha]\times[y_0-\beta,y_0+\beta].$ We assume that $\frac{\partial f}{\partial y}$ is also continuous. We consider the differential equation $$\frac{dy}{dx} = f(x,y),\quad\quad\quad y(x_0)=y_0.$$ Then there exists $\delta>0$ such that this differential equation has unique solution $y(x)$ defined over $[x_0-\delta,x_0+\delta].$

We shall not prove this in this course.

Runge-Kutta method

In the first part, we have discussed Taylor's method (including its first order version, Euler's method). To get a reasonable approximation with Taylor's method in most real life problems, we need to use a high order (say 4-th order). But then we need to compute high order derivatives, which may be quite complicated. Runge-Kutta method aims to approximate Taylor's method without this complication. Instead, it needs to evaluate $f(x,y)$ at many points.

Just like Taylor's method, we have different Runge-Kutta methods for different orders. The general form is rather complicated. In this page we shall present only the 2-nd and 4-th order methods.

We start from the given point $(x_0,y_0).$ We are given a step size $h,$ and we want to approximate $y(x)$ at the points $x_1,...,x_n,$ where $x_i = x_0+h\,i.$ Thus, the output of the Runge-Kutta methods consists of $(x_i,y_i)$ for $i=0,...,n.$

The 2-nd order Runge-Kutta method is given by $$\begin{eqnarray*} k_1 & = & h\,f(x_i,y_i)\\ k_2 & = & h\,f\left(x_i+\frac{h}{2},y_i+\frac{k_1}{2}\right)\\ y_{i+1} & = & y_i + k_2. \end{eqnarray*}$$ This method has the same order of precision as 2-nd order Taylor method. To see this use the 2-variable Taylor expansion: $$ f(x+\delta x, y+\delta y) = f(x,y)+\delta x f_1(x,y)+ \delta y f_2(x,y) +\mbox{ 2-nd order terms}. $$ where $f_1,f_2$ are partial derivatives of $f$ w.r.t. its first and second arguments, respectively.

We shall apply this to $f\left(x_i+\frac{h}{2},y_i+\frac{k_1}{2}\right).$ To simplify notation, we shall write $f,$ $f_1$ etc to mean $f(x_i,y_i),$ $f_1(x_i, y_i)$ etc, respectively. Then we see that $$\begin{eqnarray*} k_2 & =& h\left[f + \frac{h}{2} f_1 + \frac{k_1}{2} f_2 + h^2\mbox{-terms}\right]\\ & =& hf + \frac{h^2}{2} f_1 + \frac{h\,k_1}{2} f_2 + h^3\mbox{-terms}\\ & =& hf + \frac{h^2}{2} f_1 + \frac{h^2}{2} f f_2 + h^3\mbox{-terms}, \end{eqnarray*}$$ since $k_1 = hf.$

So the second order Runge-Kutta method is $$ y_{i+1} = y_i + hf + \frac{h^2}{2} f_1 + \frac{h^2}{2} f f_2 + h^3\mbox{-terms} $$ Let's compare this with the 2-nd order Taylor's method: $$\begin{eqnarray*} y_{i+1} = y_i + y_i'h + \frac{y_i''}{2}h^2 & =& y_i + hf + \frac{y_i''}{2}h^2. \end{eqnarray*}$$ Now by the multivariate chain rule, $$ y_i'' = f_1 + f_2y_i' = f_1+ff_2. $$ So the 2nd order Taylor's method becomes $$ y_{i+1} = y_i + hf + \frac{h^2}{2}(f_1+ff_2), $$ which coincides with the 2nd order Runge-Kutta method (ignoring $h^3$-terms).

The most popular Runge-Kutta method is the 4-th order version: $$\begin{eqnarray*} k_1 & = & hf(x_i,y_i)\\ k_2 & = & hf\left(x_i+\frac{h}{2},y_i+\frac{k_1}{2}\right)\\ k_3 & = & hf\left(x_i+\frac{h}{2},y_i+\frac{k_2}{2}\right)\\ k_4 & = & hf(x_i+h,y_i+k_3)\\ y_{i+1} & = & y_i + \frac{1}{6}(k_1 + 2k_2 + 2k_3 + k_4). \end{eqnarray*}$$ This has the same order of precision as the fourth order Taylor method. The proof of this fact follows the same argument as above, but is rather messy.

Here is a simple code for performing a single RK4 step:

rk4 = function(f, x,y, h) {
  k1 = h * f(x,y)
  k2 = h*f(x+h/2, y+k1/2)
  k3 = h * f(x=h/2, y+k2/2)
  k4 = h*f(x+h, y+k3)
  y + (k1+2*k2+2*k3+k4)/6
}

We may use it like this:

f = function(x,y) sin(x+y)
x = seq(0,1,len=10)
y = numeric(10); y[1] = 1

for(i in 2:length(x)) { 
   y[i] = rk4(f,x[i-1],y[i-1],x[i]-x[i-1])
}

Applications

Simple pendulum (revisited)

In the first part, we had already arrived at the following system for a simple pendulum: $$\begin{eqnarray*} \theta ' & = & \omega ~~(\equiv f_1(t,\theta,\omega))\\ \omega' & = & -\frac{g}{L} \theta ~~(\equiv f_2(t,\theta,\omega)). \end{eqnarray*}$$ The 4-th order Runge-Kutta method for this becomes $$\begin{eqnarray*} t_{i+1} & = t_i + h\\ \theta_{i+1} & = \theta_i + \frac{1}{6}(k_1+2k_2+2k_3+k_4)\\ \omega_{i+1} & = \omega_i + \frac{1}{6}(j_1+2j_2+2j_3+j_4), \end{eqnarray*}$$ where $$\begin{eqnarray*} k_1 & = & hf_1(t_i,\theta_i,\omega_i) = h\omega_i\\ j_1 & = & hf_2(t_i,\theta_i,\omega_i) = -\frac{gh}{L} \theta_i\\ k_2 & = & hf_1(t_i+\tfrac{h}{2},\theta_i+\tfrac{k_1}{2},\omega_i+\tfrac{j_1}{2}) = h[\omega_i+\tfrac{j_1}{2}]\\ j_2 & = & hf_2(t_i+\tfrac{h}{2},\theta_i+\tfrac{k_1}{2},\omega_i+\tfrac{j_1}{2}) = -\tfrac{gh}{L} [\theta_i+\tfrac{k_1}{2}]\\ k_3 & = & hf_1(t_i+\tfrac{h}{2},\theta_i+\tfrac{k_2}{2},\omega_i+\tfrac{j_2}{2}) = h[\omega_i+\tfrac{j_2}{2}]\\ j_3 & = & hf_2(t_i+\tfrac{h}{2},\theta_i+\tfrac{k_2}{2},\omega_i+\tfrac{j_2}{2}) = -\frac{gh}{L} [\theta_i+\tfrac{k_2}{2}]\\ k_4 & = & hf_1(t_i+h,\theta_i+k_3,\omega_i+j_3)=h[\omega_i+j_3]\\ j_4 & = & hf_2(t_i+h,\theta_i+k_3,\omega_i+j_3)=-\frac{gh}{L}[\theta_i+k_3]. \end{eqnarray*}$$

EXERCISE: Also try with the extensible rod case discussed earlier.

Gravity well

PROJECT:

Simulate the gravity well for the first part using 4-th order Runge-Kutta method.

Comments

To post an anonymous comment, click on the "Name" field. This will bring up an option saying "I'd rather post as a guest."

Table of contents