[<< wikibooks] Introduction to Numerical Methods/Ordinary Differential Equations
Resources:

notes on Ordinary Differential Equations
notes on Euler's Method
notes on Runge-Kutta 2nd Order Method
notes on Runge-Kutta 4th Order Method

A differential equation is a equation that relates some function with its derivatives. Such relationship is commonly found in dynamic system modeling where the rate of change of a quantity is affect by another quantity, which can be a function of other quantities. Therefore differential equations are often used to mathematically represent such relations in many disciplines including engineering, physics, economics, and biology. We will focus on ordinary differential equations, which involves no partial derivative.
The solution to a differential equation is the function or a set of functions that satisfies the equation. Some simple differential equations with explicit formulas are solvable analytically, but we can always use numerical methods to estimate the answer using computers to a certain degree of accuracy.
http://martinfowler.com/bliki/TwoHardThings.html
"Two hard things in computer science: cache invalidation, naming things, and off-by-one errors."


= Initial Value Problem =
The general form of first order differential equation: 
  
    
      
        
          y
          ′
        
        =
        f
        (
        x
        ,
        y
        )
      
    
    {\displaystyle y'=f(x,y)}
   or 
  
    
      
        
          y
          ′
        
        =
        f
        (
        x
        ,
        y
        (
        x
        )
        )
      
    
    {\displaystyle y'=f(x,y(x))}
  . The solution contains an arbitrary constant (the constant of integration) because 
  
    
      
        z
        (
        x
        )
        =
        y
        (
        x
        )
        +
        c
      
    
    {\displaystyle z(x)=y(x)+c}
   satisfies all differential equations 
  
    
      
        y
        (
        x
        )
      
    
    {\displaystyle y(x)}
   satisfies. With an auxiliary condition, e.g. y(a)=b, we can solve y'=F(x, y). This is known as the initial value problem.
An indefinite integral of a function 
  
    
      
        f
        (
        x
        )
      
    
    {\displaystyle f(x)}
   is a class of functions 
  
    
      
        F
        (
        x
        )
      
    
    {\displaystyle F(x)}
   so that 
  
    
      
        
          F
          ′
        
        (
        x
        )
        =
        f
        (
        x
        )
      
    
    {\displaystyle F'(x)=f(x)}
  . The indefinite integral is also known as the antiderivative. The definite integral of a function over a fixed interval is a number.


= Connection with Integration =
To calculate the integral 
  
    
      
        F
        (
        x
        )
        =
        
          ∫
          
            a
          
          
            x
          
        
        f
        (
        t
        )
        d
        t
      
    
    {\displaystyle F(x)=\int _{a}^{x}f(t)dt}
   is to find a 
  
    
      
        F
        (
        x
        )
      
    
    {\displaystyle F(x)}
   that satisfies the condition 
  
    
      
        
          
            
              d
              F
              (
              x
              )
            
            
              d
              x
            
          
        
        =
        f
        (
        x
        )
      
    
    {\displaystyle {\frac {dF(x)}{dx}}=f(x)}
  ,  
  
    
      
        F
        (
        a
        )
        =
        0
      
    
    {\displaystyle F(a)=0}
  . In other words,  integration can be used to solve differential equations with a initial condition.


= Euler's Method =

Euler's method uses a step-wise approach to solve ordinary differential equations with an initial condition.
Given 
  
    
      
        
          
            
              d
              y
            
            
              d
              x
            
          
        
        =
        f
        (
        x
        ,
        y
        )
      
    
    {\displaystyle {\frac {dy}{dx}}=f(x,y)}
   and 
  
    
      
        y
        (
        
          x
          
            0
          
        
        )
        =
        
          y
          
            0
          
        
      
    
    {\displaystyle y(x_{0})=y_{0}}
   (
  
    
      
        x
      
    
    {\displaystyle x}
   is an independent variable and 
  
    
      
        y
      
    
    {\displaystyle y}
   is a dependent variable) Euler's methods solves for 
  
    
      
        
          y
          
            i
          
        
      
    
    {\displaystyle y_{i}}
   for 
  
    
      
        x
        =
        
          x
          
            i
          
        
      
    
    {\displaystyle x=x_{i}}
  .
Euler's method can be derived in a number of ways. Lets look at a geographical description.
Given that 
  
    
      
        
          
            
              d
              y
            
            
              d
              x
            
          
        
        =
        f
        (
        x
        ,
        y
        )
      
    
    {\displaystyle {\frac {dy}{dx}}=f(x,y)}
   at a particular point 
  
    
      
        (
        
          x
          
            i
          
        
        ,
        
          y
          
            i
          
        
        )
      
    
    {\displaystyle (x_{i},y_{i})}
   we can approximate 
  
    
      
        f
        (
        
          x
          
            i
          
        
        ,
        
          y
          
            i
          
        
        )
      
    
    {\displaystyle f(x_{i},y_{i})}
   with 
  
    
      
        
          
            
              
                y
                
                  i
                  +
                  1
                
              
              −
              
                y
                
                  i
                
              
            
            
              
                x
                
                  i
                  +
                  1
                
              
              −
              
                x
                
                  i
                
              
            
          
        
      
    
    {\displaystyle {\frac {y_{i+1}-y_{i}}{x_{i+1}-x_{i}}}}
  . The derivative represents the instantaneous (happening continuously) change in 
  
    
      
        y
      
    
    {\displaystyle y}
   at 
  
    
      
        x
      
    
    {\displaystyle x}
  :

  
    
      
        
          lim
          
            Δ
            x
            →
            0
          
        
        
          
            
              Δ
              y
            
            
              Δ
              x
            
          
        
        =
        
          
            
              d
              y
            
            
              d
              x
            
          
        
        =
        
          y
          ′
        
        (
        x
        )
      
    
    {\displaystyle \lim _{\Delta x\to 0}{\frac {\Delta y}{\Delta x}}={\frac {dy}{dx}}=y'(x)}
  The average change in 
  
    
      
        y
      
    
    {\displaystyle y}
   over 
  
    
      
        x
      
    
    {\displaystyle x}
   is

  
    
      
        
          
            
              
                y
                
                  i
                  +
                  1
                
              
              −
              
                y
                
                  i
                
              
            
            
              
                x
                
                  i
                  +
                  1
                
              
              −
              
                x
                
                  i
                
              
            
          
        
        =
        
          
            
              Δ
              y
            
            
              Δ
              x
            
          
        
      
    
    {\displaystyle {\frac {y_{i+1}-y_{i}}{x_{i+1}-x_{i}}}={\frac {\Delta y}{\Delta x}}}
  ,which represents a discrete version of the change (changes over "steps" of time). 
If the change in 
  
    
      
        x
      
    
    {\displaystyle x}
   is small enough, either of the two can be used to approximate the other.
Assume that we know the derivative function of an unknown curve and wants to find the shape of the curve (a function that satisfies a differential equation). We can start from a known point, use the differential equation as as a formula to get the slope of the tangent line to the curve at any point on the curve, and take a small step along that tangent line up to the next point. If the step is small enough, we can expect the next point to be close to the curve. If we pretend that point is still on the curve, the same process can be used to find the next point and so on. Graphically we are using a polynomial to approximate the unknown curve. The gap between the two curves represents the error of the approximation, which can be reduced by shortening the step sizes.
We just derived the formula for Euler's method:

  
    
      
        
          y
          
            i
            +
            1
          
        
        =
        
          y
          
            i
          
        
        +
        f
        (
        
          x
          
            i
          
        
        ,
        
          y
          
            i
          
        
        )
        (
        
          x
          
            i
            +
            1
          
        
        −
        
          x
          
            i
          
        
        )
      
    
    {\displaystyle y_{i+1}=y_{i}+f(x_{i},y_{i})(x_{i+1}-x_{i})}
  


== Example ==
Given 
  
    
      
        
          
            
              d
              y
            
            
              d
              x
            
          
        
        =
        y
      
    
    {\displaystyle {\frac {dy}{dx}}=y}
  ,  and 
  
    
      
        y
        (
        0
        )
        =
        1
      
    
    {\displaystyle y(0)=1}
  , solve it for 
  
    
      
        x
      
    
    {\displaystyle x}
   at 1, i.e. y(1).
The Euler's method formula is a recursively defined function, which can be calculated iteratively. 
  
    
      
        (
        
          x
          
            i
            +
            1
          
        
        −
        
          x
          
            i
          
        
        )
      
    
    {\displaystyle (x_{i+1}-x_{i})}
   is called the step size. If we pick a step size of 1 the solution is as follows:

Another way to derive Euler's method is to use taylor expansion around 
  
    
      
        
          x
          
            i
          
        
      
    
    {\displaystyle x_{i}}
  : 

  
    
      
        y
        (
        
          x
          
            0
          
        
        +
        h
        )
        =
        y
        (
        
          x
          
            0
          
        
        )
        +
        h
        
          y
          ′
        
        (
        
          x
          
            0
          
        
        )
        +
        
          
            1
            2
          
        
        
          h
          
            2
          
        
        
          y
          ″
        
        (
        
          x
          
            0
          
        
        )
        +
        O
        (
        
          h
          
            3
          
        
        )
        .
      
    
    {\displaystyle y(x_{0}+h)=y(x_{0})+hy'(x_{0})+{\frac {1}{2}}h^{2}y''(x_{0})+O(h^{3}).}
  If we ignore higher order derivative terms we get the Euler's method formula:

  
    
      
        y
        (
        
          x
          
            0
          
        
        +
        h
        )
        =
        y
        (
        
          x
          
            0
          
        
        )
        +
        h
        
          y
          ′
        
        (
        
          x
          
            0
          
        
        )
      
    
    {\displaystyle y(x_{0}+h)=y(x_{0})+hy'(x_{0})}
  


= Runge-Kutta 2nd Order Method =
One of the Runge-Kutta 2nd order method is the midpoint method, which is a modified Euler's method (one-step method) for numerically solving ordinary differential equations:

  
    
      
        
          y
          ′
        
        (
        x
        )
        =
        f
        (
        x
        ,
        y
        (
        x
        )
        )
        ,
        
        y
        (
        
          x
          
            0
          
        
        )
        =
        
          y
          
            0
          
        
      
    
    {\displaystyle y'(x)=f(x,y(x)),\quad y(x_{0})=y_{0}}
  .The midpoint method is given by the formula

  
    
      
        
          y
          
            i
            +
            1
          
        
        =
        
          y
          
            i
          
        
        +
        h
        f
        
          (
          
            
              x
              
                i
              
            
            +
            
              
                h
                2
              
            
            ,
            
              y
              
                i
              
            
            +
            
              
                h
                2
              
            
            f
            (
            
              x
              
                i
              
            
            ,
            
              y
              
                i
              
            
            )
          
          )
        
      
    
    {\displaystyle y_{i+1}=y_{i}+hf\left(x_{i}+{\frac {h}{2}},y_{i}+{\frac {h}{2}}f(x_{i},y_{i})\right)}
  .Note that 
  
    
      
        
          y
          
            i
          
        
      
    
    {\displaystyle y_{i}}
   is the approximation of 
  
    
      
        y
        (
        
          x
          
            i
          
        
        )
      
    
    {\displaystyle y(x_{i})}
   using this method.


= Runga-Kutta 4th Order Method =
The Runga-Kutta 4th order method is often called "RK4", "classical Runge–Kutta method" or simply as "the Runge–Kutta method".

  
    
      
        
          
            
              
                
                  y
                  
                    i
                    +
                    1
                  
                
              
              
                
                =
                
                  y
                  
                    i
                  
                
                +
                
                  
                    
                      h
                      6
                    
                  
                
                
                  (
                  
                    
                      k
                      
                        1
                      
                    
                    +
                    2
                    
                      k
                      
                        2
                      
                    
                    +
                    2
                    
                      k
                      
                        3
                      
                    
                    +
                    
                      k
                      
                        4
                      
                    
                  
                  )
                
              
            
            
              
                
                  x
                  
                    i
                    +
                    1
                  
                
              
              
                
                =
                
                  x
                  
                    i
                  
                
                +
                h
              
            
          
        
      
    
    {\displaystyle {\begin{aligned}y_{i+1}&=y_{i}+{\tfrac {h}{6}}\left(k_{1}+2k_{2}+2k_{3}+k_{4}\right)\\x_{i+1}&=x_{i}+h\\\end{aligned}}}
  for i = 0, 1, 2, 3, . . . , using

  
    
      
        
          
            
              
                
                  k
                  
                    1
                  
                
              
              
                
                =
                f
                (
                
                  x
                  
                    i
                  
                
                ,
                
                  y
                  
                    i
                  
                
                )
                ,
              
            
            
              
                
                  k
                  
                    2
                  
                
              
              
                
                =
                f
                (
                
                  x
                  
                    i
                  
                
                +
                
                  
                    
                      h
                      2
                    
                  
                
                ,
                
                  y
                  
                    i
                  
                
                +
                
                  
                    
                      1
                      2
                    
                  
                
                
                  k
                  
                    1
                  
                
                h
                )
                ,
              
            
            
              
                
                  k
                  
                    3
                  
                
              
              
                
                =
                f
                (
                
                  x
                  
                    i
                  
                
                +
                
                  
                    
                      h
                      2
                    
                  
                
                ,
                
                  y
                  
                    i
                  
                
                +
                
                  
                    
                      1
                      2
                    
                  
                
                
                  k
                  
                    2
                  
                
                h
                )
                ,
              
            
            
              
                
                  k
                  
                    4
                  
                
              
              
                
                =
                f
                (
                
                  x
                  
                    i
                  
                
                +
                h
                ,
                
                  y
                  
                    i
                  
                
                +
                
                  k
                  
                    3
                  
                
                h
                )
                .
              
            
          
        
      
    
    {\displaystyle {\begin{aligned}k_{1}&=f(x_{i},y_{i}),\\k_{2}&=f(x_{i}+{\tfrac {h}{2}},y_{i}+{\tfrac {1}{2}}k_{1}h),\\k_{3}&=f(x_{i}+{\tfrac {h}{2}},y_{i}+{\tfrac {1}{2}}k_{2}h),\\k_{4}&=f(x_{i}+h,y_{i}+k_{3}h).\end{aligned}}}
  This method calculate 
  
    
      
        
          y
          
            i
            +
            1
          
        
      
    
    {\displaystyle y_{i+1}}
   by adding to 
  
    
      
        
          y
          
            i
          
        
      
    
    {\displaystyle y_{i}}
  ) a weighted average of four increments, each of which is the step size 
  
    
      
        h
      
    
    {\displaystyle h}
   multiplied by an estimated slope:

  
    
      
        
          k
          
            1
          
        
      
    
    {\displaystyle k_{1}}
   is the increment based on the estimated slope at the beginning of the interval (Euler's method) ;

  
    
      
        
          k
          
            2
          
        
      
    
    {\displaystyle k_{2}}
   is the increment based on the slope at the midpoint of the interval;

  
    
      
        
          k
          
            3
          
        
      
    
    {\displaystyle k_{3}}
   is again the increment based on the slope at the midpoint, but now using 
  
    
      
        
          k
          
            2
          
        
      
    
    {\displaystyle k_{2}}
  ;

  
    
      
        
          k
          
            4
          
        
      
    
    {\displaystyle k_{4}}
   is the increment based on the slope at the end of the interval, using 
  
    
      
        
          k
          
            3
          
        
      
    
    {\displaystyle k_{3}}
  .In averaging the four increments, greater weight is given to the increments at the midpoint. The weights are chosen such that if 
  
    
      
        f
      
    
    {\displaystyle f}
   is independent of 
  
    
      
        y
      
    
    {\displaystyle y}
  , so that the differential equation is equivalent to a simple integral.


= Case Study =
source
solution in a iPython notebook