Ch5.1.3 Richardson Extrapolation

Forward Difference Formula

Using Taylor polynomials, it is shown in Notes that

\[ f'(x) = \frac{f(x+h) - f(x)}{h} + k_1h + k_2h^2 + \ldots \]

Richardson extrapolation works with formulas of this form:

\[ M = N(h) + k_1h + k_2h^2 + \ldots \]

Using arithmetic, Richardson extrapolation converts \( O(h) \) results into ones that are \( O(h^2) \), \( O(h^3) \), \( O(h^4) \), and etc.

Richardson Extrapolation for O(h)

title

Richardson Extrapolation for O(h)

title

Richardson Extrapolation for O(h)

Consider the formula

\[ \begin{aligned} M &= N_1(h) + k_1h + k_2h^2 + \ldots \\ &= N_1(h) + O(h) \end{aligned} \]

The following improved result is derived in Notes:

\[ \begin{aligned} M & = N_k(h) + O(h^k) \\ N_k(h) & = N_{k-1}(h/2) + \frac{N_{k-1}(h/2)-N_{k-1}(h)}{2^{k-1}-1} \end{aligned} \]

Richardson Extrapolation for O(h)

title

Richardson Extrapolation for O(h)

title

R Code for Forward Difference

Forward difference formula:

\[ f'(x) \cong \frac{f(x+h)-f(x)}{h} \]

R code for finddif

findiff <- function(f,x,
  h = x*sqrt(.Machine$double.eps))
  {options(digits = 5) 

  return( ( f(x + h) - f(x) ) / h ) }

Richardson and Forward Difference

richardsonfwd <- function(f, x, h) {
#h = initial step size
#Initialize matrix R
R <- matrix(NA , nrow = 4, ncol = 5)
R[1,1] <- h
for (i in 2:4){R[i,1] <- R[i-1,1]/2 }
for (i in 1:4){R[i,2] <- findiff(f,x,h=R[i,1]) }
#Richardson extrapolation by column
for(j in 3:5) #Start with 3rd column & move right
 {k = j - 1   #O(h^k) results stored in column j
  for(i in k:4)    #Start at row k & go down
   {A <- R[i-1,k]  #Look left and up in matrix 
    B <- R[i,k]    #Look left in matrix
    R[i,j] <- B+(B-A)/(2^(k-1)-1)}} #Richardson
return(R) }

Example 1: Compare with Exact

findiff(cos,1)

[1] -0.84147

richardsonfwd(cos,1,0.2)

      [,1]     [,2]     [,3]     [,4]     [,5]
[1,] 0.200 -0.88972       NA       NA       NA
[2,] 0.100 -0.86706 -0.84440       NA       NA
[3,] 0.050 -0.85463 -0.84219 -0.84145       NA
[4,] 0.025 -0.84814 -0.84165 -0.84147 -0.84147

Ex 1: Column 1 Stores h Values

h <- c(0.2,0.1,0.05,0.025)
findiff(cos,1,h)

[1] -0.88972 -0.86706 -0.85463 -0.84814

richardsonfwd(cos,1,0.2)

      [,1]     [,2]     [,3]     [,4]     [,5]
[1,] 0.200 -0.88972       NA       NA       NA
[2,] 0.100 -0.86706 -0.84440       NA       NA
[3,] 0.050 -0.85463 -0.84219 -0.84145       NA
[4,] 0.025 -0.84814 -0.84165 -0.84147 -0.84147

Ex 1: Column 2 Stores O(h) Results

h <- c(0.2,0.1,0.05,0.025)
findiff(cos,1,h)

[1] -0.88972 -0.86706 -0.85463 -0.84814

richardsonfwd(cos,1,0.2)

      [,1]     [,2]     [,3]     [,4]     [,5]
[1,] 0.200 -0.88972       NA       NA       NA
[2,] 0.100 -0.86706 -0.84440       NA       NA
[3,] 0.050 -0.85463 -0.84219 -0.84145       NA
[4,] 0.025 -0.84814 -0.84165 -0.84147 -0.84147

Ex 1: Column 3 Stores O(h^2) Results

      [,1]     [,2]     [,3]     [,4]     [,5]
[1,] 0.200 -0.88972       NA       NA       NA
[2,] 0.100 -0.86706 -0.84440       NA       NA
[3,] 0.050 -0.85463 -0.84219 -0.84145       NA
[4,] 0.025 -0.84814 -0.84165 -0.84147 -0.84147

-0.86706 + (-0.86706 - -0.88972)/(2^(2-1) - 1)

[1] -0.8444

-0.85463 + (-0.85463 - -0.86706)/(2^(2-1) - 1)

[1] -0.8422

Ex 1: Column 4 Stores O(h^3) Results

      [,1]     [,2]     [,3]     [,4]     [,5]
[1,] 0.200 -0.88972       NA       NA       NA
[2,] 0.100 -0.86706 -0.84440       NA       NA
[3,] 0.050 -0.85463 -0.84219 -0.84145       NA
[4,] 0.025 -0.84814 -0.84165 -0.84147 -0.84147

-0.84219 + (-0.84219 - -0.84440)/(2^(3-1) - 1)

[1] -0.84145

Ex 1: Column 5 Stores O(h^4) Results

      [,1]     [,2]     [,3]     [,4]     [,5]
[1,] 0.200 -0.88972       NA       NA       NA
[2,] 0.100 -0.86706 -0.84440       NA       NA
[3,] 0.050 -0.85463 -0.84219 -0.84145       NA
[4,] 0.025 -0.84814 -0.84165 -0.84147 -0.84147

-0.84147 + (-0.84147 - -0.84145)/(2^(4-1) - 1)

[1] -0.84147

Richardson and Forward Difference

richardsonfwd <- function(f, x, h) {
#h = initial step size
#Initialize matrix R
R <- matrix(NA , nrow = 4, ncol = 5)
R[1,1] <- h
for (i in 2:4){R[i,1] <- R[i-1,1]/2 }
for (i in 1:4){R[i,2] <- findiff(f,x,h=R[i,1]) }
#Richardson extrapolation by column
for(j in 3:5) #Start with 3rd column & move right
 {k = j - 1   #O(h^k) results stored in column j
  for(i in k:4)    #Start at row k & go down
   {A <- R[i-1,k]  #Look left and up in matrix 
    B <- R[i,k]    #Look left in matrix
    R[i,j] <- B+(B-A)/(2^(k-1)-1)}} #Richardson
return(R) }

Centered Difference Formula

Using Taylor polynomials, it is shown in Notes that

\[ f'(x) = \frac{f(x+h) - f(x-h)}{2h} + k_1h^2 + k_2h^4 + \ldots \]

Richardson extrapolation also works with formulas like this:

\[ M = N(h) + k_1h^2 + k_2h^4 + \ldots \]

As before, Richardson extrapolation converts \( O(h^2) \) results into ones that are \( O(h^4) \), \( O(h^6) \), \( O(h^8) \), etc.

Richardson Extrapolation for O(h^2)

Consider the formula

\[ \begin{aligned} M &= N(h) + k_1h^2 + k_2h^4 + \ldots \\ &= N_1(h) + O(h^2) \end{aligned} \]

The following improved result is stated in Notes:

\[ \begin{aligned} M & = N_k(h) + O(h^k) \\ N_k(h) & = N_{k-1}(h/2) + \frac{N_{k-1}(h/2)-N_{k-1}(h)}{4^{k-1}-1} \end{aligned} \]

Richardson Extrapolation for O(h^2)

title

Richardson Extrapolation for O(h^2)

title

R Code for Centered Difference

Centered difference formula:

\[ f'(x) \cong \frac{f(x+h)-f(x-h)}{2h} \]

R code for symdiff

symdiff <- function(f,x,
  h = x*(.Machine$double.eps)^(1/3))
  {options(digits = 5)

  return( (f(x + h) - f(x - h)) / (2*h) ) }

Richardson and Centered Difference

richardsonsym <- function(f, x, h) {
#h = initial step size
#Initialize matrix R
R <- matrix(NA , nrow = 4, ncol = 5)
R[1,1] <- h
for (i in 2:4){R[i,1] <- R[i-1,1]/2 }
for (i in 1:4){R[i,2] <- symdiff(f,x,h=R[i,1]) }
#Richardson extrapolation by column
for(j in 3:5) #Start with 3rd column & move right
 {k = j - 1   #O(h^(2k) results stored in column j
  for(i in k:4)    #Start at row k & go down
   {A <- R[i-1,k]  #Look left and up in matrix 
    B <- R[i,k]    #Look left in matrix
    R[i,j] <- B+(B-A)/(4^(k-1)-1)}} #Richardson
return(R) }

Ex 2: Compare with Exact

symdiff(cos,1)

[1] -0.84147

richardsonsym(cos,1,0.2)

      [,1]     [,2]     [,3]     [,4]     [,5]
[1,] 0.200 -0.83587       NA       NA       NA
[2,] 0.100 -0.84007 -0.84147       NA       NA
[3,] 0.050 -0.84112 -0.84147 -0.84147       NA
[4,] 0.025 -0.84138 -0.84147 -0.84147 -0.84147

Ex 2: Column 5 Stores O(h^8) Results

      [,1]     [,2]     [,3]     [,4]     [,5]
[1,] 0.200 -0.83587       NA       NA       NA
[2,] 0.100 -0.84007 -0.84147       NA       NA
[3,] 0.050 -0.84112 -0.84147 -0.84147       NA
[4,] 0.025 -0.84138 -0.84147 -0.84147 -0.84147

-0.84147 + (-0.84147 - -0.84147)/(4^(4-1) - 1)

[1] -0.84147

Ex 3: Forward Difference

f <- function(x){x^2 + 1} 
findiff(f,3)

[1] 6

richardsonfwd(f,3,0.2)

      [,1]  [,2] [,3] [,4] [,5]
[1,] 0.200 6.200   NA   NA   NA
[2,] 0.100 6.100    6   NA   NA
[3,] 0.050 6.050    6    6   NA
[4,] 0.025 6.025    6    6    6

Ex 4: Centered Difference

f <- function(x){x^2 + 1} 
symdiff(f,3)

[1] 6

richardsonsym(f,3,0.2)

      [,1] [,2] [,3] [,4] [,5]
[1,] 0.200    6   NA   NA   NA
[2,] 0.100    6    6   NA   NA
[3,] 0.050    6    6    6   NA
[4,] 0.025    6    6    6    6

Ex 5: Forward Difference

f <- function(x)
  {x^8 + 2} 
findiff(f,0.85)

[1] 2.5646

title

Ex 5: Forward Difference

f <- function(x){x^8 + 2} 
findiff(f,0.85)

[1] 2.5646

richardsonfwd(f,0.85,0.2)

      [,1]   [,2]   [,3]   [,4]   [,5]
[1,] 0.200 6.0248     NA     NA     NA
[2,] 0.100 3.9093 1.7938     NA     NA
[3,] 0.050 3.1595 2.4098 2.6151     NA
[4,] 0.025 2.8447 2.5299 2.5700 2.5636

Ex 6: Centered Difference

f <- function(x){x^8 + 2} 
symdiff(f,0.85)

[1] 2.5646

richardsonsym(f,0.85,0.2)

      [,1]   [,2]   [,3]   [,4]   [,5]
[1,] 0.200 3.6140     NA     NA     NA
[2,] 0.100 2.8165 2.5507     NA     NA
[3,] 0.050 2.6270 2.5638 2.5646     NA
[4,] 0.025 2.5802 2.5646 2.5646 2.5646