c<html>
c<head><title>lsq3.f</title></head>
c<body>
c<pre>
program leastsq
implicit none
c
c Program designed to perform a least squares fit of a quadratic equation
c to some data. In this specific example the data is the result of
c measuring the location of a falling object at various times. The result
c of interest is gravitational acceleration obtained from the coefficient
c off t**2.
c
real t(50),z(50),coef(0:2),g
integer np,npoly
c
c Variables
c
c t - array containing times (seconds) at which positions are measured
c z - measured distance of the falling object from the nominal point of
c release (fit may let you estimate any offset in the actual release
c point
c np - number of data points
c coef - Coefficients in the second order polynomial approximating the data
c g - approximation to the gravitational acceleration deduced from data
c
c
npoly=2
call getdata(t,z,np)
call fit(t,z,np,npoly,coef(0))
g=2*coef(2)
write(6,2000) g
2000 format (' Predicted value of g = ', f6.3,'m/s**2')
write(6,2001) coef(0),coef(1)
2001 format(' Predicted initial offset = ',f7.4,'m/s'/
# ' Predicted initial velocity = ',f7.4,'m')
stop
end
subroutine getdata(t,s,n )
implicit none
c
c Get Experimental Data
c
c Input Arguments: NONE
c Output Arguments:
c t - time (sec) of measurement
c s - distance of fall measured (meters)
c n - number of data points
c
integer n
real t(*),s(*)
open(11,file='fall.data')
n=0
10 read (11,*,end=20 ) t(n+1),s(n+1)
n=n+1
go to 10
20 return
end
subroutine fit (xd,yd,ndata,npoly,c)
implicit none
c
c This version of the subroutine has been modified to improve speed
c on vector or parallel machines, at the cost of extra memory (see
c the added array xdp)
c
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(npoly+1)*x**npoly
c
real xd(ndata),yd(ndata),aa(npoly+1,npoly+1),c(*),fs,
& sqrt,rsid,xdsum(0:(npoly+npoly)),xdp(ndata,npoly+npoly)
integer ipvt(npoly+1),ir,is,ij,ndata, info, npoly,i,j,npow, neq
c
c NOTE: The use of aa(npoly+1,npoly+1), ipvt(npoly+1), etc. above is
c an example of the FORTRAN 90 AUTOMATIC ARRAY feature. It is
c a short-hand way of using the ALLOCATE statement. Space for
c aa, c, and ipvt is assigned when SUBROUTINE FIT is entered and
c is DEALLOCATED automatically when a RETURN is executed.
c Note the flexibility available in assigning the array sizes.
c A compiler without FORTRAN 90 will not permit using npoly to
c specify the sizes of these arrays unless they (aa, xdsum, ipvt)
c also appear in the argument list. Beware of use of this
c convenient FORTRAN 90 feature in frequently used subroutines.
c A fair amount of machine time is associated with each ALLOCATE
c and DEALLOCATE.
c
c Input Arguments:
c xd - x values for data points
c yd - y valuess for data point pairs
c ndata - number of data points
c npoly - order of polynomial fit (highest power of x in polynomial)
c
c Output Arguements:
c c - array containing the coefficients in the chosen 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 xdsum - array containing the sums of all needed powers of the elements
c in xd
c aa - matrix containing coefficients of the system of equations
c generated by the least squares method
c ipvt - an array used by Linpack routines to keep track of "pivoting"
c operations during Gauss Elimination.
c neq - number of linear equations to be solved
c
c Zero Sum counters
c
npow=npoly+npoly
neq=npoly+1
xdsum(0)=ndata
do 20 j=1,npow
20 xdsum(j)=0.
do 21 j=1,neq
21 c(j)=0.
c
c Generate the needed powers of xd(i)
c
do 30 i=1,ndata
xdp(i,1)=xd(i)
do 30 j=2,npow
30 xdp(i,j)=xdp(i,j-1)*xd(i)
c
c Sum the Powers of xd
c
do 40 j=1,npow
c<a name=1><font color=FF0000>
40 xdsum(j)=sum(xdp(1:ndata,j))
c</font>
c Generate the right hand sides for all equations
c
c(1)=sum(yd(1:ndata))
do 45 j=2,neq
c<a name=2><font color=FF0000>
45 c(j)=dot_product(yd(1:ndata),xdp(1:ndata,j-1))
c</font>
c DO loop 50 generates the coefficients a given equation
c
do 50 i=1,neq
do 50 j=1,neq
50 aa(i,j)=xdsum(i+j-2)
c
c Solve the Least Squares Equations
c
call sgefa(aa,neq,neq,ipvt,info)
call sgesl(aa,neq,neq,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,neq
60 fs=fs+c(is)*xd(ir)**(is-1)
65 rsid=rsid+(yd(ir)-fs)**2
rsid=sqrt(rsid/float(ndata-1))
write(6,2001) rsid
2001 format(' Fit to 2nd order polynomial has a mean error of',
$ 1p,e12.5)
return
end
c</pre>
c</body>
c</html>
Copyright ©2017 - STI - Todos os direitos reservados