-------------------------------------------------------------------------
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