12-267/Numerical Methods

Summary of Numerical Methods

Based largely off of a note available here posted by Simon1 --Twine 20:55, 25 October 2012 (EDT)

Numerical methods: ${\frac {dy}{dt}}=f(t,y)$ and $y(t_{0})=y_{0}$, $y=\Phi (t)$ is a solution.

1. Using the proof of Picard's Theorem:

$\Phi _{0}(x)=y_{0}$

$\Phi _{n}(x)=y_{0}+\int _{x_{0}}^{x}f(x,\Phi _{n-1}(x))dx$

$\Phi _{n}(x)\rightarrow \Phi (x)$

2. The Euler Method:

$y_{n+1}=y_{n}+f(t_{n},y_{n})(t_{n+1}-t_{n})=y_{n}+f_{n}h$ if h is constant

Backward Euler formula: $y_{n+1}=y_{n}+hf(t_{n+1},y_{n+1})$

Local truncation error: $e_{n+1}={\frac {1}{2}}\Phi ''(t_{n})h^{2}\leq {\frac {Mh^{2}}{2}}$ where $m=\mathrm {max} |\Phi ''(t)|$

Local error is proportional to $h^{2}$.

Global error is proportional to h.

3. Improved Euler Formula (or Heun Formula):

$y_{n+1}=y_{n}+{\frac {f_{n}+f(t_{n}+h,y_{n}+hf_{n})}{2}}h$

Local truncation error is proportional to $h^{3}$

Global truncation error is proportional to $h^{2}$

4. The Runge-Kutta Method:

$y_{n+1}=y_{n}+{\frac {k_{n1}+2k_{n2}+2k_{n3}+k_{n4}}{6}}h$

where

$k_{n1}=f(t_{n},y_{n})\quad k_{n2}=f(t_{n}+{\frac {1}{2}}h,y_{n}+{\frac {1}{2}}hk_{n1})$

$k_{n3}=f(t_{n}+{\frac {1}{2}}h,y_{n}+{\frac {1}{2}}hk_{n2})\quad k_{n4}=f(t_{n}+h,y_{n}+hk_{n3})$

Local truncation error is proportional to $h^{5}$.

Global truncation error is proportional to $h^{4}$.

Python Example of Euler's Method

In class on October 15th we discussed Euler's Method to numerically compute a solution to a differential equation. $x_{0}$ and $y_{0}$ are given as well as an increment amount $h$, $x_{n+1}=x_{n}+h$, and we use the guess $y_{n+1}=y_{n}+f(x_{n},y_{n})*h$ where f computes the derivative as a function of x and y.

Here is an example of code (written in Python) which carries out Euler's Method for the example we discussed in class, $y'=-y$:

   def f(x, y):
       return -y
   
   def euler(x, y, f, h, x_max):
       """Take in coordinates x and y, a function f(x, y) which calculates
          dy/dx at (x, y), an increment h, and a maximum value of x.
          
          Return a list containing coordinates in the Euler's Method computation
          of the solution to Phi' = f(x, Phi(x)), Phi(x) = y, with the x
          values of those coordinates separated by h, and not exceeding x_max.
       """
       if x > x_max: # we have already calculated all our values
           return []
       x_next, y_next = (x + h, y + f(x, y)*h) # calculate the next x, y values
       # return the current coordinates, and every coordinates following it, in a list
       return [(x_next, y_next)] + euler(x_next, y_next, f, h, x_max) 
   
   if __name__ == '__main__':
       print euler(0, 1, f, 0.01, 1)[-1]