Richardson Extrapolation

Leave a reply

Richardson extrapolation can produce highly accurate estimates of numerical quantities by using knowledge of the magnitude of the numerical error.  In this markdown (source here) we apply this to the finite difference calculation of a derivative.

Richardson Extrapolation

Let’s apply Richardson extrapolation to the derivative of arctan(x) cosh(x) at x = 1. The actual value of this is 1.694541176517952557683135.

First, we’ll define the function in R

func <- function(x) { atan(x) * cosh(x)}

Now let’s plot the function near 1:

N = 500
plotData = data.frame(x=0,f=0)
plotData[1:N,'x'] = seq(0.1,2.5,length.out=N)
plotData$f = func(plotData$x)
ggplot(plotData, aes(x=x, y=f)) + geom_line() + scale_x_continuous('x') + scale_y_continuous('f(x)')

The finite difference formula we want to use is

dfdx <- function(f,x,delta) {( f(x+delta/2) - f(x - delta/2))/delta}

The plot of the derivative looks like

plotData[1:N,'dfdx1'] = dfdx(func,plotData$x,0.05)
meltDF <- melt(plotData,id='x')
ggplot(meltDF, aes(x=x, y=value, color=variable)) + geom_line() + scale_x_continuous('x') + scale_y_continuous('f(x)')

The next step is to figure out how the derivative changes with respect to delta:

exact = 1.694541176517952557683135
dx = 2^seq(0,-5)
derivTable = data.frame(dx = dx, dfdx1 = dfdx(func,1,dx))
derivTable$err = exact - derivTable$dfdx1
ggplot(derivTable, aes(x = dx, y = abs(err))) + geom_point() + scale_x_log10() + scale_y_log10()

summary(lm( log(abs(err)) ~ log(dx), data = derivTable))

## 
## Call:
## lm(formula = log(abs(err)) ~ log(dx), data = derivTable)
## 
## Residuals:
##         1         2         3         4         5         6 
##  0.003711 -0.002189 -0.002500 -0.001297  0.000294  0.001981 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -2.362114   0.001999   -1182  3.1e-12 ***
## log(dx)      2.002481   0.000952    2103  3.1e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.00276 on 4 degrees of freedom
## Multiple R-squared:     1,   Adjusted R-squared:     1 
## F-statistic: 4.42e+06 on 1 and 4 DF,  p-value: 3.07e-13

We see the expected second-order accuracy. We then apply Richardson extrapolation

RichExtrap <- function(f1, f2, p,n){
  (n^p*f1 - f2)/(n^p-1)
}
derivTable$dfdx2[2:6] = RichExtrap(derivTable$dfdx1[2:6], derivTable$dfdx1[1:5], 2, 2)
derivTable$err2 = exact - derivTable$dfdx2
meltDF <- melt(derivTable[,c('dx','err','err2')], id=dx)
ggplot(meltDF, aes(x=dx,y=abs(value),color=variable)) + geom_point() + scale_x_log10() + scale_y_log10()

## Warning: Removed 1 rows containing missing values (geom_point).

summary(lm( log(abs(err2)) ~ log(dx), data = derivTable))

## 
## Call:
## lm(formula = log(abs(err2)) ~ log(dx), data = derivTable)
## 
## Residuals:
##        2        3        4        5        6 
## -0.02280  0.01874  0.01632  0.00233 -0.01459 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -5.5604     0.0223    -249  1.4e-07 ***
## log(dx)       3.9742     0.0097     410  3.2e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.0213 on 3 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:     1,   Adjusted R-squared:     1 
## F-statistic: 1.68e+05 on 1 and 3 DF,  p-value: 3.2e-08

Fourth-order accuracy!

Now if we re-apply Richardson we can do better.

derivTable$dfdx3[3:6] = RichExtrap(derivTable$dfdx2[3:6], derivTable$dfdx2[2:5], 4, 2)
derivTable$err3 = exact - derivTable$dfdx3
meltDF <- melt(derivTable[,c('dx','err','err2','err3')], id=dx)
ggplot(meltDF, aes(x=dx,y=abs(value),color=variable)) + geom_point() + scale_x_log10() + scale_y_log10() + geom_line()

## Warning: Removed 3 rows containing missing values (geom_point).
## Warning: Removed 3 rows containing missing values (geom_path).

summary(lm( log(abs(err3)) ~ log(dx), data = derivTable))

## 
## Call:
## lm(formula = log(abs(err3)) ~ log(dx), data = derivTable)
## 
## Residuals:
##       3       4       5       6 
## -0.0209  0.0259  0.0108 -0.0159 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -5.5629     0.0446    -125  6.4e-05 ***
## log(dx)       5.9549     0.0175     340  8.7e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.0271 on 2 degrees of freedom
##   (2 observations deleted due to missingness)
## Multiple R-squared:     1,   Adjusted R-squared:     1 
## F-statistic: 1.16e+05 on 1 and 2 DF,  p-value: 8.65e-06

Sixth-order FTW!

Last time:

derivTable$dfdx4[4:6] = RichExtrap(derivTable$dfdx3[4:6], derivTable$dfdx3[3:5], 6, 2)
derivTable$err4 = exact - derivTable$dfdx4
meltDF <- melt(derivTable[,c('dx','err','err2','err3','err4')], id=dx)
ggplot(meltDF, aes(x=dx,y=abs(value),color=variable)) + geom_point() + scale_x_log10() + scale_y_log10() +geom_line()

## Warning: Removed 6 rows containing missing values (geom_point).
## Warning: Removed 6 rows containing missing values (geom_path).

summary(lm( log(abs(err4)) ~ log(dx), data = derivTable))

## 
## Call:
## lm(formula = log(abs(err4)) ~ log(dx), data = derivTable)
## 
## Residuals:
##       4       5       6 
##  0.0453 -0.0906  0.0453 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)   -3.896      0.320   -12.2    0.052 . 
## log(dx)        8.004      0.113    70.7    0.009 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.111 on 1 degrees of freedom
##   (3 observations deleted due to missingness)
## Multiple R-squared:     1,   Adjusted R-squared:     1 
## F-statistic: 5e+03 on 1 and 1 DF,  p-value: 0.00901

8th order accuracy might be overkill.

This entry was posted in R Markdowns, Verification on September 4, 2014 by .