# Numerical differentiation using the y-intercept

I have a routine for determining the derivative of a function using the y-intercept (B) to infer the finite difference step (h). Instead of directly using smaller and smaller h, we can solve for the slope m in y=mx+B with a given B. That B represents a specific h after some substituting and simplifying where B = f(x) - x((f(x)-f(x-h))/h).

I would like to know where & how to best increase the precision. I am currently using uniroot to solve for h while an exact solution would be preferred. Thanks for other suggestions and insights!

Finite Step Routine:

Finite.step = function(point,h){

f.x = f(point)
f.x.h.min = f(point - h)
f.x.h.pos = f(point + h)

neg.step = (f.x - f.x.h.min)/h
pos.step = (f.x.h.pos - f.x)/h

return((c("f(x-h)"=neg.step,"f(x+h)"=pos.step,mean(c(neg.step,pos.step)))))

}


Numerical Differentiation Routine:

Numerical.differentiation = function(point,h,tol=1e-10,print.trace="FALSE"){

options(digits=20)
Bs = numeric()
Bl = numeric()
Bu = numeric()

### Step 1 initalize the boundaries for B

### Use the initial step size h
f.x = f(point)
f.x.h = f(point - h)

### Y = mX + B
Y = f.x
m = (f.x - f.x.h)/h
mX = ((f.x - f.x.h)/h)*point
B = f.x - ((f.x - f.x.h)/h)*point

### Initial interval for B given inputted h-step-value

f.x.h.lower = f(point - h)
f.x.h.upper = f(point + h)

B1 = f.x - ((f.x - f.x.h.lower)/h)*point
B2 = f.x - ((f.x.h.upper-f.x)/h)*point

low.B = min(c(B1,B2))
high.B = max(c(B1,B2))

lower.B = low.B
upper.B = high.B

## Return "N/A" if lower.B and upper.B are identical to 20 digits
if(lower.B==upper.B){
return(c("No Derivative Exists"))}

new.B = mean(c(lower.B,upper.B))

i=1
while(i>=1L){
Bl[i] = lower.B
Bu[i] = upper.B
new.f = function(x) -f.x + ((f.x - f(point - x))/x)*point + new.B

###  Step #2 solve FOR h.  This can be better served with an exact solution

inferred.h = uniroot(new.f, c(-2*h,2*h))\$root

if(print.trace=="TRUE") {print(c("h"=inferred.h,"Lower B" = lower.B,"Upper B" = upper.B))}

Bs[i]=new.B
### Step 3 Stop if the inferred h is within the tolerance level

if(abs(inferred.h)<tol) {
final.B = mean(c(upper.B,lower.B))
slope = solve(point, f.x-final.B)
z = complex(real = point, imaginary = h)
par(mfrow=c(1,3))

### Plot #1
plot(f, xlim = c(min(c(point-(100*h), point+(100*h)),0),max(c(point-(100*h), point+(100*h)),0)),
col='azure4',
ylab='f(x)',
lwd=2,
ylim = c(min(c(min(c(B1,B2)),min(na.omit(f((point-(100*h)):(point+(100*h))))))),
max(c(max(na.omit(f((point-(100*h)):(point+(100*h))))),
max(c(B1,B2))))),
main='f(x) and initial y-intercept range')

abline(h=0,v=0,col='grey')
points(point,f.x, pch=19,col='green')
points(point-h,f.x.h.lower,col=ifelse(B1==high.B,'blue','red'),pch=19)
points(point+h,f.x.h.upper,col=ifelse(B1==high.B,'red','blue'),pch=19)
points(x=rep(0,2),y=c(B1,B2),col=c(ifelse(B1==high.B,'blue','red'),ifelse(B1==high.B,'red','blue')),pch=1)
segments(0,B1,point-h,f.x.h.lower,col=ifelse(B1==high.B,'blue','red'),lty=2)
segments(0,B2,point+h,f.x.h.upper,col=ifelse(B1==high.B,'red','blue'),lty=2)

### Plot #2
plot(f,col='azure4',ylab='f(x)',lwd=3,main='f(x) narrowed range and secant lines',
xlim = c(min(c(point-h,point+h,0)),max(c(point+h,point-h,0))),
ylim= c(min(c(B1,B2,f.x.h.lower,f.x.h.upper)),max(c(B1,B2,f.x.h.lower,f.x.h.upper))))

abline(h=0,v=0,col='grey')
points(point,f.x, pch=19,col='green')
points(point-h,f.x.h.lower,col=ifelse(B1==high.B,'blue','red'),pch=19)
points(point+h,f.x.h.upper,col=ifelse(B1==high.B,'red','blue'),pch=19)
segments(0,B1,point-h,f.x.h.lower,col=ifelse(B1==high.B,'blue','red'),lty=2)
segments(0,B2,point+h,f.x.h.upper,col=ifelse(B1==high.B,'red','blue'),lty=2)
points(x=rep(0,2),y=c(B1,B2),col=c(ifelse(B1==high.B,'blue','red'),ifelse(B1==high.B,'red','blue')),pch=1)

### Plot #3
plot(Bs,ylim=c(min(c(Bl,Bu)),max(c(Bl,Bu))),
xlab="Iterations",
ylab="y-inetercept",
col='green',pch=19,
main='Iterated range of y-intercept')
points(Bl,col='red',ylab='')
points(Bu,col='blue',ylab='')

legend('topright',c("Upper y-intercept","Lower y-intercept","Mean y-intercept"),col= c('blue','red','green'),pch=c(1,1,19))

return(as.matrix(c("Value of f(x) at point"=f(point),
"Final y-intercept (B)" = final.B,
"DERIVATIVE"=slope,
"Inferred h" = inferred.h,
"iterations"=i,
Finite.step(point,h)[1:2],
"Averaged Finite Step Initial h "=Finite.step(point,h)[3],
"Inferred h"=Finite.step(point,inferred.h)[1:2],
"Inferred h Averaged Finite Step"=Finite.step(point,inferred.h)[3],
"Complex Step Derivative (Initial h)" = Im(f(z))/Im(z))))
}

### Step 4 Narrow the range of B based on the sign of the inferred.h
if(B1==high.B){
if(sign(inferred.h) < 0) {
lower.B = new.B
upper.B = upper.B
}
else {
upper.B = new.B
lower.B = lower.B
}}  else {
if(sign(inferred.h) < 0) {
lower.B = lower.B
upper.B = new.B
}
else {
upper.B = upper.B
lower.B = new.B
}}

new.B = mean(c(lower.B,upper.B))

i = i+1
}
}


Image above generated with:

f=function(x) sin(x)/x
Numerical.differentiation(4.1,.1)
[,1]
Value of f(x) at point              -1.9957978318644154e-01
Final y-intercept (B)                1.7566181876316217e-01
DERIVATIVE                          -9.1522341938927748e-02
Inferred h                           3.7079016940166909e-11
iterations                           2.8000000000000000e+01
f(x-h)                              -1.0379159359459517e-01
f(x+h)                              -7.9382578644127122e-02
Averaged Finite Step Initial h      -9.1587086119361144e-02
Inferred h.f(x-h)                   -9.1521721250617208e-02
Inferred h.f(x+h)                   -9.1523218354796668e-02
Inferred h Averaged Finite Step     -9.1522469802706938e-02
Complex Step Derivative (Initial h) -9.1522966669957975e-02