Curve and surface algorithm are an important topic in SimplexNumerica for geometric modelling and visualization courses. In addition to manually set Bezier and B-Splines from previous Pulldown-menu, interpolation and approximation provides another level of sophistication. In many situations such as surface re-engineering and facial movement animation, you may specify a set of data points that describes a desired shape (e.g., surface model) through any probing or scanning, and obtain a surface that contains all data points. Interpolation is also important in computer animation. An animator may specify a number of key camera positions and orientations (i.e., key frames), interpolate these positions with a B-spline curve (i.e., camera path), and interpolate the key frames with additional frames. While interpolation can produce a curve/surface that follows the shape of the data points, it may oscillate or wiggle its way through every point. Approximation can overcome this problem so that the curve/surface still captures the shape of the data points without containing all of them.

• Regression
  
• Interpolation
  
• Spline Interpolation
  
• Approximation
  
• Fourier Analysis of Time Series
  
• Bezier & Spline Approximation
  
• Numerical Integration
  
• Numerical Differentiation
  

 Consider the nature of most experimental data. Typically such data include noise due to many different effects. The noisy data from an experiment might appear as shown in the following Table and Figure. We assume that the x values are accurate. Visual inspection of the data suggests a positive relationship between x and y = f(x), i.e., higher values of y are associated with higher values of x. One strategy for deriving an approximating function for this data might be to try to fit the general trend of the data without necessarily matching the individual points. A straight line could be used to generally characterize the trend in the data without passing through any particular point. The line in next Figure has been sketched through the points. Although this approach may work well in many cases, it does not provide us with any quantitative measure of how good the fit of the line is to the data. We need a criterion with which to measure the goodness of fit of the line to the data. One way to do this is to derive a curve that minimizes the discrepancy between the data points and the curve. The technique for accomplishing this is called least-squares regression.

Often data are available at discrete points and we require estimates at points between the discrete values. In this section we will discuss techniques to fit curves to data in order to estimate intermediate, or fitted, values. Two methods of curve fitting are generally considered, depending on the amount of error in the data. When the data are known to be precise, the method of interpolation is used. The primary purpose of interpolation is to provide information between tabular data, and, as accurately as possible, to force the approximating function to assume exactly the value provided at each of the points where the data is supplied. For significantly “noisy” data, a single curve representing the general trend of the data is derived by the method of least-squares regression.

In statistics, Regression means finding a description of a data set. For example, if a data set fits into a normal distribution, the whole data set can then be described by two numbers: its mean and standard deviation.

The next functions (simple form) are available in SimplexNumerica:

Linear Least Squares:                 y = a + b x

Exponential function:                   y = a e^bx; a > 0

Logarithmic function:                   y = a + b ln(x)

Power function:                           y = a x^b; a > 0

Power function with weighting:    y = a x^b; a > 0

Parabolic Regression:                 y = a + bx + cx²

Cubic Regression:                       y = a + bx + cx² + dx³

Polynomial 0-9ten of degree       y = a₀ + a₁ x + a₂ x² + a₃ x³ +. . . . . + a_n xⁿ

The coefficients are determined by the solution of the Gauss' Eliminationsverfahren with partial Pivotierung. 

Error messages result: 

During logarithmic computation by values for x_i < = 0;

During exponential computation for values of y_i < = 0;

During computation to a power function; x_i and y_i must be positive values. 

The results will be shown in the Report Window.

We will show you now only a few examples with the mathematical background behind.

Linear Least Squares

Consider fitting a straight line to a set of data such as those shown the previous chart. Let the data be represented by the set of n data points.

The equation of a straight line through the data is

where a₀ is the vertical intercept and a₁ is the slope of the line. If the relationship between x and y were indeed truly linear and there was no noise in the data, then the slope and intercept could be estimated such that the line passed through all of the data points. However, as can be seen from Figure 1, this is not the case. This is a discrepancy, or residual, between the true value of y and the linear approximation. This residual is denoted by e and is defined by

The approximating function, the straight line, must now be chosen such that, in some sense, the discrepancy e is minimized over the entire range of x where the approximation is to be applied.

The sense in which e is minimized is clearly a vital factor in determining the character of the approximation. We could minimize the maximum value of e for all data points. However, this criterion is usually not an effective one to use in selecting a continuous functional approximation of noisy data, simply because it permits individual points which may be badly in error to exert overpowering influence on the approximating function. That is, a single point can force the approximating function to shift drastically toward it in order to minimize the maximum error which would tend to occur at that point. A much more favourable condition to apply to minimize e for this type of approximation is the least-squares criterion.

If we denote the x coordinates at which data are available as xi, then the i-th residual ei is

and if there are n such coordinates, then the sum of squared residuals over all the data points is

In order to determine the values of the coefficients a₀ and a₁, we can minimize S_r. The Minimization is accomplished by setting the partial derivatives of S_r with respect to each

Coefficient equal to zero:

or, dividing by -2 and summing term by term, we have

Similarly, the second equation is

or, dividing by -2 and summing term by term, we have

Now, Equations (1) and (2) represent of two simultaneous linear equations in two unknowns (a0 and a1):

These are called the normal equations. The solution to these equations is

Example:

Given the following noisy data, fit a straight line to this data by using least squares.

Each element of these equations can now be computed

The solution of the normal equations is

The required straight line is thus

Polynomial Regression

Previously, we fit a straight line to noisy data

using the least-squares criterion. As we have seen, some data are poorly represented by a straight line and for these cases a curve is better suited to fit the data. The most commonly used function for this purpose is the polynomial such as a parabola

or a cubic

or in general an mth degree polynomial:

where a₀, a₁, a₂, … , a_m-1, a_m are the constant coefficients of the polynomial.

If the relationship between x and y were indeed truly m-th degree polynomial and there was no noise in the data, then the coefficients could be estimated such that the polynomial passed through all of the data points. However, this is hardly ever the case. As in the linear case, the discrepancy (residual) between the true value of y and the polynomial approximation is

and if there are n such pairs of points (xi, yi), then the sum of squared residuals over all the data points is

In order to determine the values of the coefficients a₀, a₁, a₂, … , a_m-1, a_m, we can minimize S_r . The minimization is accomplished by setting the partial derivatives of S_r with respect to each coefficient equal to zero:

Now, dividing by -2 and summing term by term, we have

Similarly, the second equation is

Dividing by -2 and summing term by term, we have

It can now be inferred that the complete set of simultaneous linear equations in the coefficients (the normal equations) of the polynomial is

Example:

Given the following data, choose the most suitable low order polynomial and fit it to this data using the least-squares criterion.

The data are plotted in the following Figure. The data appear to have a maximum near x

= 5 and a minimum near x = 15. The lowest order polynomial which can reproduce such

behaviour is a cubic. The least-squares equations (normal equations) for this set of data (n = 24) and for m = 3 are

Thus the equations for the interpolating polynomial is

Linearization of Nonlinear Relationships

In order to apply the techniques of linear least-squares regression, the function whose coefficients are being approximated must be linear in the coefficients. Many relationships among independent and dependent variables in engineering are not linear. However, in many cases a transformation can be applied to the relationships to render them linear in the coefficients. Consider an exponential relationship,

y(x) = a e^bx

where the base is the number e , and a and b are constants. If we take the natural logarithm of both sides of the equation, we have

which is a linear relationship between ln(y) and x. The coefficients to be determined in this

expression are ln(a) and b. A power law relationship be written as

If we take the natural logarithm of both sides of this equation, we have

which is a linear relationship between ln(y) and ln(x). Again, the coefficients to be determined in this expression are ln(a) and b.

Similar to curve fitting, interpolation methods attempt to find a function (or often a set of functions) that represents a set of discrete data points (SampleData).  However, in this case, the data points are considered to be very accurate, such that the approximating relationship must not only pass exactly through each point, but also allow accurate interpolation between them. These methods are most commonly used as a means of obtaining a continuous representation of accurately measured experimental data (e.g., equilibrium data) for subsequent computational purposes.  They are also particularly useful in situations where a tedious numerical method has been employed to generate a fairly sparse set of discrete information, and this data must then be used in continuous form as a starting point for some other type of analysis on the computer.  For example, when a mathematical function like the sine function is evaluated on the computer, it is not determined from scratch using, e.g., its infinite series representation.  Rather, a call to sin(x) invokes an interpolation method that determines functional values at positions intermediate to a small number of very accurately determined values of the sine function.  Interpolation methods are also used extensively in graphical routines, to produce smooth curves connecting a set of discrete points.

Thus the general point of view in this chapter can be stated as follows:

Given a set of “accurate” data points, (x_i, y_i), i = 1, 2 , ×××, n, construct a function f(x) such that f(x_i) = y_i for each i and f(x) assumes “reasonable” values for all x between the data points.

One approach to the problem of approximating a set of data according to the above criteria is to simply find the (n-1)-order polynomial that passes through all n points.  This approach certainly meets the first condition that f(x_i) = y_i, but, in most cases, because higher-order polynomials tend to oscillate, it usually won’t meet the second.  Thus, most interpolating methods attempt to find many separate functions that fit only portions of the data, often just the region between each pair of neighbouring x values.  Furthermore, these interpolating functions are usually lower-order polynomials, e.g., quadratics, cubics, quartics, etc.

There are two different approaches that can be taken to find a lower-order polynomial which passes through two neighbouring data points and allows accurate interpolation at all intervening x values.  One method is to determine the polynomial that passes exactly through several points surrounding the interpolation interval.  This approach is illustrated below, where 3^rd-order polynomials are passed through the 4 points (2 on either side) surrounding each interval. For instances, this is the approach adopted by the Lagrange Interpolation Method.

In the second approach, a low order polynomial is found that passes only between the two points bracketing the interpolation interval, and satisfies conditions of continuity and smoothness at the two end points.  That is, not only must the value of two adjoining functions be the same at a given point, but also their various derivatives must match. The Method of Splines is the interpolation technique that adopts this second approach.  The divided difference (initiated by Newton) and the Lagrange interpolation methods have largely been superseded by the method of Splines in today’s computing environment.  However, because the former methods provide very convenient ways to determine the coefficients of the m-1 order polynomial that passes exactly through discrete points, a task that must be carried out in integration and differential equation solving methods, SimplexNumerica will examine the best interpolation techniques.

	y

 

Note that these methods become equivalent if the interpolating polynomial has only 2 terms, since there is only one straight line that can pass through a given pair of points.

Polygonal Curve

The simple polygonal function connects the individual measuring points by straight lines (Pulldown-menu Display à Polygonal Lines). Here in this menu, these straight lines are also present as CurveData. The individual measuring points are not only simply connected; in each subinterval a linear curve is computed. Thus one can integrate or differentiate a polygonal curve.

For a given set of data containing n points, there are n-1 intervals and hence n-1 separate Polygonal Curve that must be determined.  Assume that each pair of neighbouring points, (x_i, y_i) and (x_i+1, y_i+1), is joined by a linear Polygonal Curve (i.e., a 1^st-order polynomial) of the form:

where p_i and q_i (i = 1, 2, ···, n-1) are 2n-2 constants that must be determined from the following conditions:

,             i = 1, 2, ···, n-1                                               n-1 equations

,       i = 1, 2, ···, n-1                                              n-1 equations

2n-2 equations

 The first condition yields:

                    ,               \      

The second condition gives:

                    \      

Additive Segmentation

The procedure for the curve with additive segmenting is primarily for series of measurements with a large number of measuring data. The number of SampleData should be more than 50, so that one gets a smooth course of the curve.

The algorithm of the Additive Segmentation contains a calculation specification, with their assistance the progression of a function is described by a finite number of values, laying on a curve; their curvature are closely as possible at the point process.

The procedure divides the existing segments (places of the measuring points) with the help of auxiliary points into smaller segments. The number of the added segments depends on the curvature of the function. Within the range of a strong curvature the points are nearer together.

Caution:

For the production of the auxiliary points the number of CurveData are adapted!

Re-establish

You can set it back with the Pulldown-menu <Functions>, menu <Surrogate>!

(n-1)-Polynomal

The interpolation algorithm (n-1)-Polynomal computes the polynomial, which runs exactly through a number of given measured points (SampleData). If n is the number of SampleData, then this is a polynomial of (n-1) of degree.

Lagrange Interpolation

The Lagrange interpolation method uses a divided difference method, i.e., a (d-1)-order polynomial is fitted to the d points that lie closest to the interpolation interval.  However, in the Lagrange method, a different form of the interpolating polynomial is selected so that its coefficients can be determined even more easily than was the case for the divided difference method.  In the Lagrange method, the (d-1)-order polynomial that passes through the points (x_i, y_i), (x_i+1, y_i+1), ···, (x_i+d-1, y_i+d-1) has the following general form:

Note that, when all of the factors are multiplied together, each r.h.s. term is, in itself, a polynomial of order d-1.  When the polynomial is written in this form, the unknown coefficients a_i, i = 1, 2, ···, d, become very easy to determine, e.g.,

When

                        \        

Similarly,              

Thus, the Lagrange interpolation polynomial for the argument  becomes

where:                        

Note that the second form of f() is cheaper to compute; but, if , a divide by zero will occur.  However, it should be remembered that, if , then interpolation is not really necessary because, automatically, .

Newton Interpolation

Assume that a polynomial having d terms is to pass exactly through d points (x_i, y_i), (x_i+1, y_i+1), ×××, (x_i+d-1, y_i+d-1), chosen such that the interpolation argument, , is as close as possible to the centre of the points.

We will discuss in a moment how i is to be determined for a given value of d (this of course depends on how close  is to one end or other of the set of points).  The question now is how do we determine the d coefficients of the (d-1)-order polynomial that passes through these d points.  One way of determining these coefficients is to simply solve the set of simultaneous equations that result when f(x_i) is set equal to y_i for each of the d points.  However, this numerical task is not conveniently carried out when a computer is not available, which was the case in Newton’s time.  To allow easier determination of the needed coefficients, Newton started off by writing the interpolation polynomial as follows:

Note that, in this polynomial, the highest power of x is d-1, i.e., it is a polynomial of order d-1.  The advantage of writing it in this fashion is that it now becomes a straightforward task to analytically determine the coefficients, a_i, i = 1, 2, ×××, d.

When   ,    

          \      

where  is the symbol used to represent what Newton called the zeroth divided difference.

 When   ,    

When   ,    

Similarly,

Thus, Newton’s divided difference interpolation polynomial of order d-1 becomes:

where

             and

In order to use this formula for interpolation, it is first necessary to construct a matrix of divided difference values (the so-called difference table) where the first column contains the zeroth divided differences, the second column the first divided differences, etc. Each interpolation is then carried out by finding 2 neighbouring x values that lie on either side of the interpolation argument , choosing a set of points x_i to x_i+d-1 such that  lies as close as possible to the centre of that interval, and then applying the formula.  In order to illustrate how the method works, let’s consider the following example where we are given values of sin(x) at 6 equally spaced points between x = 0 and x = p/2, and we wish to find the interpolated values at  = p/4 and  = 9p/20 when using a 3^rd-order interpolating polynomial (i.e., d = 4).

The Spline interpolation used a method in which the functional behaviour between each pair of neighbouring data points (x_i, y_i) and (x_i+1, y_i+1) is represented by a lower-order polynomial of the form

where the a_ki (k = 1, 2, ···, m; i = 1, 2, ···, n-1) are constants determined by satisfying conditions of both continuity and smoothness at the end points of the interval, i.e., at x_i and x_i+1. For example, if m = 2, the splines are called Linear Splines (àPolygonal Curve), and the constants a_1i  and a_2i are determined by ensuring that the straight line segment passes through (x_i, y_i) and (x_i+1, y_i+1).  Note that all three methods, divided differences, Lagrange interpolation and splines yield identical results when a first-order interpolating polynomial is employed. The most commonly used splines are Cubic Splines, obtained by setting m = 4.  In this case, the a_ki coefficients are determined by ensuring that the value, the first derivative and second derivative of each neighbouring spline match at every interior point. 

Cubic Splines

For the cubic spline interpolation of the set of data points, (x_i, y_i), i = 1, 2, ···, n, it is assumed that each pair of neighbouring points, (x_i, y_i) and (x_i+1, y_i+1), is joined by a 3^rd-order polynomial of the form:

where the p_i, q_i, r_i and s_i are 4n-4 coefficients that must be determined from the following conditions:

1.              ,                      i = 1, 2, ···, n-1                          n-1 equations

2.              ,                 i = 1, 2, ···, n-1                         n-1 equations

3.              ,         i = 1, 2, ···, n-2                         n-2 equations

4.              ,         i = 1, 2, ···, n-2                         n-2 equations

4n-6 equations

5.              Plus one condition at x₁ and one condition at x_n that must be specified arbitrarily.

Note that, since the derivatives of neighbouring splines can only be matched at the interior points (x_i, i = 2, 3, ···, n-1), then, for cubic splines, two additional conditions must be specified in order that all 4n-4 coefficients may be determined.  In fact, there are 3 different types of arbitrary choices that are standardly made, and each yields a different set of cubic interpolating formulae.  We will examine each case and discuss their relative merits later.  For the moment, let’s see what information can be derived from the first 4 conditions.

 Condition 1:

                                                                                               (i)

Condition 2:

or                                                                           (ii)

where         

Now            

and             

Condition 3:

or                                                                                 (iii)

Condition 4:

or                

or                                                                                                (iv)

Substituting Eq. (iv) into Eq. (iii) yields:

or                               

Also,                                                                           (v)

Substituting Eq. (iv) into Eq. (ii) gives:

or                                                                           (vi)

Also,                                                      (vii)

Finally, substituting Eqs. (vii) and (vi) into Eq. (v) yields:

or                         (viii)

Note that an equation of this form can be written for i = 2, 3, ···, n-1, and consequently forms part of a tridiagonal set of equations that can be solved to determine the unknown r_i values.  To use the Thomas algorithm to solve this set, where each equation has the general form:

then we must set:           

Once the r_i values have been determined, the q_i and s_i coefficients can be calculated from Eqs. (vi) and (iv), respectively.

Approximation or Fitting means finding a smooth curve that describes the underlying pattern without necessary to touch it data points.

Gauss/Cholesky

The function Gauss/Cholesky computes the coefficients of a fit polynomial of n^th degree with the discrete least squares method from Gauss. The simultaneous equation will be solved with the method of Cholesky.

Put in the degree of the Fit polynomial in the next dialogbox.

Please enter now the priority factors, either for each individual measurement point or global for all.

The individual Buttons means:

Local

Input of the smoothing parameters separately in each case for two neighbouring sample points.  

Global

Input of the smoothing parameters for all sample points at the same time.

Load/Save

Save or load the parameter array on disk.

Afterwards the function will be calculated and then the curve is plotting.

If you activated the Report Window, then the results of the computation are noted there.

When solving systems of equations, we want Ax = b. Thus, we define the residual as r = b-Ax, and note that the residual gives us some measure of the error. Of course the goal is to obtain r=0 if possible.

Residuals are differences between the observed values and the corresponding values that are predicted by the model and thus they represent the variance that is not explained by the model. The better the fit of the model, the smaller the values of residuals.  

If a matrix is both symmetric and positive definite (which happens quite bit), we can obtain a Cholesky factorization of A =.

Suppose that we write down in order to represent A= for a

2x2 matrix. Then by multiplying out the right hand side we obtain

.

Thus we can solve for , , and .

a.    Then general algorithm is as follows:

                                                   i.           for(j=1,n) {

                                                  ii.                 for(k=1,j-1) for(i=j,n) ;

                                                iii.                  ; for(k=j+1,n) }

b.    In words this is the same as

                                                   i.          For each column j of the matrix.

                                                  ii.          Loop over all previous columns k, and subtract a multiple of column k from the current column j.

                                                iii.          Take the square root of the diagonal entry, and scale column j by that value.

c.     One important thing to notice here is that the algorithm above factors the matrix “in place”.  This dramatically reduces the storage necessary for the matrix factorization.   It is not necessary to store both A and L.  Rather, the entries of L overwrite the entries of A when they are no longer needed.  The linear systems that arise in practice are often enormous and this saving in storage cost is often critical.

d.    In the first class we discussed different types of errors (e.g., round off, truncation, modelling, etc.).  It is very important to be aware of what leading sources of error are when using numerical algorithms.  For example, in a seminal cloth simulation paper for computer graphics by Baraff and Witkin they encountered a matrix equation, Ax=b, which was difficult to solve.  They simply discarded the non-symmetric part of A and symmetrized it reducing their problem to one which can use more effective techniques like Cholesky decomposition.

2.    Both Gaussian Elimination and the Cholesky factorization method are examples of direct methods. The other kind of method is called an iterative method where one starts with an initial guess  and iterates through a sequence, ,  … ending up with a final guess of .

a.    Issues here include both finding an initial guess  and deciding on a stopping criterion that says that , for some m, is good enough.

b.    One important thing to mention here is the use of a preconditioner. The “good enough” final guess x_m can be obtained with dramatically fewer iterations with an appropriate preconditioner.  Basically we multiply both sides of Ax=b by a matrix P such that P*A~I, and then solve PAx=Pb. Preconditioners can be helpful in direct methods too, but are more often seen in iterative solvers. 

3.    So far we have discussed solving Ax = b for square  matrices A. For more general  matrices, there are a variety of scenarios. 

a.    When m < n, the problem is underdetermined since there is not enough information to determine a unique solution for all the variables. Usually m<n implies that there are infinite solutions. However, in some cases, there may be contradictory equations leading to the absence of any solutions.

b.    When m > n, the problem is overdetermined although this in itself doesn’t tell us everything about the nature of the solution. For example, if enough equations are linear combinations of each other, there can still be a unique solution or infinite solutions.

                                                   i.          We can use the rank of the matrix to enumerate the possibilities. Recall that the rank of a matrix is the number of linearly independent columns that it has. Thus a  matrix has at most a rank of n.

                                                  ii.          If the rank < n some columns are linear combinations of others and we say that the matrix is rank-deficient and there may be an infinite number of solutions.

                                                iii.          On the other hand, if the rank = n, i.e. all the columns are linearly independent and we are guaranteed either a unique solution or no solution. In the case of no solution, there is the notion of the “closest fit” in a least squares sense. That is, the least squares solution finds the closest possible solution.

4.    The least squares method finds the best-fit or best approximation to the solution in the sense that the least squares solution x minimizes the norm of the residual, i.e. it minimizes.

5.    As an example, consider interpolating a given a set of m data points, , with a straight line . Here, each data point leads to a new equation, for example with three data points we obtain .

a.    Here, if one gave the exact same pair  three times, then there is really only one point and there are an infinite number of lines that go through it.

b.    If the three points all line on the same line, then that line is a unique solution to the problem.

c.     If the three points do not all lie on a line, the problem is overdetermined and there is no solution, i.e. no straight line that passes through the points. In this case, we can look for the best least squares solution, or the line that passes as close as possible to the three points.

6.    Normal equations. If the columns of A are linearly independent, has a unique solution.

a.    In fact,  is a square matrix that is symmetric & positive definite so we can use a Cholesky factorization, i.e. or.

b.    How do we know that the matrix A^TA is symmetric and positive definite?  The symmetry is immediate from (A^TA)^T = A^TA^TT=A^TA.  There are important classes of matrices that come up often in practice that have the form of A^TA. Consider any non-zero Vector x.   Then x^TA^TAx = y^Ty >=0 with equality only when y=Ax is zero, which can only occur if A does NOT have full rank.

c.     This sounds too good to be true, and it is. The condition number of  is equal to the condition number of A squared, which can lead to real problems if A is already poorly conditioned.

Simplex-Fit

The Simplex-Algorithm or shorter Simplex-Fit is a very effective procedure for a nonlinear Fit. The Simplex algorithm function performs non-linear approximation of all the parameters that have been in the formula marked to be floated.  As always in non-linear approximation, good starting guesses are essential (for the nonlinear parameters).  Therefore, before executing the fit command, it is essential that parameter estimates are entered into the Formula parser. A similar fitting strategy implemented in SimplexNumerica is the Gauß-Algorithm. 

This program originally appeared in the May 1984 issue of Byte Magazine. It was originally written in Pascal by M. S. Caceci & W. P. Cacheris at Florida State University.  This program is based upon the Simplex curve fitting algorithm. For a detailed description of this program and its workings see the above mentioned article. I acknowledge the work of Marco Caceci and William Cacheris for writing the original Pascal program from which this is derived.  The original authors explicitly stated ''no copyright''.

The Simplex-Fit is quite fast and uses no analytical derivative of the fit function. The draw-back is that there is no direct error estimation of the optimized parameters. To get an idea of the error, one would need to repeat the fit very often starting from random initial conditions and then evaluate the statistics of fit results. Although the Simplex-Fit is not always the most efficient approach, it is known to be very stable. Starting from the initial parameter estimates, a set of random variations of the parameters are generated. The sum of squared residuals is calculated for each of the variations, and the simplex algorithm performs a series of optimization steps, in each of which the worst parameter set is eliminated and a better one is introduced. The algorithm is known to be susceptible to 'circling' the optimum, i.e. it can get trapped in a path of parameter values close to the optimum without converging.  As convergence criterion, the sum of squared residuals and all of the parameters have to be within the pre-determined Tolerance value.  After convergence, the simplex is repeatedly restarted (using randomly selected variations around the previously found optimal values), until the parameters and the values are within the tolerance.  Two consecutive simplex procedures are required to converge to the same values and to the same parameter values for the fit to stop.  After the simplex has converged, SimplexNumerica will show you the results. 

The text output during the fitting procedure shows the number of simulations of the sedimentation process, and the simplex number, the simplex step # and the finite element gridsize. The remaining lines show the values and the parameter values from the current simulation. The last line shows residuals from the current simulation.

Convergence can be observed if the value does change only very little, and when the parameter values stay virtually the same.

The Simplex algorithm used four mechanisms of vertex movements. Vertices are moved toward the minimum point by the four mechanisms, reflection, expansion, contraction, and shrinkage. Reflection moves the highest (worst) vertex to the opposite side of the center of the other vertices. Expansion moves the highest vertex to twice the distance from the center, which is used when a reflected vertex is lower than the lowest vertex and the expanded vertex becomes the lowest one. Contraction moves to the middle between the highest vertex and the center, which is used when a reflected vertex is higher than the highest vertex. When the contracted vertex is still higher than the highest vertex, all the vertices except the lowest vertex move toward the lowest one by half, which is called shrinkage.  

After the call off the Simplex-Fit you can immediately stop it.

After the Next-Button, put the function equation (formula) with the associated parameters into the Simplex Formula parser.

The operation of the above formula box facilitates the input of the function equation for the user. One can arrange the formula both over the keyboard and by activating the small buttons with the mouse. Go to the Gauß-Fit in this manual for getting more info about the formula parser dialog and the right parameter order.

Die Parameter dienen vorerst als Startwerte für den Algorithmus. The program varies these parameters and tries to find a minimum of the sum of the discrepancy squares.

Next step, decide where the output goes in:

Now a command line window appears on the screen. Enter the data in old DOS manner.

Max. Number of Iteration: Maximal number of fit iterations. When this number of iterations is reached, the fit will be terminated.

Number of parameters to fit: This is normally the number of parameters in the input formula. Be carefully to setup the right amount of the number.

Start coordinate (1 .. n): The start coordinate is the start value for each individual parameter.

Max, error (1 .. n): The error value for each individual parameter.

If the algorithm found a minimum, then one says: it converges. It can be possible, that there are more than one minimum are available. Then, the algorithm does not have necessarily found the absolute minimum. As abort criteria will be the maximum number of iterations and the difference of the remaining variance between two iterations. When it fall short of these criteria, thus the Fit becomes cancelled. Sometimes it is difficult to find suitable initial values and increments to a convergence of the algorithm.

For this some references:

The initial values should be oriented at the orders of magnitude of the measured values. The increments should amount to approximately a third step of the data, whereby they can be also more largely selected at the beginning. If the algorithm a minimum, it reduces the increment independently. The results of a not handy Fit should be used as new initial values for the next try.

If you have activated the report window, then the results stand there in short form.

The adapted function is drawn afterwards in the chart.

The output on disk was saved in a file and has the following content:

*************** S I M P L E X - algorithm ***************

Date: 23.12.2006

Time: 14:15:08

Fit Funkcion: f(x)= b+f*x+g*x^2+h*x^3

Max. iteration: 50

No. of parameters: 4

Start coord.:

Start coord.(1)= -7.30000000020914

Start coord.(2)= 8.56274281291185

Start coord.(3)= -2.00058275061708

Start coord.(4)= 0.13364413364612

Start_Step(1)= 0.1

Start_Step(2)= 0.1

Start_Step(3)= 0.1

Start_Step(4)= 0.1

Error(1)= 0.0001

Error(2)= 0.0001

Error(3)= 0.0001

Error(4)= 0.0001

Error(5)= 0.0001

Exit after 50 iterations

The final simplex is:

                    -7.30000000020914

8.56274281291185

-2.00058275061708

0.13364413364612

6.29230769230795

-7.2999682539064

8.56263415934837

-2.00069978598065

0.133651400687777

6.29242318002935

-7.30014329738388

8.56251458954117

-2.00068239819608

0.133661116175654

6.29234377881553

-7.29942926051024

8.56324741922046

-2.00062851357081

0.133642830329596

6.29232577588895

    -7.30058318261815

8.56227966156391

-2.00061895183689

0.133647192955672

6.29246130424079

The estimated fractional error is:

-0.000158083881182603

                    0.000113012927126811

-0.0000585006361428703

0.000136807521751745

0.0000244120583999725

The mean is:

Simplex(1)= -7.300024798925