In addition to electronic differentiation, this project explores methods for numerically differentiating the I-V data.  The first method employed was a cubic spline interpolation which fits a different third degree polynomial between every two data points.  It also imposes the constraints that the derivative must be smooth at zone boundaries and the second derivative must be continuous.   The solution requires solving a tridiagonal system of equations.  The code is implemented in C++ with an option to average the results to eliminate oscillations due to noise.  Numerical differentiation was also performed by fitting a high order polynomial to the entire data set using a least squares method. The resulting function can be exactly differentiated.   Both Cramer's rule and Gauss-Jordan Elimination were used to solve the linear system of equations that minimize the error.  The implementation consists of a C++ matrix class as well as a driver program.

The cubic spline method fits cubic polynomials to between consecutive data points, called a "zone" for this discussion.  The picture illustrates that a given zone is scaled from x_m to x_p (short for x_minus to x_plus).  At the zone boundaries (x_m and x_p) the polynomial must equal the given values (ym and yp).  This constraint would produce a continuous function but the derivatives would not be continuous.  By adding the constraints that the first derivative at xm equals the first derivative in the zone to the left at xp coupled with a similar constraint for the second derivatives, one may guarantee a smooth derivative and continuous second derivative. 

Figure 1: Diagram of a single zone in the cubic spline method. 

Figure 2: The basic equations defining the cubic spline method in a single zone.

The basic equations necessary to construct a cubic spline method are given in figure 2.  However, one notes that the known variables are xm, xp, ym, and yp, while there are 8 unknown quantities: A, B, C, D, ym', yp', ym'', and  yp''.  Since only six equations are specified it is not a trivial system of linear equations.   However, the situation could be worse, for if one is given n data points there are n-1 zones in the entire function.  This equates to 6(n-1) equations and 6(n-1)+2 unknowns since ym' and ym'' in one zone equal yp' and yp'' in the zone to the right.  The extra 2 unknowns arise since the left most boundary and right most boundary are not repeated. 

To solve this tridiagonal system of equations one must assume two values.  A logical choice is to assume the value of ym'' in the left most zone of the entire domain and yp'' in the right most zone.  Then, starting in the left most zone, one can solve for yp'' in terms of only one unknown variable.  In this case, the method solves for yp'' = AA*yp' + BB, where AA and BB are constants.  Since the first and second derivatives are equal at zone edges, in the next zone to the right we have successfully eliminated ym'' and the exact same procedure may be used to tabulate the constants AA and BB in each zone.  Once the right most zone is reached, one assumes the value of yp'' which allows direct calculation of the coefficients for the cubic polynomial in that zone.   Furthermore, this allows calculation of ym'' which equals yp'' of the zone to the left.  Therefore, the same procedure is used to calculate all the coefficients from the right most zone to the left most zone. 

Cubic Spline Implementation Issues

A cubic spline method was implemented using Borland C++ version 5.0 compiler.  The left-most and right-most second derivatives which need to assumed are actually computed by fitting a cubic polynomial to 4 data points at either edge of the domain.  The input/output uses the fstream.h library and expects the input to be a text file with each line contains one data point (x,y) separated by a tab.  The output is a text file with every line containing the four values x, y, y', y'' also tab delimited. 

Solving the linear series of equations was done with the Mathematica software from Wolfram.  The resulting functions were than converted into C++ syntax using another C++ program.  (I did not include the code because it is limited to certain types of functions and could cause catastrophic disasters in other situations.)  This machinery prevented simple errors that would inevitably occur in such calculations.   Keep in mind that this program works for arbitrarily spaced x data, instead of the typical uniform grid that simplifies most interpolation algorithms.

One problem experienced with this algorithm is osscilations due to random noise.   This becomes a greater problem as the x data values become closer together.   Therefore, the algorithm works better with less data input.  To combat this problem I've implemented a full set averaging system.  The program prompts for the number of averages to perform.  It than uses only part of the data to do the cubic spline and averages the results with the cubic spline on the remaining data points.  For example, if 2 averages are specified it completes a cubic spline interpolation on every other data point and averages the results with the cubic spline on the remaining points.  While this fixes some the eradicate behavior it is not extremely successful and it ruins the benefits of the algorithm. These problems are discussed in on the results page, with illustrations.

This method fits any order polynomial to the entire data set by minimizing the chi^2, which is the goodness of fit parameter.  The can be done directly by solving a system of degree+1 linear equations, where degree is the order of the fitting polynomial. Two methods for solving the system of equations were explored: Cramer's method and Gauss-Jordan elimination. 

Cramer's rule is easier to code, for it simply requires one to take determinants of matrices. In this case, the determinant is obtained using a recursive algorithm that constructs the cofactors of the original matrix.  It was implemented as a matrix class in C++.  This method is highly susceptible to numerical errors, which are discussed on the results page.

One may also solve the matrix equations for the coefficients of the fitting polynomial.   In particular, this method requires the inverse operation in the matrix class, which is implemented with Gauss-Jordan elimination.  Since the inverse of a matrix is the unity matrix divided by the original matrix, this method converts the original matrix to the unity and performs the same operations on the unity matrix.  This results in the unity matrix being converted into the inverse of the original matrix. 

These methods are combined in a program called polyfit, which is currently set up to use Gauss-Jordan elimination.  The program reads the input data in a form exactly the same of the previously discussed cubic spline program.  It prompts for the number of points to read and the number of points to output.  The output is also the same format as earlier, consisting of the x and y data, as well as dy/dx and d2y/dx2.  The exact derivative values are found by simply differentiating the fitting polynomial.

Note: This section is rather brief, so feel free to email any questions to enzx0002@tc.umn.edu .  Also, one may consult Bevington, Philip and D. Keith Robinson. Data Reduction and Error Analysis for the Physical Sciences - Second Edition. McGraw-Hill, Inc. 1992.

Please send comments, criticisms, queries or congratulations to Michael Enz at enzx0002@tc.umn.edu.   This page was created 5/4/98 and last modified on 5/18/98.  It will be updated as work progress through June 1998.