c
clsq.f
c
c
module arrays
real, allocatable :: x(:), y(:)
integer np
c
c x - array containing independent variable for data pairs
c y - array containing the dependent variable for data pairs
c point
c np - number of data points
c
end module
c
program leastsq
use arrays
implicit none
c
c Program designed to perform a least squares fit of a cubic equation
c to some data.
c
c John Mahaffy 4/6/97
c
c Variables
c
c x - array containing independent variable for data pairs
c y - array containing the dependent variable for data pairs
c np - number of data points
c coef - Coefficients in the third order polynomial approximating the data
c y = coef(0) + coef(1)*x + coef(2)*x**2 + coef(3)*x**3
c
real coef(0:3)
integer i
c
call getdata
call cubicfit(x,y,np,coef(0))
write(*,*) 'For the equation:'
write(*,*) 'y = coef(0) + coef(1)*x + coef(2)*x**2',
& ' + coef(3)*x**3'
write(*,2000) (i,coef(i),i=0,3)
2000 format (/,' coef(',i1,') = ',1p,e15.7)
stop
end
subroutine getdata
use arrays
implicit none
c
c Get Experimental Data
c
c Variables
c
c x - time (sec) of measurement
c s - distance of fall measured (meters)
c np - number of data points
c fname - input file name
c
integer i, iend
character*32 :: fname='visc.tab'
logical lexist
c
c===============================================================
c Insert lines here to prompt for and read the input file name
c===============================================================
write(*,2000, advance='no')
2000 format (' Data File Name: ')
read *, fname
c
c Inquire optional for this homework
c
inquire (file=fname,exist=lexist)
if (.not.lexist) then
print *, trim(fname),' does not exist.'
stop
endif
open(11,file=fname)
c
c Count Data Pairs in the input file
c
np=-1
do while (iend.eq.0)
read(11,*,iostat=iend)
np = np + 1
enddo
if (iend.gt.0.or. np.lt.1) then
print *, 'Empty or Bad File, check fall.data'
stop
endif
c
c Rewind, Allocate Arrays, and Read data
c
rewind (11)
print *, np, ' data points read.'
allocate (x(1:np),y(1:np))
do i = 1,np
read (11,*) x(i),y(i)
enddo
c
end
c
subroutine cubicfit (xd,yd,ndata,c)
implicit none
c
c Take corresponding data points from the arrays xd and yd and fit
c them with the following equation:
c
c y = c(1) + c(2) * x + c(3) * x**2 + c(4)*x**3
c
real xd(ndata),yd(ndata),aa(4,4),c(4),fs,suma,sumb,sqrt,rsid
integer ipvt(4),ir,is,ij,ndata, info
c
c Input Arguments:
c xd - x values for data points
c yd - y values for data point pairs
c ndata - number of data points
c
c Output Arguements:
c c - array containing the three coefficients in the 2nd order
c polynomial that provide the "best" fit to the data from a
c least squares method. Note that it also temporarily holds
c the values of the right hand sides of each Least Squares
c equation.
c Other Key Variables:
c aa - matrix containing coefficients of the system of equations
c generated by the least squares method
c suma - a variable that tallies the sum that is needed to generate
c each element of aa.
c sumb - a variable that tallies the sum that is needed for the right
c hand side of each Least Squares equation.
c ir - an index that is used to keep track of the equation number
c within the system of equations .
c is - an index used to track the coefficient number within a given
c Least Squares equation.
c
c
c DO loop 55 generates terms for each of the 4 Least Squares Equations
c
do 55 ir=1,4
c
c DO loop 45 generates the right hand side for a given equation
c
sumb=0.
do 45 ij=1,ndata
45 sumb=sumb+yd(ij)*xd(ij)**(ir-1)
c(ir)=sumb
c
c DO loop 50 generates the coefficients a given equation
c
do 50 is=1,4
suma=0.
do 52 ij=1,ndata
52 suma=suma+xd(ij)**(is-1)*xd(ij)**(ir-1)
50 aa(ir,is)=suma
55 continue
c
c Solve the Least Squares Equations
c
call sgefa(aa,4,4,ipvt,info)
call sgesl(aa,4,4,ipvt,c ,0)
c
c Calculate a measure of mean error between the data and the curve
c
rsid=0.
do 65 ir=1,ndata
fs=0.
do 60 is=1,4
60 fs=fs+c(is)*xd(ir)**(is-1)
65 rsid=rsid+(yd(ir)-fs)**2
c
rsid=sqrt(rsid/float(ndata-1))
c
write(6,2001) rsid
2001 format(' Fit to 3rd order polynomial has a mean error of',
$ 1p,e12.5)
return
end
c
c
c