---
title: "Richardson Extrapolation"
author: "Ryan McClarren"
date: "September 4, 2014"
output: html_document
---

```{r, echo=FALSE,message=FALSE}
require(plyr)
require(ggplot2)
require(reshape2)
```
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
```{r}
func <- function(x) { atan(x) * cosh(x)}
```

Now let's plot the function near 1:

```{r}
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

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

The plot of the derivative looks like

```{r}
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:

```{r}
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))
```

We see the expected second-order accuracy. We then apply Richardson extrapolation
```{r}
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()

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

Fourth-order accuracy!

Now if we re-apply Richardson we can do better.
```{r}
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()

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

Sixth-order FTW!

Last time:


```{r}
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()

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

8th order accuracy might be overkill.