Subsections

• Runge-Kutta Methods
  • Euler's Method
  • Improvement of Euler's Method
  • Runge-Kutta Method
  • Systems of Equations
  • Adaptive Runge-Kutta Method
  
  
• Stiffness and Multistep Methods
  • Stiffness
  • Multistep Methods
  
  
• Boundary-Value and Engenvalue Problems
  • General Methods of Boundary-Value Problems
  • ODEs and Eigenvalues with Libraries and Packages
  
  
• Engineering Applications: Ordinary Differential Equations

Ordinary Differential Equations

Differential equation is

• differential equations : composed of an unknown function and its derivatives.
• rate equations : it expresses the rate of change of a variable as a function of variables and parameters.

Variables are divied by

• dependent variables
• independent variables

Differential equation is classified as

• ordinary differential equation : one independent variable
• partial differential equation : more than one independent variables

A differential equation is usually accompanied by auxiliary conditions to specify the solution completely. For first-order ODEs an initial value is required to determine the constant and obtain a unique solution.

• initial-value problem : all conditions are specified at the same value of the independent variable.
• boundary-value problem : specification of conditions occurs at different values of the independent variable.

Runge-Kutta Methods

Ordinary differential equation is

$$\frac{dy}{dx} = f(x,y)$$

The solution is

$$\textrm{New value} = \textrm{old value} + \textrm{slope} \times \textrm{step size}$$

or, in mathmatical terms,

$$y_{i+1} = y_i + \phi \times h$$ (7.1)

The slope estimate of $\phi$ is used to extrapolate from an old value $y_i$ to a new value $y_{i+1}$ over a distance $h$.

Figure 7.1: Graphical depiction of a one-step method.

Euler's Method

Euler's method is

$$y_{i+1} = y_i + f(x_i,y_i) h$$ (7.2)

More smaller step-size gives more accurate solution but further step-size reduction requires much computation times.

Improvement of Euler's Method

A funcdermental source of error in Euler's method is that the derivative at the beginning of the interval is assumed to apply across the entire interval.

Heun's Method : average two derivatives which are obtained at the initial point and the end point. The slope at the beginning of an interval

$$y'_i = f(x_i,y_i)$$ (7.3)

is used to extrapolate linearly to $y_{i+1}$

$$y^0_{i+1} = y_i + f(x_i,y_i)h$$ (7.4)

This equation is called a predictor equation. The slope at the end of the interval

$$y'_{i+1} = f(x_{i+1},y_{i+1})$$ (7.5)

Thus, the two slopes can be combined to obtain an average slope for the interval

$$\bar{y}'=\frac{y_i+y'_{i+1}}{2} = \frac{f(x_i,y_i)+f(x_{i+1},y^0_{i+1})}{2}$$ (7.6)

This average slope is then used to extrapolate linearly from $y_i$ to $y_{i+1}$

$$y_{i+1} = y_i + \frac{f(x_i,y_i)+f(x_{i+1},y^0_{i+1})}{2}h$$ (7.7)

which is called a correct equation. See the figure 25.9 at p 688.

The Heun method is a predictor-corrector approach.

The Midpoint Method : use Euler's method to predict a value of $y$ at the midpoint of the interval.

$$y_{i+1/2} = y_i + f(x_i,y_i)\frac{h}{2}$$ (7.8)

This slope is then used to extrapolate linear form from $x_i$ to $x_{i+1}$

$$y_{i+1} = y_i + f(x_{i+1/2},y_{i+1/2})\frac{h}{2}$$ (7.9)

Runge-Kutta Method

Runge-Kutta methods achieve the accuracy of a Taylor series approach without requiring the calculation of higher derivatives.

$$y_{i+1} = y_i + \phi(x_i, y_i, h) h$$ (7.10)

where $\phi(x_i, y_i, h)$ is called an increment function, which can be interpreted as a representative slope over the interval. The increment function is

$$\phi = a_1 k_1 + a_2 k_2 + \cdots + a_n k_n$$ (7.11)

where the $a$'s are constants and the $k$'s are

Notice that the $k$'s are recurrence relationship. Because each $k$ is a functional evaluation, this recurrence makes RK methods efficient for computer calculations.

Second-order Runge-Kutta Methods

The second-order version of RK method:

$$y_{i+1} = y_i + (a_1 k_1 + a_2 k_2) h$$ (7.13)

where

To determine values for the constant $a_1$, $a_2$, $p_1$, and $q_{11}$, use a Taylor's series for $y_{i+1}$ in terms of $y_i$ and $f(x_i,y_i)$

$$y_{i+1} = y_i + f(x_i,y_i) h + \frac{f'(x_i,y_i)}{2!}h^2$$ (7.15)

where $f'(x_i,y_i)$ is

$$f'(x_i,y_i) = \frac{\partial f(x,y)}{\partial x} + \frac{\partial f(x,y)}{\partial y} \frac{dy}{dx}$$ (7.16)

Then,

$$y_{i+1} = y_i + f(x_i,y_i) h + \frac{\partial f}{\partial x} + \frac{\partial f}{\partial y} \frac{dy}{dx}\frac{h^2}{2!}$$ (7.17)

The basic strategy underlying Runge-Kutta methods is to use algebraic manipulations to solve for values of $a_1$, $a_2$, $p_1$, and $q_{11}$ that make eq (7.13) and eq (7.17) equivalent.

The Taylor's series for a two-variable function is

$$f(x_i + p_1 h, y_i + q_{11}k_1 h) = f(x_i,y_i) + p_1h\frac{\partial f}{\partial x} + q_{11}k_1h \frac{\partial f}{\partial y} +O(h^2)$$ (7.18)

This gives

$$y_{i+1} = y_i + a_1 h f(x_i,y_i) + a_2 h f(x_i,y_i) + a_2p_1h^2 \... ...}{\partial x} + a_2 q_{11}h^2 f(x_i,y_i) \frac{\partial f}{\partial y} + O(h^3)$$ (7.19)

By correcting terms

$$y_{i+1} = y_i + \left[a_1 f(x_i,y_i) + a_2 f(x_i,y_i)\right]h + \... ...ial x} + a_2 q_{11} f(x_i,y_i) \frac{\partial f}{\partial y} \right]h^2+ O(h^3)$$ (7.20)

Comparing this equation with eq (7.17)

$$a_1 + a_2 = 1$$ (7.21)

$$a_1 p_1 = \frac{1}{2}$$ (7.22)

$$a_2 q_{11} = \frac{1}{2}$$ (7.23)

Because these three equations contain the four unknown constants, we must assume a value of one of the unknowns to determine the other three. Suppose that we specify a value for $a_2$.

$$a_1 = 1 - a_2$$ (7.24)

$$p_1 = q_{11} = \frac{1}{2a_2}$$ (7.25)

Also we can choose an infinite number of values for $a_2$, there are an infinite number of second-order RK methods.

Heun Method with a Single Corrector($a_2 = 1/2$)

Assume $a_2 = 1/2$

$$a_1 = 1/2, \quad p_1 = q_{11} = 1$$ (7.26)

These parameters yield

$$y_{i+1} = y_i + \left( \frac{1}{2}k_1 + \frac{1}{2}k_2 \right) h$$ (7.27)

where

$$k_1 = f(x_i,y_i)$$ (7.28)

$$k_2 = f(x_i + h, y_i + k_1h)$$ (7.29)

Midpoint Method($a_2 = 1$)

$$y_{i+1} = y_i + k_2 h$$ (7.30)

where

$$k_1 = f(x_i,y_i)$$ (7.31)

$$k_2 = f\left(x_i + \frac{1}{2}h, y_i + \frac{1}{2}k_1h \right)$$ (7.32)

Ralston's Method ($a_2 = 2/3$)

$$y_{i+1} = y_i + \left( \frac{1}{3}k_1 + \frac{2}{3}k_2 \right) h$$ (7.33)

where

$$k_1 = f(x_i,y_i)$$ (7.34)

$$k_2 = f\left(x_i + \frac{3}{4}h, y_i + \frac{3}{4}k_1h \right)$$ (7.35)

See the figure 25.14 at p 699.

Fourth-order Runge-Kutta Methods

The classical fourth-order RK method

$$y_{i+1} = y_i + \frac{1}{6}(k_1 + 2k_2 + 2k_3 + k_4) h$$ (7.36)

where

See the figure 25.16 at p 704.

Systems of Equations

The procedure for solving a system of equations simply involves applying the one-step technique for every equation at each step before proceeding to the next step.

Adaptive Runge-Kutta Method

• Automatically adjust the step size to avoid overkills
• Need an estimate of the local truncation error be obtained at each step

Strategies for adaptive RK

• Use different step sizes to calculate local error
• Use different order of RK method to calculate local error

Stiffness and Multistep Methods

Stiffness

A stiffness system is one involving repidly changing components together with slowly changing ones. In many cases, the rapidly varying components die away quickly, after which the solution becomes dominated by the slowly varying components.

An example of a single stiff ODE is

$$\frac{dy}{dt} = -1000y + 3000 - 2000e^{-t}$$ (7.38)

If $y(0) = 0$,

$$y = 3 - 0.998 e^{-1000t} - 2.002 e^{-t}$$ (7.39)

The solution is initially dominated by the fast exponential term ( $e^{-1000t}$). After a very short period $t < 0.005$, this transient dies out and the solution becomes dictated by the slow exponential ($e^{-t}$).

Step size consideration: Insight into the step size required for stability for a solution.

$$\frac{dy}{dt} = -ay$$ (7.40)

If $y(0) = y_0$,

$$y = y_0 e^{-at}$$ (7.41)

Thus, the solution starts at $y_0$ and asymptotically approaches zero.

Use Euler's method

$$y_{i+1} = y_i + \frac{dy_i}{dt}h$$ (7.42)

Substituting $dy/dt$

$$y_{i+1} = y_i - ay_i h$$ (7.43)

or

$$y_{i+1} = (1- ah)y_i$$ (7.44)

The stability of this formula clearly depend on the step size $h$. That is, $\vert 1-ah\vert$ must be less than 1.

For the fast transient part, the step size to maintain stability must be very small. In addition, an even smaller step size is required to obtain an accurate solution.

Implicit method: developed by evaluating the derivative at the future time

$$y_{i+1} = y_i + \frac{dy_{i+1}}{dt}h$$ (7.45)

Substituting the derivative term

$$y_{i+1} = y_i + ay_{i+1}h$$ (7.46)

which can be solved for

$$y_{i+1} = \frac{y_i}{1+ah}$$ (7.47)

For this case, regardless of the size of the step, $\vert y_i\vert \rightarrow 0$ as $i \rightarrow \infty$. See the figure 26.2 at p 722.

Multistep Methods

The one-step methods utilize information at a single point $x_i$ to predict a value of the dependent variable $y_{i+1}$ at a future point $x_{i+1}$. Alternative approaches, called multistep methods, are base on information of the previous points. The curvature of the lines connecting these previous values provides information regarding the trajectory of the solution. The multistep methods exploit this information to solve ODEs.

The non-self-starting Huen method

The Heun method use Euler's method as predictor

$$y^0_{i+1} = y_i + f(x_i,y_i)h$$ (7.48)

and the trapezoidal rule as a corrector

$$y_{i+1} = y_i + \frac{f(x_i,y_i)+f(x_{i+1},y^0_{i+1})}{2}h$$ (7.49)

Thus, the predictor and the corrector have local truncation errors of $O(h^2)$ and $O(h^3)$, repectively. Consequently, one way to improve Heun's method is to develop a predictor that has a local error of $O(h^3)$. This can be accomplished by using Euler's method and the slope at $y_i$, and extra information from a previous point $y_{i-1}$, as in

$$y^0_{i+1} = y_{i-1} + f(x_i,y_i)2h$$ (7.50)

Notice that eq. (7.50) attains $O(h^3)$ at the expense of employing a larger step size, $2h$. In addition, eq. (7.50) is not self-starting because it involves a previous value of the dependent variable $y_{i-1}$. Because of the fact it is called the non-self-starting Heun method.

Derivation of non-self-starting Heun method

Consider the general ODE

$$\frac{dy}{dx} = f(x,y)$$ (7.51)

Integrating between limits at $i$ and $i+1$

$$\int^{y_{i+1}}_{y_i} dy = \int^{x_{i+1}}_{x_i} f(x,y) dx$$ (7.52)

Integrated value is

$$y_{i+1} = y_i + \int^{x_{i+1}}_{x_i} f(x,y) dx$$ (7.53)

Use trapezoidal rule to integrate the second term at right hand side

$$y_{i+1} = y_i + \frac{f(x_i,y_i) + f(x_{i+1},y_{i+1})}{2}h$$ (7.54)

which is the corrector equation for the Heun method and the trapezoidal rule gives the local truncation error of $O(h^3)$.

A similar approach can be used to derive the predictor. For this case, the integration limits are from $i-1$ to $i+1$.

$$\int^{y_{i+1}}_{y_{i-1}} dy = \int^{x_{i+1}}_{x_{i-1}} f(x,y) dx$$ (7.55)

which can be integrated and rearranged to yield

$$y_{i+1} = y_{i-1} + \int^{x_{i+1}}_{x_{i-1}} f(x,y) dx$$ (7.56)

Use the first Newton-Cote open integration formula

$$\int^{x_{i+1}}_{x_{i-1}} f(x,y) dx = 2hf(x_i,y_i)$$ (7.57)

which is called the midpoint method.

$$y_{i+1} = y_{i-1} + 2hf(x_i,y_i)$$ (7.58)

which is the predictor for the non-self-starting Heun.

Integration formulas

The non-self-starting Heun method employs an open integration formula (the midpoint method) to make an initial estimate. This predictor step requires a previous data point. Then, a closed integration formula (the trapezoidal rule) is applied iteratively to improve the solution.

Newton-Cotes formulas estimate the integral over an interval spanning several points. In constrast, the Adam formulas use a set of points from an interval to estimate the integral solely for the last segment in the interval. See the figure 26.7 p 734.

• Newton-Cotes Formulas
  • Open formulas

    $$y_{i+1} = y_{i-n} + \int^{x_{i+1}}_{x_{i-n}} f_n(x) dx$$ (7.59)

    
    
    if $n=1$

    $$y_{i+1} = y_{i-1} + 2hf_i$$ (7.60)

    
    
    which is referred to as the midpoint method and was used previously as the predictor in the non-self-starting Heun method.
  • Closed formulas

    $$y_{i+1} = y_{i-n+1} + \int^{x_{i+1}}_{x_{i-n+1}} f_n(x) dx$$ (7.61)

    
    
    if $n=1$

    $$y_{i+1} = y_i + \frac{h}{2}(f_i + f_{i+1})$$ (7.62)

    
    
    which is equivalent to the trapezoidal rule.
• Adams Formulas : Many popular computer algorithms for multistep solution of ODEs are based on these methods.
  • Open formulas(Adams-Bashforth) : start with a forward Taylor series expansion at $x_i$

    $$y_{i+1} = y_i + f_i h + \frac{f'_i}{2}h^2 + \frac{f''_i}{6} h^3 + \cdots$$ (7.63)

    
    
    which can also be written as

    $$y_{i+1} = y_i + h\left( f_i + \frac{h}{2}f'_i + \frac{h^2}{3!}f''_i + \cdots \right)$$ (7.64)

    
    
    Use a backward difference

    $$f'_i = \frac{f_i - f_{i-1}}{h} + \frac{f''_i}{2}h + O(h^2)$$ (7.65)

    
    
    Then

    $$y_{i+1} = y_i + h\left(\frac{3}{2}f_i - \frac{1}{2}f_{i-1} \right) + \frac{5}{12}h^3 f''_i + O(h^4)$$ (7.66)

    
    

  • Closed formulas(Adams-Moulton) : start with a backward Taylor series around $x_{i+1}$

    $$y_i = y_{i+1} - f_{i+1} h + \frac{f'_{i+1}}{2}h^2 + \cdots$$ (7.67)

    
    
    Solving for $y_{i+1}$ yields

    $$y_{i+1} = y_i + h \left( f_{i+1} - \frac{h}{2}f'_{i+1} + \cdot \right)$$ (7.68)

    
    
    Use a difference to approximate the first derivative

    $$f'_{i+1} = \frac{f_{i+1} - f_i}{h} + \frac{f''_{i+1}}{2}h +h(O^2)$$ (7.69)

    
    
    Then

    $$y_{i+1} = y_i + h \left(\frac{1}{2}f_{i+1} + \frac{1}{2}f_i \right) - \frac{1}{12}h^3f''_{i+1} - O(h^4)$$ (7.70)

    
    

Boundary-Value and Engenvalue Problems

Boundary-value : which is specified at the extreme points or boundaries of a system.

Classification of boundary condition

• Dirichlet condition : the value of independent variable is specified at a boundary
• Neumann condition : the value of the derivative of independent variable is specified at a boundary

General Methods of Boundary-Value Problems

The shooting method: based on converting the boundary-value problem into an equivalent initial-value problem. A trial-and-error approach is then implemented to solve the initial-value version.

Finite-difference methods: finite divided differences are substituted for the derivatives in the original equation. Thus, a linear differential equation is transformed into a set of simultaneous algebraic equations.

ODEs and Eigenvalues with Libraries and Packages

• Matlab
  • ode23
  • ode45
• IMSL
  • IVPRK
  • IVPAG
  • BVPFD
  • BVPMS

Engineering Applications: Ordinary Differential Equations

See the textbook