Exponential Curve Fitting

download exerciser program
download complete Delphi-7 project
see source code of curve fitting


This new article describes the exponential curve fitting method implemented in Graphics-Explorer,
my equations grapher program.
It replaces the old article, which can be found [here].
New is an exerciser program allowing step by step observation of the curve fitting process.

The curve fitter calculates the best fitting exponential function given a set of points.
This function is
    y = a.bx + c
where a,b,c are called the parameters.


Exerciser description

Below is a reduced picture of the exerciser at work
    a,b,ccalculated parameters. (black:valid; red:invalid)
    preset paramsclick radio button to preset parameters
    preset x1..x5click radio button to preset points
    x, ylist of points, maximum is 10
    [OK] validateclick to check/organize points prior to calculation
    [GO]start calculation of parameters a,b,c
    [ ] testmodemark to watch step by step calculation


Points can be added by:
    1. selecting 1 of 4 preset point sets (click radio button)
    2. mouseclick on the coordinate system (2nd click removes point)
    3. typing into the [x,y] points list
After changing points in any way, click the [OK] button to check and to sort the points.
To watch the process step by step, mark the [testmode] checkbox.
Then click the [GO] button to calculate the a,b,c parameters and paint the graph.
To clear parameters or points, click the appropriate [X] button or [RESTART].

Curve fitting theory

Input to the curve fitter is a set of points [x1,y1]....[xn,yn]
The minimal required number of points is 3.
The maximum number of points is 10.
The fitter calculates parameters a,b,c such that the curve y = a.bx + c has the smallest distance to these points.

For internal computations y = a.ebx + c is used (eb is reported for b)
The method has a numerical and an analytical part.
These are the steps:
    1. obtain an approximation of b
    2. using b, obtain an approximation of a
    3. using a,b, obtain an approximation of c
    4. discard a,b and replace y = aebx + c ....by.... y - c = aebx
    5. ln(y - c) = ln(a) + bx ...........a linear function
    6. call the least square method to find a,b of the best fitting line

approximation of b

array dx holds the x differences between adjacent points
dx1 = x2 - x1....dx3 = x4 - x3

array dy holds the y differences
dy1 = y2 - y1....dy3 = y4 - y3

array cx holds the centers between points
cx1 = (x1+x2)/2....cx3 = (x3+x4)/2

array dq holds the differential quotients
dq1 = dy1/dx1....dq3 = dy3/dx3

Note: in the Delphi language, cx1 is written as cx[1]...etc.

Starting with
    y = aebx + c..........differentiation
    y' = abebx
which eliminates c.
So, we may write
    y1' = abebx1
    y2' = abebx2
    y3' = abebx3
Parameter a can be eliminated by division of subsequent y':

with above arrays in mind:

The center values (cx) are used, because the tangent (differential quotient dq) is most accurate
between two x values.
For subsequent points, the average value of b is calculated.
This is the approximation of b.

approximation of a

Starting with
    y1 = aebx1 + c
    y2 = aebx2 + c
we may eliminate c by subtraction of subsequent rows:
    dy1 = y2 - y1 = a(ebx2 - ebx1)

An average value of a is calculated.
This is the approximation of a.

approximation of c

    y = aebx + c .....so
    c = y - aebx
We may write
    c = y1 - aebx1
    c = y2 - aebx2
An avarage value of c is calculated.
This is the approximation of c.


In equation
    y = gx
the value of y will increase by a factor g if x increases by 1.
g is called the growth factor.

Exponential functions (having c=0, the x-axis is the horizontal asymptote), are the result
of constant relative growth.

In equation
    y = ebx
we notice that g = eb.
b is called the relative growth speed because

So, if a function has a constant relative growth speed of b, it's growth factor g = eb
where e is the base of the natural logarithm. (e = 2.71828....)

Before, we have obtained an approximation of c so we may write ...y - c = aebx
If the points (x1,y1)...(xn,yn) are on a true exponential function,
we will find the same growth factor, comparing subsequent points.
In the next part, fine tuning of c, we are going to vary c to obtain the smallest
deviation in growth factors.

fine tuning c

At this point we forget the a and b values and continue with c and the points [x1,y1]...[xn,yn].

For a certain value of c we write
    y1-c = aebx1
    y2-c = aebx2
division eliminates a

Two successive points produce a value of ee.
The values of ee are summed in variable AVG, for now AVG = ee1 + ee2 + ....
The growthfactors are stored in array growthfactor:
growthfactor1 = eee/dx1..........

By summing the ee values, we actually sum values b.dx1 + b.dx2 + b.dx3... = b(last X - first X)
So, b = AVG / (last X - first X)
Average growth factor is eb
    AVG = eAVG/(lastX-firstX)
Next, we have to compare the values in the growthfactor array with AVG.
Deviation = abs(AVG-growthFactor1) + abs(AVG-growthfactor2) + ......
Deviation is the result of procedure deviation which is called for a certain value of c.

least square method

The least-square-method is an algorithm to calculate the best fitting polynomial y = c0 + c1x + c2x2+ ....
given a set of points.
With a final value of c, we may write

which represents a linear function like y = a + bx.
We present the points (ln(yi-c),xi) to the least square procedures to obtain the best a,b values.
Resulting a is replaced by ea to cancel the ln( ) function.

The least square method is a nice application of linear algebra.
Look [here] for a description.


The deviation procedure produces an abort if c is between two y values.
In that case y = c intersects the curve.
c must be outside the range y1, y2,y3.....
    ee = (y2 - c)/(y1 - c)
    abort = ee <= 0
A negative or zero ee value does not allow for an ln function.
Below is the flowchart for the fine tuning of c
DeltaC is the increment value of c.
The fine tuning ends when deltaC has become very small.

Watchdog timer

Procedure morepointsexp controls the calculation of the a,b,c parameters.
Floating point operation are between try...except directives.
There are 2 cases for an abort:
1.Illegal Floating point operation
This situation occurs if the points are too far off to make an exponential function
an example being (1,1) , (5,5) , (8,2)

2.C runaway
If C is too far out, an even further position may yield a smaller growth factor deviation.
So, C "runs away".
The watchdog counter w ends this process.
Fine tuning of c must be completed in maximal 100 steps, otherwise the process is aborted.
This situation occurs in case of sharp bends.

The Delphi-7 project

There are 4 units and 1 form:
    form1,unit1: paintbox, events, controls
    unit2: data and procedures for the curve fitting algorithm
    unit3: procedures for testmode reporting
    unit4: least square procedures


Graphs are painted in bitmap map, points are painted on paintbox1.
The dimension is 800*800 pixels, coordinate (0,0) is at map position (400,400).
x domain is -10< x < 10,.........y range is -10< y <10.
The scale is 40:1.
There are 2 important flags:
    pointsOK : all points in the list are OK and sorted in ascending order.
    paramsOK : parameters are valid
If paramsOK = true, pressing [GO] just plots the function.
If paramsOK = false, pressing [GO] calculates a,b,c and then plots the function.
if parameters a,b,c do not represent the actual situation, they are painted in red.

Please refer to the source code for details.


Main call is made to procedure makeExpFunction.
Points are collected from unit1 and stored in array Epoints[ ] of type Tpnt.
type TPnt = record
             x,y : double;

procedure deviation calculates the variation of growthfactors while fine tuning c.
These values are stored in variables Bvar1,Bvar2.


Unit2 calls procecures in unit3 to display the steps in the curve fitting process.
Information is written in form1 paintbox1.
The waitflag : boolean (in unit1) controls the step-by-step action.
While the waitflag is true, the program executes the processmessages procedure.


Has the necessary procedures for the least square method.
Main call is to procedure buildexp12(var pp : Tpoints; maxp : byte),
where Tpoints is the array of points and maxp the number of points supplied.

This concludes the description.