-------------------------------------------------------------------------
PHYS210 2009-10-29 LAB SUMMARY

1) Interactive application of O(h) forward finite difference approximation
   (FDA) of first derivative

2) Writing an octave function to implement the FDA formula

3) Extending the octave function to handle vector of deltax's

4) Examining convergence behaviour of approximation using scaled errors

NOTE: In the following, '>>' denotes the octave prompt

+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1) Interactive application of O(h) forward finite difference approximation
   (FDA) of first derivative
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Change to your ~/octave directory

Start octave, and use the formula for the O(h) forward FDA of the 
first derivative

      f'(x) ~ f(x + delx) - f(x)
              -----------------       (1)
                    delx

to compute an approximation to the first derivative of sin(x) at 
x=0.2.  Use long format so that numbers are displayed to full
precision.

Start with delx = 0.1

Note that the exact value of the derivative is cos(0.2)

% cd ~/octave
% octave

>> format long
>> % Assign the exact value of the derivative to xct
>> xct = cos(0.2)
xct =  0.980066577841242

>> % Evaluate the formula for delx = 0.1 and assign the result to d
>> d = (sin(0.2 + 0.1) - sin(0.2)) / 0.1
d =  0.968508758662784

>> % Compute the error and assign the result to e
>> e = xct - d
e =  0.0115578191784578

>> Using the command-recall and editing functions re-evaluate for 
>> delx = 0.01, 0.001, 0.00001, 0.000001 and 0.0000001

>> d = (sin(0.2 + 0.01) - sin(0.2)) / 0.01
d =  0.979056905103837

>> e = xct - d
e =  0.00100967273740493

>> d = (sin(0.2 + 0.001) - sin(0.2)) / 0.001
d =  0.979967079839716

>> e = xct - d
e =  9.94980015252001e-05

>> d = (sin(0.2 + 0.0001) - sin(0.2)) / 0.0001
d =  0.980056642741201

>> e = xct - d
e =  9.93510004110298e-06

>> d = (sin(0.2 + 0.00001) - sin(0.2)) / 0.00001
d =  0.980065584479939

>> e = xct - d
e =  9.93361302881191e-07

>> d = (sin(0.2 + 0.000001) - sin(0.2)) / 0.000001
d =  0.980066478528663

>> e = xct - d
e =  9.93125788273375e-08

>> d = (sin(0.2 + 0.0000001) - sin(0.2)) / 0.0000001
d =  0.980066567901616

>> e = xct - d
e =  9.93962534501236e-09

Observe how the error decreases by very nearly a factor of 10 each time
that delx is reduced by a factor of 10.  

This provides "empirical verification" that the FDA (1) *is* first 
order accurate, i.e. that the leading order error term is *linear* 
in delx.

+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2) Writing an octave function to implement the FDA formula
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
In your octave directory, ~/octave, create an octave function file
fd.m that defines a function fd with the following header

function d = fd(fcn, x, delx)

The input arguments to fd are as follows
   1) fcn:  A FUNCTION HANDLE, which will typically be a handle of one 
            of octave's built in math functions.  In octave (and MATLAB)
            one creates a function handle by applying the @ operator to 
            the name of a function.
   2) x:    The scalar value of x at which the derivative is to be 
            approximately evaluated.
   3) delx: The finite spacing, delx appearing in (1)

The output argument from fd is as follows
   1) d:   The approximate value of d fcn(x) / dx evaluated using formula 
           (1)

SOLUTION

   function d = fd(fcn, x, delx)
   % Implements O(h) forward FDA of 1st derivative of fcn.
      d = (fcn(x+delx) - fcn(x)) / delx;
   end

IMPORTANT!!  Remember that output arguments in octave/MATLAB functions
are to be treated as variables that MUST be assigned values before 
the end of the function is reached (or an explicit return statement
is encountered)

Test the function.

Ensuring that you are in your ~/octave directory, start octave, and 
try the following---note that '@sin' passes the function handle to 
'sin' to 'fd', so that when 'fcn' is invoked in 'fd', it is actually
'sin' that is being called.

% cd ~/octave
% octave

>> format long 
>> xct = cos(0.2)
xct =  0.980066577841242

>> d = fd(@sin, 0.2, 0.1)
d =  0.968508758662784

>> e = xct - d
e =  0.0115578191784578

>> % Use a for loop to evaluate the formula using the same set of 
>> % delx's as previously, use 'disp('')' to generate an empty line
>> for delx = [0.1 0.01 0.001 0.0001 0.00001 0.000001 0.0000001]
       d = fd(@sin, 0.2, delx)
       e = xct - d
       disp('')
   end

d =  0.968508758662784
e =  0.0115578191784578

d =  0.979056905103837
e =  0.00100967273740493

d =  0.979967079839716
e =  9.94980015252001e-05

d =  0.980056642741201
e =  9.93510004110298e-06

d =  0.980065584479939
e =  9.93361302881191e-07

d =  0.980066478528663
e =  9.93125788273375e-08

d =  0.980066567901616
e =  9.93962534501236e-09

% Use 'fd' to compute an approximation to the derivative of another 
% function: d(sqrt(x))/dx at x = 4
%
% The exact answer in this case is (1/2)x^{-1/2} evaluated at x = 4
% which is 1/4

>> format long
>> xct = 0.25

>> for delx = [0.1 0.01 0.001 0.0001 0.00001 0.000001 0.0000001]
       d = fd(@sqrt, 4, delx)
       e = xct - d
       disp('')
   end

d =  0.248456731316584
e =  0.00154326868341581

d =  0.249843945007866
e =  1.56054992134003e-04

d =  0.249984376953005
e =  1.56230469947616e-05

d =  0.249998437520382
e =  1.56247961768941e-06

d =  0.249999843759952
e =  1.56240048482248e-07

d =  0.249999984269778
e =  1.57302224579325e-08

d =  0.249999998480632
e =  1.51936774273054e-09

So again we observe the expected linear convergence of the error with 
delx.
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
3) Extending the octave function to handle vector of delx's
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Now, modify fd.m so that fd will return a row vector of approximations
if the input argument delx is a row vector

SOLUTION

   function d = fd(fcn, x, delx)
   % Implements O(h) forward FDA of 1st derivative of fcn.
      d = (fcn(x+delx) - fcn(x)) ./ delx;
   end

Note that the modification was extremely easy.  As long as x is a 
scalar, then if delx is a vector, fcn(x+delx) - fcn(x) will also be a 
vector of the same length (provided that fcn returns a vector of values
when supplied with a vector as its argument).  We then simply have 
to make the division by delx an element-by-element operation (./) and 
then the return argument d will be a row vector that contains 
various approximations to the derivative, according to the specific
values supplied in delx.

Usage example

>> format long

>> % First, define a row vector defining the delx's we have been using
>> % so far. 
>> vdelx = 10.^[-1:-1:-7]

vdelx =

 Columns 1 through 3:

   1.00000000000000e-01   1.00000000000000e-02   1.00000000000000e-03

 Columns 4 through 6:

   1.00000000000000e-04   1.00000000000000e-05   1.00000000000000e-06

 Column 7:

   1.00000000000000e-07


>> % IMPORTANT!! Note that we MUST use .^ in the above, ^ won't work 
>> % as you can verify

>> % Compute the vector of approximations
>> vd = fd(@sin, 0.2, vdelx)
vd =

 Columns 1 through 4:

   0.968508758662784   0.979056905103837   0.979967079839716   0.980056642741201

 Columns 5 through 7:

   0.980065584479939   0.980066478528663   0.980066567901616


>> % ... and the vector of errors
>> e = cos(0.2) - vd

e =

 Columns 1 through 3:

   1.15578191784578e-02   1.00967273740493e-03   9.94980015252001e-05

 Columns 4 through 6:

   9.93510004110298e-06   9.93361302881191e-07   9.93125788273375e-08

 Column 7:

   9.93962534501236e-09


+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
4) Examining convergence behaviour of approximation using scaled
   errors
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

>> format long

>> % Define another vector of delx's that are reciprocal powers of 2
>> vdelx = 2.^(-1*[1:20])
vdelx =

 Columns 1 through 3:

   5.00000000000000e-01   2.50000000000000e-01   1.25000000000000e-01

 Columns 4 through 6:

   6.25000000000000e-02   3.12500000000000e-02   1.56250000000000e-02

 Columns 7 through 9:

   7.81250000000000e-03   3.90625000000000e-03   1.95312500000000e-03

 Columns 10 through 12:

   9.76562500000000e-04   4.88281250000000e-04   2.44140625000000e-04

 Columns 13 through 15:

   1.22070312500000e-04   6.10351562500000e-05   3.05175781250000e-05

 Columns 16 through 18:

   1.52587890625000e-05   7.62939453125000e-06   3.81469726562500e-06

 Columns 19 and 20:

   1.90734863281250e-06   9.53674316406250e-07

>> % Compute vector of approximations ...
>> vd = fd(@sin, 0.2, vdelx)

vd =

 Columns 1 through 4:

   0.891096712885260   0.945184813264676   0.965115640495518   0.973222242391746

 Columns 5 through 8:

   0.976803113904581   0.978474626747447   0.979280559992660   0.979676059861639

 Columns 9 through 12:

   0.979871941775130   0.979969415562408   0.980018035643297   0.980042316478034

 Columns 13 through 16:

   0.980054449593581   0.980060514325942   0.980063546236124   0.980065062076392

 Columns 17 through 20:

   0.980065819971060   0.980066198906570   0.980066388379782   0.980066483112751

>> % ... and vector of errors.
>> evd = cos(0.2) - vd
evd =

 Columns 1 through 3:

   8.89698649559820e-02   3.48817645765656e-02   1.49509373457237e-02

 Columns 4 through 6:

   6.84433544949514e-03   3.26346393666044e-03   1.59195109379451e-03

 Columns 7 through 9:

   7.86017848582121e-04   3.90517979602878e-04   1.94636066111586e-04

 Columns 10 through 12:

   9.71622788334958e-05   4.85421979441458e-05   2.42613632075450e-05

 Columns 13 through 15:

   1.21282476606144e-05   6.06351529985893e-06   3.03160511738731e-06

 Columns 16 through 18:

   1.51576484985760e-06   7.57870181944398e-07   3.78934671418918e-07

 Columns 19 and 20:

   1.89461459187967e-07   9.47284910512991e-08

>> % Finally, scale these errors by powers of 2 which increase as 
>> % delx decreases.  The fact that the numbers that result are almost 
>> % the same, especially near the end of the vector makes it easier to 
>> % see that the convergence is very close to linear in delx as 
>> % expected 

>> sevd = evd .* 2.^[1:20]
sevd =

 Columns 1 through 3:

   0.1779397299119641   0.1395270583062622   0.1196074987657898

 Columns 4 through 6:

   0.1095093671919223   0.1044308459731340   0.1018848700028485

 Columns 7 through 9:

   0.1006102846185115   0.0999726027783367   0.0996536658491323

 Columns 10 through 12:

   0.0994941735254997   0.0994144213896107   0.0993745436981044

 Columns 13 through 15:

   0.0993546048357530   0.0993446346728888   0.0993396364865475

 Columns 16 through 18:

   0.0993371652002679   0.0993355604878161   0.0993354505044408

 Columns 19 and 20:

   0.0993323695147410   0.0993300222326070