# Rcjp's Weblog

## February 13, 1991

### BSc Fortran Numerics Course – Fitting

Filed under: physics — rcjp @ 11:18 am

I remember well this next program, I used it to fit a high order polynomial to some simple data – the result was a very wiggly line that ran through all the points! Though, I *think* this program is wrong – the original was in Propero and this looks like a half finished port of that program.

``````
CCCC  ******************************************************************
CCCC  **      DRAWS A GRAPH OF DATA HELD IN FILE 'PDATA.DAT'          **
CCCC  **         AND DRAW LEAST SQUARES FIT TO POLYNOMIAL             **
CCCC  ******************************************************************

CCCC  Data file includes the number of points
CCCC  and X,Y coordinates

INTEGER DEG
DIMENSION XCOOR(20),YCOOR(20),COEFF(0:20,0:20),CONST(0:20)

OPEN (5,FILE='PDATA.DAT',STATUS='OLD')
DO 10 I=1,NPTS
10   CONTINUE
PRINT *,'ENTER DEGREE OF POLYNOMIAL '
CALL FITQ (XCOOR,YCOOR,NPTS,DEG,COEFF,CONST)
CALL GAUSPC(A,B,DEG+1,1E-3,G)
CALL SCRSET
CALL GRAFIT (X,Y,NPTS,DEG,MX,MY,MNX,MNY,A,B,G)
CALL GET_KEY@(K)
CALL TEXT_MODE@
CLOSE (5)
END

SUBROUTINE FITEQNS(XCOOR,YCOOR,NPTS,DEG,
*     MATRIX,YPXSUM)

CCCC  ==================================================================
CCCC  Creates a set of simultaneous equations which when solved give the
CCCC  best fit of a polynomial y = a + a(1) x + a(2) x**2 + ...
CCCC  of degree DEG to the coordinates (XCOOR,YCOOR).
CCCC  Returned is a matrix of the coefficients of the equations -
CCCC  MATRIX(I,J) and the associated constants - YPXSUM(J).
CCCC  Where I and J take values 0 -> DEG.
CCCC  ==================================================================
CCCC  DESCRIPTION OF VARIABLES
CCCC  ------------------------
CCCC  INTEGER   DEG             Degree of polynomial fit
CCCC  NPTS        Number of data points
CCCC  REAL      XCOOR()         X - coordinates
CCCC  YCOOR()         Y - coordinates
CCCC  MATRIX(,)       Matrix of simultaneous eqn coeff
CCCC  YPXSUM()        Sum of Y times X to power n (n=0->DEG)
CCCC  XPWSUM()        Sum of X to power n (n=0->DEG)
CCCC  ==================================================================
INTEGER           DEG
REAL              MATRIX
DIMENSION         XCOOR(*),YCOOR(*),
*     MATRIX(0:100,0:*),YPXSUM(0:*),XPWSUM(0:*)
C
C     Find the sum of the x-coordinate to the power of n
C     where n varies between 0 and 2 times the degree of polynomial
C     for all coordinates.
C
XPWSUM(0)=NPTS
DO 20 I=1,2*DEG
DO 10 J=1,NPTS
XPWSUM(I)=XPWSUM(I)+XCOOR(J)**I
10      CONTINUE
20   CONTINUE
C
C     Find the sum of (x-coordinate to the power n) times by the
C     y-coordinate for all coordinates.
C
DO 110 I=0,DEG
DO 100 J=1,NPTS
YPXSUM(I)=YPXSUM(I)+YCOOR(J)*XCOOR(J)**I
100     CONTINUE
110  CONTINUE
C
C     Create matrix of coefficients of simultaneous equations
C     of the form a(0,0) + a(1,0) + ... + a(i,0) = const(0)
C     a(0,1) + a(1,1) + ... + a(i,1) = const(1)
C     :        :              :       :
C     a(0,j) + a(1,j) + ... + a(i,j) = const(j)
C
C     Note the const() is YPXSUM() and is already evaluated.
C
DO 210 J=0,DEG
DO 200 I=0,DEG
MATRIX(I,J)=XPWSUM(I+J)
200     CONTINUE
210  CONTINUE
END

CCCC  Solve the best fit equations

SUBROUTINE GAUSPC (A,B,N,TOL,G)

INTEGER I,J,K,N,L
REAL A,B,X,M,F,TOL,C,D,G
DIMENSION A(20,20),G(20),B(20)
CCCC  PARTIAL PIVOTAL CONDENSATION

DO 320 J=1,N
L=J
DO 300 I=1,N
IF (A(I,J).GT.A(L,J)) L=I
300     CONTINUE
DO 310 K=1,N
C=A(J,K)
A(J,K)=A(L,K)
A(L,K)=C
310     CONTINUE
D=B(J)
B(J)=B(L)
B(L)=D
320  CONTINUE

CCCC  Gaussian elimination method

DO 450 K=1,N-1
DO 400 I=K+1,N
IF (ABS(A(K,K)).LT.TOL) STOP
M=A(I,K)/A(K,K)
DO 350 J=K,N
A(I,J)=A(I,J)-M*A(K,J)
350        CONTINUE
B(I)=B(I)-M*B(K)
400     CONTINUE
450  CONTINUE

CCCC  Back substitution

G(N)=B(N)/A(N,N)
DO 110 J=N-1,1,-1
f=0
DO 100 L=J+1,N
F=F+A(J,L)*G(L)
100     CONTINUE
G(J)=1/A(J,J)*(B(J)-F)
110  CONTINUE
END

CCCC  Initialise screen and clear it.

SUBROUTINE SCRSET
CALL VGA@
CALL CLEAR_SCREEN@
END

CCCC  Plots the data, draws the axis
CCCC  and the line of best fit

SUBROUTINE GRAFIT (X,Y,NPTS,DEG,MX,MY,MNX,MNY,A,B,G)
INTEGER I,NPTS,XP,YP,Y1,Y2,OLDX,OLDY,DEG
REAL X,Y,SCALEX,SCALEY,MX,MY,MNX,MNY,A,B,G,C,D
DIMENSION X(20),Y(20),A(20,20),B(20),G(40)

IF (MNX.GT.0.0) MNX=0.0
IF (MNY.GT.0.0) MNY=0.0
SCALEX=600.0/(MX-MNX)
SCALEY=310.0/(MY-MNY)

CCCC  Draw axis
XP=INT(20-MNX*SCALEX)
YP=INT(330+MNY*SCALEY)
CALL DRAW_LINE@(XP,20,XP,330,15)
CALL DRAW_LINE@(20,YP,620,YP,15)
CCCC  Mark points
DO 500,I=1,NPTS
XP=INT(X(I)*SCALEX+20-MNX*SCALEX)
YP=INT(330+MNY*SCALEY-Y(I)*SCALEY)
CALL DRAW_LINE@(XP-2,YP,XP+2,YP,15)
CALL DRAW_LINE@(XP,YP-2,XP,YP+2,15)
500  CONTINUE

CCCC  Draw Curve

DO 620 C=MNX,MX+0.005,0.01
D=0
DO 610 I=1,DEG+1
IF (C.EQ.0.0) THEN
D=G(1)
GOTO 610
ENDIF
D=D+G(I)*C**(I-1)
610     CONTINUE
XP=INT(C*SCALEX+20-MNX*SCALEX)
YP=INT(330+MNY*SCALEY-D*SCALEY)
IF (OLDX.EQ.0.0) GOTO 615
CALL DRAW_LINE@(OLDX,OLDY,XP,YP,15)
615     OLDX=XP
OLDY=YP
620  CONTINUE
PRINT *,'DEGREE OF POLYNOMIAL FIT ',DEG
END
``````

## Iterative Root Finding

``````
CCCC  ******************************************************************
CCCC  **               EVALUATION OF SQUARE ROOT                      **
CCCC  **                  BY ITERATIVE METHOD                         **
CCCC  **                                                              **
CCCC  **   Improved guess= 0.5 (guess+number/guess)                   **
CCCC  **                                                              **
CCCC  ******************************************************************

CCCC  Initialise variables

REAL NUM,GU,TOL,VAL
INTEGER FLAG,MNIT,N

CCCC  Enter data

PRINT*,'ENTER NUMBER, GUESS, TOLERANCE, MAX.ITERATIONS    '

CCCC  Call square root routine

CALL SQ(NUM,TOL,MNIT,GU,FLAG,VAL,N)
IF (FLAG.GE.2) THEN
IF (FLAG.EQ.2) PRINT *,'NUMBER TOO SMALL OR NEGATIVE'
IF (FLAG.EQ.3) PRINT *,'GUESS TOO SMALL OR NEGATIVE'
STOP
ELSE
IF (FLAG.EQ.1) PRINT *,'WARNING:MAXIMUM NO. OF ITERATIONS',
1  ' WAS REACHED'
ENDIF
PRINT*,'VALUE IS ',VAL,'  NO. OF ITERATIONS',N
END

CCCC  Square root routine

SUBROUTINE SQ(NUM,TOL,MNIT,GU,FLAG,VAL,N)
REAL NUM,TOL,GU,VAL
INTEGER MNIT,FLAG,N
IF (NUM.LT.TOL) THEN
FLAG=2
RETURN
ELSE
IF (GU.LT.TOL) THEN
FLAG=3
RETURN
ENDIF
ENDIF
DO 10 N=0,MNIT
VAL=(GU+NUM/GU)*0.5
IF (ABS(GU-VAL).LT.TOL) RETURN
GU=VAL
10  CONTINUE
FLAG=1
RETURN
END
``````

But I think writing our own routines was just an exercise; we used the NAG library too

``````
cccc*******************************************************************
cccc           Calculates the best fit of data in requested file
cccc           to equation given in lsfun1, uses NAG E04FDF
cccc*******************************************************************
cccc===================================================================
c     Equation is the number of counts from parent-daughter decay
c     Find best fit of data using constants lambda a & b the decay
c     rates and c0 the initial decay rate.
c     filenm - filename holding experimental data
c     t()    - time
c     ct()   - counts at time t
c     m      - number of data points
c     n      - number of unknown constants in equation
cccc===================================================================

integer t(50), ct(50)
parameter (m=41, n=3, liw=1, lw=7*n+n*n+2*m*n+3*m+n*(n-1)/2)
double precision cttmp(50), iw(liw), w(lw), x(n)
character filenm*8
external e04fdf
common t, ct

cccc  get filename for experimental data

print *,'Enter filename containing experimental data  '
open (5,file=filenm)

cccc  read in data file and convert counts to integer

do 10 i=1,m
ct(i)=int(cttmp(i))
10 continue

cccc  initial guesses

x(1)=1.0d-3
x(2)=3.0d-3
x(3)=10040.0d0
ifail=1

call e04fdf(m, n, x, fsumsq, iw, liw, w, lw, ifail)

print *,'Minimum point is ',(x(j),j=1,n)
c      print 1001,'Minimum point is ',x(1),x(2),x(3)
close (5)
1000 format (a)
1001 format (1x,'lambda a ',d9.3,'  lambda b ',d9.3,
*         '  initial count ', d9.3)
end

subroutine lsfun1(m, n, xc, fvecc)
cccc*******************************************************************
cccc  Calculates the difference between the predicted count rate
cccc  using the estimates constants lambda a,b and the initial count
cccc*******************************************************************
cccc===================================================================
c     la    - lambda a
c     lb    - lambda b
c     c0    - initial count
c     fvecc - difference between calculated and experimental values
cccc===================================================================
integer t, ct
double precision fvecc(50), xc(3), la, lb, c0
dimension t(50),ct(50)
common t, ct

la=xc(1)
lb=xc(2)
c0=xc(3)
do 10 i=1,m
fvecc(i)=((lb*2-la)*exp(-la*t(i))-
*               lb*exp(-lb*t(i)))*c0/(lb-la)-ct(i)
10 continue
end
``````

## Interpolation, Extrapolation and Bisection

``````
CCCC  ******************************************************************
CCCC  **  PROGRAM TO DEMONSTRATE A VARIETY OF METHODS FOR SOLVING     **
CCCC  **             TRANSCENDENTAL EQUATIONS                         **
CCCC  ******************************************************************
CCCC
CCCC  THIS PROGRAM GIVES THE ANGLE OF A GALVONOMETER (THI) FOR A GIVEN
CCCC  CURRENT THROUGH IT BY SOLVING THE TRANSCENDENTAL EQUATION
CCCC  THI=(NAB/C)ICOS(THI)
CCCC  THIS PROGRAM HOLDS DATA FOR THE VARIABLES BUT MAY BE USED WITH
CCCC  OTHER VALUES BY REMOVING THE COMMENT MARKS IN THE MARGIN BELOW
CCCC  AND THE DATA STATEMENT.

CCCC  Initialise variables

REAL C,A,B,TOL,D,THI,X1,X2,I,PIB2,E,THI1,THI2,THI3,THI4
INTEGER N,STEP,MNIT,J,J1,J2,J3,J4
PARAMETER (PIB2=3.141592/2)

CCCC  Enter variable values

C     PRINT *,'ENTER C,N,A,B,LENGTH OF TABLE,MAX.NO ITERATIONS AND '
C     1  ,'TOLERANCE'
DATA C,N,A,B,STEP,MNIT,TOL/5,200,0.15,0.05,5,200,0.001/
D=C/(N*A*B)

PRINT99
PRINT100
PRINT101
PRINT 99

CCCC  Submit current values to subroutines
CCCC  and display solutions

DO 10 I=0,4*D,4*D/STEP
E=I/D
X1=0
X2=PIB2
IF (F(X1,E)*F(X2,E).GT.0.0) THEN
PRINT *,I,'NO ROOTS'
GOTO 10
ENDIF
CALL BISECT (E,MNIT,X1,X2,J,THI)
J1=J
THI1=THI
X1=0
X2=PIB2
CALL NEWT (E,MNIT,X1,J,THI,TOL)
J2=J
THI2=THI
X1=0
X2=PIB2
CALL INTERP (E,MNIT,X1,X2,J,THI,TOL)
J3=J
THI3=THI
X1=0
X2=PIB2
CALL EXTRAP(E,MNIT,X1,X2,J,THI,TOL)
J4=J
THI4=THI
PRINT102,I,J1,THI1,J2,THI2,J3,THI3,J4,THI4
10   CONTINUE
PRINT99
99   FORMAT (1X,62('-'))
100  FORMAT (1X,'CURRENT',3X,'BISECTION',5X,'NEWTON',3X,
1     'INTERPOLATION',3X,'EXTRAPOLATION')
101  FORMAT (11X,'N',3X,'THI',6X,'N',2X,'THI',7X,'N',4X,'THI',
1     8X,'N',3X,'THI')
102  FORMAT (1X,F7.4,1X,I3,1X,F7.4,2X,I3,1X,F7.4,2X,I3,1X,F7.4,
1     5X,I3,1X,F7.4)
END

CCCC  Solution by bisection method

SUBROUTINE BISECT(E,MNIT,X1,X2,J,X3)
REAL D,X1,X2,X3,E,OLD,TOL
INTEGER MNIT,J
DO 200 J=0,MNIT
X3=0.5*(X1+X2)
IF (F(X1,E)*F(X3,E).LT.0.0) THEN
X2=X3
ELSE
X1=X3
ENDIF
IF (J.NE.0) THEN
IF (ABS(F(OLD,E)-F(X3,E)).LE.TOL) RETURN
ENDIF
OLD=X3
200  CONTINUE
END

CCCC  Solution by Newton's method

SUBROUTINE NEWT(E,MNIT,X1,J,X2,TOL)
INTEGER J,MNIT
DO 300 J=0,MNIT
IF (J.NE.0) THEN
IF (ABS(F(X2,E)-F(X1,E)).LE.TOL) RETURN
ENDIF
X1=X2
300  CONTINUE
END

CCCC  Solution by interpolation

SUBROUTINE INTERP(E,MNIT,X1,X2,J,X3,TOL)
INTEGER J,MNIT
DO 400 J=0,MNIT
IF (F(X1,E)*F(X3,E).LT.0.0) THEN
X2=X3
ELSE
X1=X3
ENDIF
IF (J.NE.0) THEN
IF (ABS(F(OLD,E)-F(X3,E)).LE.TOL) RETURN
ENDIF
OLD=X3
400  CONTINUE
END

CCCC  Solution by extrapolation

SUBROUTINE EXTRAP(E,MNIT,X1,X2,J,X3,TOL)
INTEGER J,MNIT
DO 500 J=0,MNIT
IF ((X2-X1).EQ.0.0) THEN
J=0
X3=0
RETURN
ENDIF
IF (ABS(X3-X1).LT.ABS(X3-X2)) THEN
X2=X3
ELSE
X1=X3
ENDIF
IF (J.NE.0) THEN
IF (ABS(F(OLD,E)-F(X3,E)).LE.TOL) RETURN
ENDIF
OLD=X3
500  CONTINUE
END

CCCC  Form of transcendental equation

FUNCTION F(X,E)
REAL F,X,E
F=X-E*COS(X)
END
``````

There were perhaps another dozen or so other programs for e.g. calculating a certain branch of arctan by a formula etc. but not very interesting by today’s coding standards.