Differentiation with a Computer

• The derivative of a function, $f(x)$ at, $x=a$, is defined as
  
  $$\frac{df}{dx} \vert_{x=a}=lim_{\Delta x \rightarrow 0 } \frac{ f(a+\Delta x)-f(a)}{\Delta x}$$
  
  
  This is called a forward-difference approximation. The limit could be approached from the opposite direction, giving a backward-difference approximation.
• It should be clear that a computer can calculate an approximation to the derivative from
  
  $$\frac{df}{dx}\vert_{x=a}=\frac{f(a+\Delta x)-f(a)}{\Delta x}$$
  
  
  if a very small value is used for $\Delta x$.
• What if we do this, recalculating with smaller and smaller values of $x$ starting from an initial value that is not small? We should expect to find an optimal value for $x$ because round-off errors in the numerator will become great as $x$ approaches zero, and these are magnified by the small value in the denominator.
• When we try this for $f(x)=e^x sin(x)$ at $x = 1.9$, starting with $\Delta x=0.05$ and halving $\Delta x$ each time, we find that the errors of the approximation decrease as $\Delta x$ is reduced until about $\Delta x = 0.05/128$.
• The analytical answer is 4.1653826. Table 1 gives the results. Notice that each successive error is about 1/2 of the previous error as $\Delta x$ is halved until $\Delta x$ gets quite small, at which time round off affects the ratio.
• At values for $\Delta x$ smaller than $0.05/128$, the error of the approximation increases due to round off. In effect, the best value for $\Delta x$ is when the effects of round-off and truncation errors are balanced.
  

  Table 1: Forward-difference approximations for $f(x)=e^x sin(x)$.

  
  
  
• If a backward-difference approximation is used; similar results are obtained.
• It should be clear that a computer can calculate an approximation to the derivative from
  
  $$\frac{df}{dx}\vert_{x=a}=\frac{f(a)-f(a-\Delta x)}{\Delta x}$$
  
  

• With MATLAB,
  • it can get the analytical answer to the function of Table 1.

    >> f='exp(x)*sin(x)'
    >> df=diff(f,'x')
    >> numeric(subs(df,1.9,'x'))
    

  • it can compute numerically;

    >> x=[1.9 1.9 1.9 1.9 1.9 1.9 1.9 1.9 1.9];
    >> del=[.05 .05/2 .05/4 .05/8 .05/16 .05/32
            .05/64 .05/128 .05/256];
    >> xplus=x+del;
    >> f=exp(x).*sin(x);
    >> fplus=exp(xplus).*sin(xplus);
    >> num=fplus-f;
    >> deriv=num./del;
    

• It is not by chance that the errors are about halved each time. Look at this Taylor series where we have used $h$ for $\Delta x$:
  
  $$f(x+h)=f(x)+f'(x)*h +f''(\xi)*h^2/2$$
  
  
  where the last term is the error. The value of $\xi$ is at some point between $x$ and $x + h$. If we solve this equation for $f'(x)$, we get
  

  $$f'(x)=(f(x+h)-f(x))/h-f''(\xi)*h^2/2$$ (1)

  
  
  which shows that the errors should be about proportional to $h$, precisely what Table 1 shows. If we repeat this but begin with the Taylor series for $f(x - h)$, it turns out that
  

  $$f'(x)=(f(x)-f(x-h))/h+f''(\zeta)*h^2/2$$ (2)

  
  
  where $\zeta$ is between $x$ and $x - h$, so the two error terms are not identical though both are $O(h)$.
• If we add Eqs. 1 and 2, then divide by 2, we get the central-difference approximation to the derivative:
  

  $$f'(x)=(f(x+h)-f(x-h))/(2h)-f'''(\xi)h^2/2$$ (3)

  
  
  We had to extend the two Taylor series by an additional term to get the error because the $f''(x)$ terms cancel.
• This shows that using a central-difference approximation is a much preferred way to estimate the derivative; even though we use the same number of computations of the function at each step, we approach the answer much more rapidly.
• Table 2 illustrates this, showing that errors decrease about four fold when $\Delta x$ is halved (as Eq. 3 predicts) and that a more accurate value is obtained.
  

  Table 2: Central-difference approximations for $f(x)=e^x sin(x)$.