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')
          READ (5,*) NPTS
          DO 10 I=1,NPTS
             READ (5,*) XCOOR(I),YCOOR(I)
     10   CONTINUE
          PRINT *,'ENTER DEGREE OF POLYNOMIAL '
          READ *,DEG
          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    '
          READ*,NUM,GU,TOL,MNIT

    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  '
          read (*,1000) filenm
          open (5,file=filenm)

    cccc  read in data file and convert counts to integer                 

          do 10 i=1,m
            read (5,*) t(i), cttmp(i)
            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'
    C     READ*,C,N,A,B,STEP,MNIT,TOL
          DATA C,N,A,B,STEP,MNIT,TOL/5,200,0.15,0.05,5,200,0.001/
          D=C/(N*A*B)

    CCCC  Print table heading

          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)
          REAL E,X1,X2,GRAD,TOL
          INTEGER J,MNIT
          DO 300 J=0,MNIT
             GRAD=(F(X1+0.01,E)-F(X1,E))/0.01
             X2=X1-(F(X1,E)/GRAD)
             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)
          REAL X1,X2,X3,E,GRAD,OLD,TOL
          INTEGER J,MNIT
          DO 400 J=0,MNIT
             GRAD=(F(X2,E)-F(X1,E))/(X2-X1)
             X3=X1-F(X1,E)/GRAD
             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)
          REAL X1,X2,X3,E,GRAD,OLD,TOL
          INTEGER J,MNIT
          DO 500 J=0,MNIT
             IF ((X2-X1).EQ.0.0) THEN
                J=0
                X3=0
                RETURN
             ENDIF
             GRAD=(F(X2,E)-F(X1,E))/(X2-X1)
             X3=X1-F(X1,E)/GRAD
             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.

Advertisements

Leave a Comment »

No comments yet.

RSS feed for comments on this post. TrackBack URI

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Create a free website or blog at WordPress.com.

%d bloggers like this: