lab 10 everything

Harish Sridharan CE 311K Lab 10 Lab 10 Fortran Code ! CE311K Introduction to Computer Methods ! Lab10 - Least Square Method ! 3/25/08 Program leastsquare Implicit none ! Declaration of variables Real time(100), flux(100), x(100), y(100) Real r2,height,v,syx,a0,a1 Integer n, i ! Open I/O files Open (unit=5, file= 'input.txt' , status= 'old' ) Open (unit=6, file= 'out.dat' , status= 'unknown' ) Read (5, *) n print *,n Write (6,110) Write (6,100) ! Initialize the arrays Do i=1, n Read (5,*) time(i),flux(i) x(i) = log (time(i)) y(i) = log (flux(i)) Write (6,*) time(i),flux(i),x(i),y(i) Enddo close (5) Call LSquares(x, y, n,a0,a1,r2,syx) y(n+1)= exp (a0+a1*( log (5.))) y(n+2)= exp (a0+a1*( log (50.))) y(n+3)= exp (a0+a1*( log (500.))) write (6,*) 'flux at 5hrs' , y(n+1) write (6,*) 'flux @ 50hrs' , y(n+1) write (6,*) 'flux @ 500hrs' , y(n+3) 100 Format ( '====================================================' ) 110 Format (3x, 'time(h)' , 6x, 'flux(mg/h)' , 6x, 'ln(Time(h))' , 4x, 'ln(Flux(mg/h))' ) End Subroutine LSquares(x, y, n,a0,a1,r2,syx) Implicit none Real x(100), y(100) Real a0, a1, r2

Real sumx, sumy, sumxy, sumx2, xm, ym,syx Real sr, st Integer i, n sumx=0 sumy=0 sumxy=0 Do i=1,n sumx=sumx+x(i) sumy=sumy+y(i) sumxy=sumxy+x(i)*y(i) sumx2=sumx2+x(i)**2 End Do xm=sumx/n ym=sumy/n a1=((n*sumxy)-(sumx*sumy))/((n*sumx2)-(sumx**2)) a0=ym-(a1*xm) Do i=1,n Sr=Sr+(y(i)-a0-a1*x(i))**2 St=St+(y(i)-ym)**2 End do r2=(St-Sr)/St syx= SQRT (sr/(n-2)) write (6,*) 'a0=' , a0, 'a1=' , a1, 'r2=' , r2, 'sr=' , sr, 'st=' , st Return End Lab 10 Fortran Output

maxt = max(t); % c. With the max value of t, create a vector called xfit that starts at 0, % steps by 0.01, and ends at maxt. xfit = linspace(0,maxt,20); % d. Create a vector yfit that uses to the coefficients you found earlier with polyfit()to % determine the value of the polynomial at each of the xfit values. You % will be using the polyval() function to do this. Again, there is an example % program for you to compare with or you can use Matlab's built-in help % function. yfit = polyval(coef,xfit); %% 4. Plotting % a. Generate a figure using the figure() function. Simply place any number in % the paranthesis. figure(1) % b. Plot the input data as blue circles and the fitted data as a red line with a % linewidth of 2. If you forgot how to do this, try typing "help plot" on % the Command Window. plot(t,flux, 'o' ,xfit,yfit, 'r-' , 'Linewidth' ,2) % c. As always, give your plot axis labels. Remember the data has been % log-transformed so your axis labels should have some mention of this. You % will be using the xlabel() and ylabel() functions. xlabel( 'ln(Time)' ) ylabel( 'ln(Flux)' ) % d. Give your figure a legend so people can understand it. Each title is % given in single quotes and the first title you give corresponds to the % first set of data you plotted in the legend() function. legend( 'data' , 'polynomial fit' ) %% 5. Summary Data % a. Get the r-values matrix by using the corrcoef() function with you % transformed input data. Store this matrix in r. R12 = corrcoef([t flux]); % b. Now find r^2. Take the element in row 1, column 2 of r and square it. % Store this in a new varible called r2.

r12 = R12(1,2); r12sq = r12^2; %% 6. Output of Summary Data % a. Print out the equation of the line using the values from polyfit(). % Remember the first value given after using polyfit() corresponds to the % coefficient in front of the x with the greatest degree. For example, if % your coefficients matrix was [6.5 3.4] then your slope would be 6.5 and % your y-intercept 3.4. Use fprintf() to write out a statement that looks % like an equation for a line in slope-intercept form: y = mx + b. a0 = coef(1,2); fprintf( 'a0 is %7.4f\n' , a0) a1 = coef(1,1); fprintf( 'a0 is %7.4f\n' , a0) fprintf( 'y = %7.4f(t) + %7.4f\n' , a1, a0) % b. Finally, output the value of r2 again using fprintf(). fprintf( 'R^2 = %7.4f' , r12) %% Assignment % 1. Reading in from Excel file % a. On line 23, read in the flux values from Excel using xlsread() % 2. Transforming the Data % a. Take the natural log of your time and flux values on lines 31 and 32 % b. Determine the number of entries n on line 37 % 3. Fitting % a. Get the coefficients for your polynomial fit on line 47 % b. Use the max() function to get the maximum value of your transformed % time variable on line 52 % c/d. Get your x- and y-fit on lines 57 and 65 % 4. Plotting % a/b. Generate a figure and plot on it on lines 72 and 78 % c/d. Give your plot x and y labels as well as a legend on lines 84, 85, % and 91 % 5. Summary Data % a. Find the r-values matrix on line 98 % b. Locate the correct element in the r-values matrix and square it on % line 103 % 6. Output of Summary Data % a. Print out the equation of the linear fit on line 114 % b. Print out the correlation coefficient on line 118 %% To Turn-In % 1. Copy and paste this code into your master Word Document % 2. Copy and paste your output in the Command Window to your master Word % Document % 3. Copy and paste the figure into your master Word Document