ProfitPress Mega CDROM2 …eeware (MSDOS)(1992)(Eng)/ProfitPress-MegaCDROM2.B6I/MISC/MATH/DELP.ZIP/DELP.PAS

Pascal/Delphi Source File | 1991-09-30 | 13.9 KB | 334 lines

{**************************************************************************}
{*                                                                        *}
{*                @@@@@     @@@@@@  @        @@@@@                        *}
{*                @    @    @       @        @    @                       *}
{*                @     @   @@@@    @        @    @                       *}
{*                @    @    @       @        @@@@@                        *}
{*                @@@@@     @@@@@@  @@@@@@   @                            *}
{*                                                                        *}
{ Differential Equation Legendre Polynomial - method of determining the    }
{ best fit linear homogeneous differential equation with constant coeffi-  }
{ cients,  by solving explicitly for the coefficients of its characteristic}
{ equation. It is demonstrated here to find the (real) rate constants of a }
{ function that is a sum of exponentials.                                  }
{ MUCH ABRIDGED DEMONSTRATION VERSION-full version on request to DJM       }
{**************************************************************************}

{ Devised by Dr J Martin, Dept Physics, Kings College, Strand,             }
{ LONDON WC2R 2LS,  and implemented by DJ Maconochie, Washington University}
{ BOX 8054, 660 S Euclid, St Louis MO 63110 (to whom correspondence may be }
{ initially addressed.                                                     }
{ E-mail [e-mail]                                            }
{                                                                          }
{ Feel free to use the procedures in this program, but please reference    }
{ Martin and Maconochie (1989) J Physiol 418:9P7 in any publications that  }
{ analysed data using this algorithm.                                      }
{                                                                          }
{ The following are available on application to DJM  a full mathematical   }
{ treatment of method in the form of a paper,  a data analysis suite       }
{ writen in Borland PASCAL Turbo 5 that incorporates this method in an op- }
{ timised form,  (limited) free advice on how best to use this algorithm   }
{to fit sums  of closely spaced exponentials to your data.                 }
{**************************************************************************}

PROGRAM DELP;

{**} Uses crt, graph, mygraph;  { Routines specific to Turbo Pascal or DJM }

CONST maxData             = 2000;        { Maximum number of data points   }
      maxSeriesSize       = 50;          { Maximum size of Legendre series }
      maxIntegrationLevel = 8;           { Max. number exponential cmpts.  }

TYPE  Tdata          = array[0..maxData] of real;
      TdataPTR       = ^Tdata;
      Tcomplex       = RECORD Re : real; Im : real; END;


VAR   LegendreSeries   : array[0..maxSeriesSize] of real;
      SeriesSize       : integer;
      IntegratedSeries : array[0..maxIntegrationLevel] of
                         array[0..maxSeriesSize] of real;
      IntegrationLevel : integer;

      CalcMatrix      : array[0..maxIntegrationLevel] of
                         array[0..maxIntegrationLevel] of real;
      Exponents        : array[1..maxIntegrationLevel] of Tcomplex;

      dataPTR          : TdataPTR;  { Points to memory containing data     }
      solutionPTR      : TdataPTR;  { Memory used to store polynomial fit  }
      StartFit         : integer;   { Start using data from here           }
      EndFit           : integer;
      FitPoints        : integer;

      i, j             : integer;

PROCEDURE even_number_points; { Data "x" range must be symmetrical about 0 }
begin
  FitPoints := EndFit-StartFit;
  if (FitPoints mod 2 = 1) then
    EndFit := EndFit-1;
  FitPoints := EndFit-StartFit;
end;

PROCEDURE dummy_data(L1, L2, G : real); { Create some artificial data: two }
var m : integer;                     { exponential components plus noise   }
begin
                   { First the noise, partly correlated from point to point}
  dataPTR^[startFit] := 0;
  for m := StartFit+1 to EndFit do
             dataPTR^[m] :=
               0.1*dataPTR^[m-1] + G*exp(random*2)*cos(random*PI);
                                        { Now add the exponential function }
  for m := StartFit to EndFit do
    dataPTR^[m] := dataPTR^[m] + 5000*exp(-L1*2*(m/FitPoints))
           - 5000*exp(-L2*2*(m/FitPoints));
end; { of procedure dummy_data }

PROCEDURE Get_Legendre_series{(SeriesSize : integer)}; { Fits the first    }
{ SeriesSize members of the Legendre series of polynomial functions to the }
{ data.  Considerable economy of computational time may be achieved by     }
{ declaring a number of local variables, and doing many of the calculations}
{ only once.                                                               }
var r,m : integer;
    x,s : real;
    p : array[0..maxSeriesSize] of real;

begin
  for r := 0 to maxSeriesSize do
    LegendreSeries[r] := 0;
  s := 8/3;
  for m := startFit to endFit do begin    { Integrate over all data points }
    x := -1+2*(m-startFit)/(FitPoints);   { Normalise "x" to -1 to +1      }
    s := 4-s;                             { Toggles between 8/3 and 12/3   }
    p[0] := s/FitPoints;                  { First two Legendre polynomials }
    p[1] := x*p[0];                       { must be given explicitly       }

    if (m=startFit) or (m=endFit) then begin   { These two points are      }
      p[0] := 2/(3*FitPoints);   { special: Simpson's rule for integration }
      p[1] := x*p[0];
    end;

    for r := 2 to SeriesSize  do         { Find the other polynomials from }
      p[r] := ((2*r-1)*x*p[r-1]-(r-1)*p[r-2])/r; { this recurrence relation}

{ Standard formula for finding the expansion of F as an infinite orthogonal}
{ series.  (This is similar to taking a Fourier transform)                 }
    for r := 0 to SeriesSize do
      LegendreSeries[r] := LegendreSeries[r]+         {      ··⌠+1····     }
                           (r+0.5)*p[r]*(dataPTR^[m]) {  (r+½).│···F.Pr.dx·}
  end;                                                {    ····⌡-1····     }
END; { of procedure Get_Legendre_Series }

PROCEDURE Legendre_fit{(SeriesSize : integer)};
{ This procedure creates  a curve corresponding to the truncated Legendre  }
{ series, to overlay the data. It is stored in the memory at solutionPTR.  }

var r,m : integer;
    f,x,s : real;
    p : array[0..maxSeriesSize] of real;

begin

  p[0] := 1;                             { The first Legendre polynomial   }
  for m := startFit to endFit do begin

    x := -1+2*(m-startFit)/FitPoints;    { Normalise "x" range to -1 to +1 }
    p[1] := x;                           { The second Legendre polynomial  }

    for r := 2 to SeriesSize do          { Find the other polynomials from }
      p[r] := ((2*r-1)*x*p[r-1]-(r-1)*p[r-2])/r; { recurrence relation     }

    f := 0;
    for r := 0 to SeriesSize do     { Add up the contribution, at each "x" }
      f := f+LegendreSeries[r]*p[r];{ of each  Legendre coefficent p(r)    }

    solutionPTR^[m] := f;           { Write to memory location solutionPTR }
  end;
END;{ of procedure Legendre_fit }

PROCEDURE Integrate_Legendre_Series; { The integral of a Legendre series   }
{ another Legendre series. The coefficients of the integrated series are   }
{ related by a recurrence relation to the un-integrated series.            }
var l,r : integer;

begin
  for r := 0 to maxSeriesSize do            { Zero the matrix containing   }
    for l := 0 to maxIntegrationLevel do    { the integrated series        }
      IntegratedSeries[l,r] := 0;

  for r := 0 to SeriesSize do             { enter the un-integrated series }
    IntegratedSeries[0,r] := LegendreSeries[r];

  for l := 1 to IntegrationLevel do   { integrate "IntegrationLevel" times }
    for r := l to SeriesSize - l do
      IntegratedSeries[l,r] := IntegratedSeries[l-1,r-1]/(2*r-1)
                             - IntegratedSeries[l-1,r+1]/(2*r+3);

END;{ of procedure Integrate_Legendre_Series }

PROCEDURE accumulate_matrix(k : integer);
{ This procedure builds up the matrix representing "IntegrationLevel"      }
{ simultaneous equations, whose solution is the set of coefficients of the }
{ auxiliary equation of the Integral/Differential equation, that is the    }
{ best least squares fit to the data.                                      }

var i,j : integer;

begin

  for i := 0 to IntegrationLevel do
    for j := 0 to IntegrationLevel do
      CalcMatrix[i,j] := CalcMatrix[i,j]
            +IntegratedSeries[i,k]*IntegratedSeries[j,k]/(k+0.5);

END;{ of procedure accumulate_matrix }

PROCEDURE Solve(n : integer);
{ This procedure solve a set of "n" simultaneous equations decribed by the }
{ matrix "CalcMatrix"                                                      }

var i,j,k        : integer;
    temp         : real;

begin
    j := n;
    for j := n downto 1 do begin
      i := j-1;
      for i := j-1 downto 0 do begin
        temp := CalcMatrix[i,j]/CalcMatrix[j,j];
        if j=1 then
          CalcMatrix[i,0] := CalcMatrix[i,0]-CalcMatrix[j,0]*temp
        else
          for k := 0 to j-1 do
            CalcMatrix[i,k] := CalcMatrix[i,k]-CalcMatrix[j,k]*temp;
      end;
    end;

    for i := 1 to n do begin
      for j := 0 to i-1 do
        CalcMatrix[i,j] := CalcMatrix[i,j]/CalcMatrix[i,i];
    end;

    for j := 1 to n-1 do begin
      for i := j+1 to n do
        CalcMatrix[i,0] := CalcMatrix[i,0]-CalcMatrix[j,0]*CalcMatrix[i,j];
    end;
{ The solutions are: 1, and the coefficients "CalcMatrix[i,0]"            }
END;{ of procedure solve }

PROCEDURE Find_Integral_Equation;

var l, i, j : integer;

begin
  for i := 0 to maxIntegrationLevel do
    for j := 0 to maxIntegrationLevel do begin
      CalcMatrix[i,j] := 0;
    end;

  for l := IntegrationLevel +1 to SeriesSize - IntegrationLevel do
    accumulate_matrix(l);        { Get together the simultaneous equations }
  solve(IntegrationLevel);       { Solve them.                             }
END; { of procedure Find_Integral_Equation }

PROCEDURE find_roots; { This routine is only for demonstration of real     }
var a, b, ac : real;  { single and double exponential fits. For anything   }
                      { else, you will need a routine for finding the      }
                      { roots of a general polynomial.                     }
begin
  if IntegrationLevel = 1 then begin
    Exponents[1].Re := CalcMatrix[1,0];
    Exponents[1].Im := 0;
    writeln('Single root only:',Exponents[1].Re:6:2);
  end;
  if IntegrationLevel = 2 then begin
    b  := CalcMatrix[1,0];
    a  := -1;
    ac := a*CalcMatrix[2,0];
    Exponents[1].Re := -0.5*(-b+sqrt(sqr(b)-4*ac));
    Exponents[2].Re := -0.5*(-b-sqrt(sqr(b)-4*ac));
    Exponents[1].Im := 0;
    Exponents[2].Im := 0;
    writeln('Real roots only, L1:',Exponents[1].Re:6:2,
                           '  L2:',Exponents[2].Re:6:2);
  end;
  if
  IntegrationLevel > 2 then begin
    writeln('You need a routine for finding the "n" roots of an equation');
    writeln('this one substitutes a Borland Turbo Pascal routine');
  end;
END; { of procedure find_roots. Don't forget to replace this with a more   }
{ general numerical routine.                                               }

BEGIN; { main program body }
  randomize;
  New(dataPTR);           { Get memory allocations at "dataPTR" and        }
                          { solutionPTR (to hold data and solution).       }
  New(solutionPTR);

  StartFit := 1; EndFit := 1000;  { These parameters can be varied.        }
  SeriesSize := 6;                { SeriesSize rarely needs to be more than}
  IntegrationLevel := 2;          { ten, and for well behaved data can be  }
                                  { the equal to 2*(n+1) where n is the    }
                                  { number of exponentials being fitted -  }
                                  { defined here as IntegrationLevel.      }
  even_number_points;
  dummy_data(3,2,10);             { Creates a double exponential function  }
                                  { with constants -3 and -2, and arbitr-  }
                                  { ary noise factor 1.                    }
  Get_Legendre_Series;            { Get finite Legendre representation of  }
                                  { data.                                  }

  clrscr;
  for i := 0 to SeriesSize do
     writeln('LegendreSeries[',i,']:',LegendreSeries[i]);
  readln;

  Legendre_Fit;
  Integrate_Legendre_Series;

  clrscr;
  writeln('Integrated Legendre Series');
  for i := 0 to SeriesSize do begin
    for j := 0 to IntegrationLevel do begin
      GoToXY(1+j*12,4+i);
      write(IntegratedSeries[j,i]:6:2);
    end;
  end;
  readln;

  startgraph(true,0,0,0,0); { Substitute appropriate graphics              }
                        {initialisation, and graph plotting routines here. }
  For i := Startfit to EndFit do
       PutPixel(i,round(DataPTR^[i]/10)+screen_height div 2,15);
  For i := Startfit to EndFit do
       LineTo(i,round(SolutionPTR^[i]/10)+screen_height div 2);
  readln;
  restoreCRTmode;

  Find_Integral_Equation;

  writeln('This is the matrix used for finding the coefficients of the integral');
  writeln('equation: the coefficients are -1 and the first row of the matrix.');
  writeln;
  writeln('These are also necessarily the coefficients of the characteristic');
  writeln('equation of the linear homogeneous differential equation with constant');
  writeln('coefficients whose time dependent solution is the sum of exponentials.');
  for i := 0 to IntegrationLevel do begin
    for j := 0 to IntegrationLevel do begin
      GoToXY(1+j*24,8+i);
      write(CalcMatrix[j,i]:6:4);
    end;
  end;
  readln;

  find_roots;

  dispose(dataPTR);
  dispose(solutionPTR);
  solutionPTR := NIL;
  dataPTR     := NIL;
  readln;
end.