2
\$\begingroup\$

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!

enter image description here

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
\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.