Exponential Curve Fitting


y = 1.5e0.3x - 2
Introduction
This article describes the exponential curve fitting
method implemented in Graphics-Explorer.
At the end of this page, the complete Delphi unit is listed.

Given a set of points (xi , yi) , (i = 1..n) in a plane,
the user may select an exponential function of type crossing these points with a,b,c calculated for minimal deviation.

Finding a,b,c requires the following steps: Initialization
Procedure "Expfunc(n)" is called. n indicates the type of function requested.
This description will focus on y = aebx + c.
In this case procedure "maakexpfunctie" is called to If there are 3 or more points, "meerpuntenexp" is called.
From table Epunten new tables "dx" , "dy" , "dxy" are built (notice, that dxy is a record) The idea is to differentiate y = aebx + c .
In the derivative y' = abebx ,c is eliminated.
We guess, that the slope dy/dx between two adjacent points is most accurate in the middle
between the x-coordinates, so x is modified accordingly.

Approximation of b
Working with the derivative we can write the system of equations: for 2 adjacent points i and j = i+1 in the dxy table: An average value of b is calculated.

Approximation of a
Starting with we may write: An avarage value of a is calculated.

Approximation of c
Starting with Again, an average value of c is calculated.

Fine tuning c
Considering between 2 adjacent points, for this value of c, the growth factor is calculated
and placed in a table (groeifactor).
From this table the average growth factor is calculated.
Then the absolute deviation of the factors is made.
If c is accurate, the deviation of growth factors from the avarage value
will be minimal, since exponential functions are the result of constant
growth factors.

The growth factor between points xi and xj is: If (yj - c) / (yi - c) <= 0, an abort is generated.
This means, that y = c crosses the curve.

The following steps take place, while generating the deviation of growth factors for each c: A "watchdog" counter w terminates the process if 100 steps are taken.
An error message is generated in this case.
If y = c gets too far away from the curve, for each step further away
less deviation between the growth factors will be found (as they are near 1),
which is erroneously seen as a better result, while c runs away.
If the exponential curve has a sharp bend, there is only a small interval
for c to find it's best value.

Finalization
Having an accurate value of c, now, temporarily
- replacing ln(yi - c) by y and
- ln(a) by a, a set of equations can be written The least-squares method is called to find the best values for a and b.
(procedure buildexp12).
The values for a and b are returned in array Epunten.
Finally, the formula string is made.
Procedure expPlot places the string in the edit-window and calls
procedure addform (outside unit) to translate and plot the function.

This concludes the description.

David E. Dirkse
Castricum
the Netherlands
david.math@hccnet.nl

Below, I present the complete Graphics-Explorer unit, that implements the exponential curve fitting.

(note: I noticed the missing of an occasional '/' character in the listing below.
An error of my conversion program for Delphi pascal units to HTML.
I believe to have these errors fixed manually.)

remark:
var scaledY : double;//is the distance between 2 adjacent pixels on the Y scale.
//it is calculated outside the listed unit

unit expfuncUnit;
{
 calculates exponential functions

 expfunc 1 : y = ab^x+c
         2 : y = ae^(bx)+c
         3 : y = e^(ax^2+bx+c)
         4: ab^x
         5: ae^(bx)
}

interface

uses sysutils,xlateUnit,puntenUnit,formUnit,scaleUnit,dialogs;

type puntarray = array[1..40] of TPunt;//punt = point

procedure expfunc(n : byte);

implementation

uses controlUnit,klkwUnit;

const FlString : string = '####0.0#####';

var Epunten : puntarray;
    a,b,c : double;       //y = a*e(bx)+c
    emacht : boolean;     //type 1,2 te berekenen
    dx,dy : array[1..40] of double;
    dxy : array[1..40] of TPunt;

procedure ExpPlot(s : string);
begin
 controlform.formedit.text := s;
 addform(formNr);         //addsformule nr formNr to screen
end;

procedure Getpunten;     //from points Unit
//sort for x small to large
var i,j : byte;
    h : Tpunt;
begin
 for i := 1 to maxpunt do
  begin
   Epunten[i].x := punten[i].x;
   Epunten[i].y := punten[i].y;
  end;
 for i := 1 to maxpunt-1 do
  for j := i+1 to maxpunt do
    if Epunten[j].x < Epunten[i].x then
     begin
      h.x := Epunten[i].x;
      h.y := Epunten[i].y;
      Epunten[i].x := Epunten[j].x;
      Epunten[i].y := Epunten[j].y;
      Epunten[j].x := h.x;
      Epunten[j].y := h.y;
     end;
end;

function Nietgenoegpunten(n : byte) : boolean;//if not-enough-points:false
begin
 if maxpunt < n then
  begin
   controlform.msgtext.caption := inttostr(n)+' dots needed';
   result := true;
  end
 else result := false;
end;

procedure tweepuntenExp;//2 points 
//calculate a and b, if just 2 points
var teken1,teken2 : boolean;
    s : string;
    i : byte;
begin
 teken1 := Epunten[1].y < 0;//teken means sign
 teken2 := Epunten[2].y < 0;
 if teken1 <> teken2 then
  begin
   controlform.msgtext.caption := 'no exp.function possible: unlike signs';
   exit;
  end;
 if teken1 then
  for i := 1 to 2 do Epunten[i].y := -Epunten[i].y;//a negatief, keer om
 for i := 1 to 2 do Epunten[i].y := ln(Epunten[i].y);
 b := (Epunten[2].y - Epunten[1].y)/(Epunten[2].x - Epunten[1].x);
 a := Epunten[1].y - b*Epunten[1].x;
 a := exp(a);
 if not emacht then b := exp(b);
 if teken1 then a := -a;
//a,b,c berekend
 s := 'y = ';
 s := s + formatfloat(flstring,a) + ' * ';
 if emacht then s := s + 'e^(' + formatfloat(Flstring,b) +'x)'
  else s := s + formatfloat(Flstring,b) + '^x';
 expPlot(s);  //invoeren en plotten
end;

function Afw(cc : double; var af : boolean) : double;
//calculate deviation in growthfactors after asympt. c correction
//called for fine tuning c in ae^(bx)+c formula
var i : byte;
    avG,ee : double;
    groeifactor : array[1..40] of double;
begin
 avG := 0;
 for i := 1 to maxpunt-1 do
  begin
   ee := (Epunten[i].y-cc);//noemer
   if ee <> 0 then ee := (Epunten[i+1].y-cc)/ee;
   if ee <= 0 then
    begin
     af := true;//abort flag
     exit;
    end;
   ee := ln(ee);
   avG := avG + ee;
   groeifactor[i] := exp(ee/dx[i]);//save
  end;
 af := false;
 avG := exp(avG /(Epunten[maxpunt].x-Epunten[1].x));//exp toegevoegd
 ee := 0;
 for i := 1 to maxpunt-1 do
  ee := ee + abs(avG-groeiFactor[i]);
  result :=  ee;
end;

procedure maakXYarrays;
var i : byte;
begin
 for i := 1 to maxpunt-1 do
  begin
   dx[i] := Epunten[i+1].x - Epunten[i].x;
   dy[i] := Epunten[i+1].y - Epunten[i].y;
   dxy[i].x := (Epunten[i].x + Epunten[i+1].x)/2;
   dxy[i].y := dy[i]/dx[i];
  end;
end;

procedure ExpAbort(m : byte);
const abortmessage : array[1..2] of string =
      ('no exp.function',
       'out of range');
begin
 controlform.msgtext.caption := abortmessage[m];
 controlform.formedit.text := '';
end;

procedure meerpuntenExp;
//called more than 2 points available
var Som,deltaC,Bvar1,Bvar2 : double;
    i : byte;
    w : word;
    s : string;
    abort : boolean;
begin
 maakXYarrays;
//find approximation of b...growthfactor e^b
 try
  som := 0;//som = sum
  for i := 1 to maxpunt-2 do
   Som := Som + ln(dxy[i+1].y/dxy[i].y)/(dxy[i+1].x-dxy[i].x);
  b := Som / (maxpunt - 2);
//find approximation of a 
  Som := 0;
  for i := 1 to maxpunt-1 do
   Som := Som + dy[i]/(exp(b*Epunten[i+1].x) - exp(b*Epunten[i].x));
  a := Som / (maxpunt - 1);
//find approximation of c 
  som := 0;
  for i := 1 to maxpunt do
   som := som + Epunten[i].y-a*exp(b*Epunten[i].x);
  c := som / maxpunt;
 except
  expAbort(1);
  exit;
 end;
 if a > 0 then deltaC := -scaledY //loop away from asymptoot at first
  else deltaC := scaledY;
 repeat
  Bvar1 := Afw(c,abort);         //get deviation in growth factors
  if abort then c := c + deltaC;
 until not abort;
 w := 0;
 deltaC := -deltaC;//first towards asymtoot
 repeat
  inc(w);
  Bvar2 := Afw(c + deltaC,abort);
  if abort then deltaC := deltaC/2  //halve distance, same direction
  else
   if Bvar2 < Bvar1 then
    begin
     Bvar1 := Bvar2;  //if better
     c := c + deltaC;
    end
   else deltaC := -deltaC/2;
 until (abs(deltaC) <= 1e-8) or (w > 100);
 if abort then
  begin
   expAbort(1);
   exit;
  end;
 if w > 100 then
  begin
   expAbort(2);
   exit;
  end;  
 //use least squares method to make best a, b
 for i := 1 to maxpunt do
  Epunten[i].y := ln(abs(Epunten[i].y-c));
 buildexp12(Epunten);
 if a < 0 then a := -exp(Epunten[1].y)
 else a := exp(Epunten[1].y);
 b := Epunten[2].y;
//make formula
 s := 'y = ' + formatfloat(Flstring,a)+' * ';
 if emacht then s := s + 'e^(' + formatfloat(Flstring,b) + 'x)'
 else begin
       b := exp(b);
       s := s + formatfloat(Flstring,b) + '^x';
      end;
 if c < 0 then begin
                c := abs(c);
                s := s + ' - ';
               end
 else s := s + ' + ';              
 s := s + formatfloat(Flstring,c);
 expPlot(s);
end;

procedure maakExpFunctie;//make exponential function
begin
 if nietgenoegpunten(2) then exit;
 GetPunten;
 if maxpunt = 2 then
  begin
   c := 0;
   tweepuntenExp;//just 2 points
   exit;
  end;
 meerpuntenexp;//more than 2 points
end;

procedure exp3;//function '3'
var i : byte;
    s : string;
begin
 if nietgenoegpunten(3) then exit;
 getPunten;
 for i := 1 to maxpunt do
  if Epunten[i].y <= 0 then
   begin
    controlform.msgtext.caption := 'error: y <= 0';
    exit;
   end
   else Epunten[i].y := ln(Epunten[i].y);
 buildexp3(Epunten); //call klkw unit,return a,b,c at Epunten[1..3].y
 s := 'y = e^(';
 s := s + formatfloat(Flstring,Epunten[3].y) +'x^2';
 if Epunten[2].y < 0 then
  begin
   s := s + ' - ';
   Epunten[2].y := abs(Epunten[2].y);
  end
  else s := s + ' + ';
 s := s + formatfloat(Flstring,Epunten[2].y) +'x';
 if Epunten[1].y < 0 then
  begin
   s := s + ' - ';
   Epunten[1].y := abs(Epunten[1].y);
  end
  else s := s + ' + ';
 s := s + formatfloat('####0.0#####',Epunten[1].y);
 s := s + ')';
 ExpPlot(s);
end;

procedure exp45;//functions '4,5'
var i : byte;
    s : string;
begin
 if nietgenoegpunten(2) then exit;
 getPunten;
 for i := 1 to maxpunt do
  if Epunten[i].y <= 0 then
   begin
    controlform.msgtext.caption := 'error: y <= 0';
    exit;
   end
   else Epunten[i].y := ln(Epunten[i].y);
 buildexp12(Epunten);
 a := exp(Epunten[1].y);
 b := Epunten[2].y;
 if not emacht then b := exp(b);
//a,b,c berekend
 s := 'y = ';
 s := s + formatfloat(flstring,a) + ' * ';
 if emacht then s := s + 'e^(' + formatfloat(Flstring,b) +'x)'
  else s := s + formatfloat(Flstring,b) + '^x';
 expPlot(s);  //invoeren en plotten
end;

procedure expFunc(n : byte);
//n 1..5 : type of function requested
begin
 case n of
  1 : begin
       emacht := false;//no power of e
       maakexpfunctie;
      end;
  2 : begin
       emacht := true;//power of e
       maakexpfunctie;
      end;
  3 : exp3;
  4 : begin
       emacht := false;
       exp45;
      end;
  5 : begin
       emacht := true;
       exp45;
      end;
 end;//case
end;

end.