Chapter 19 The Normal Equations Can Be Derived as

CHAPTER 19 19.1 The normal equations can be derived as  11  2.416183 2.018098 2.416183 2.018098  A0   83.9  6.004565 0.017037  A1  = 15.43934  0.017037 4.995435  B1  10.81054    which can be solved for A0 = A1 = B1 = 7.957538 -0.6278 -1.04853 The mean is 7.958 and the amplitude and the phase shift can be computed as C1 = ( −0.6278) 2 + ( − 104853 . ) 2 = 1.222 12 hrs  − 1.04853  θ = tan −1  = 8.06 hr  + π = 2.11 radians ×  − 0.6278  π Thus, the final model is  2π  f ( t ) = 7.958 + 1.222 cos (t + 8.06)  24  The data and the fit are displayed below: 10 8 6 0 12 Note that the peak occurs at 24 − 8.06 = 15.96 hrs. 19.2 The normal equations can be derived as 24 1890 127 .279 5  0 −1  0 − 568.187   A0   350,265  − 1   A1  = − 381.864  5   B1   156.281  which can be solved for A0 = 195.2491 A1 = -73.0433 B1 = 16.64745 The mean is 195.25 and the amplitude and the phase shift can be computed as C1 = ( −73.0433) 2 + (16.6475) 2 = 74.916 180 d  16.6475  θ = tan −1  = 192 .8 d  + π = 3.366 radians ×  − −73.0433  π Thus, the final model is  2π  f ( t ) = 195.25 + 74.916 cos (t + 192.8)  360  The data and the fit are displayed below: 300 200 100 0 0 90 180 270 360 19.3 In the following equations, ω0 = 2π/T ∫0 sin( ω 0 t ) dt = − ω 0 [ cos( ω 0 t ) ] 0 T T T T ∫0 cos( ω 0 t ) dt = ω 0 [ sin( ω 0 t ) ] 0 T T T T = = −ω0 ω0 T T ( cos 2π − cos 0) = 0 ( sin 2π − sin 0) = 0  t sin( 2ω 0 t )  T sin 4π T  −  2 − − 0+0 sin ( ω 0 t ) dt  2 4 ω  0 2 4ω 0 1  0 0 = = = T T T 2 T ∫  t sin ( 2ω 0 t )  T sin 4π T  +  2 + −0−0 cos ( ω 0 t ) dt  2 4ω 0  2 4ω 0 1  0 0 = = = T T T 2 T ∫ ∫ T 0 cos( ω 0 t ) sin( ω 0 t ) dt T T 2  sin 2 ( ω 0 t )  sin 2π = −0=0  = 2ω 0 T  2Tω 0  0 19.4 a0 = 0 ak = 2 T ∫ T/2 −T / 2 − t cos( kω 0 t ) dt T /2  2 1 t  =−  cos k ω t + sin k ω t ( ) ( ) 0 0 T  ( kω ) 2 kω 0  0   T /2 bk = 2 T T/2 ∫−T /2 − t sin ( k ω 0 t ) dt T /2  2 1 t =−  sin k ω t − cos k ω t ( 0 ) k ω ( 0 )  T  ( kω 0 ) 2 0   T /2 On the basis of these, all a’s = 0. For k = odd, bk = 2 kπ For k = even, bk = − 2 kπ Therefore, the series is f (t) = − 2 1 2 1 sin( ω 0 t ) + sin ( 2ω 0 t ) − sin ( 3ω 0 t ) + sin ( 4ω 0 t ) + ⋅ ⋅ ⋅ π π 3π 2π The first 4 terms are plotted below along with the summation: 1 0 -2 0 2 -1 19.5 a0 = 0.5 ak = 2 2  0  1 ∫−1 − t cos( kπt ) dt + ∫0 t cos( kπt ) dt  0 1    cos( kπt ) t sin ( kπt )     cos( kπt ) t sin( kπt )    = 1 − − + + 2 2 kπ kπ       ( kπ )  −1  ( kπ ) 0   2 = ( cos kπ − 1) ( kπ ) 2 bk = 0 Substituting these coefficients into Eq. (19.17) gives f (t) = 1 12 12 12 − 2 cos( πt ) − cos( 3πt ) − cos( 5πt ) + ⋅ ⋅⋅ 2 2 π 9π 25π 2 This function for the first 4 terms is displayed below: 1 0.5 0 -2 -1 -0.5 0 1 2 19.6 0.7 0 2 4 2 4 2 4 2 4 -0.7 π 19.7 0.7 0 -0.7 -1.4 π 19.8 1.2 0.6 0 0 5 10 -0.6 19.9 0.4 0 -0.4 4 π 4 19.10 Here’s a Fortran 90 code that implements the DFT. It is set up to solve Prob. 19.11. PROGRAM DFourier IMPLICIT NONE INTEGER i,N REAL f(0:127),re(0:127),im(0:127),omega,pi,t,Tp,dt pi=4.*atan(1.) N=32 omega=2.*pi/N t=0. Tp=2.*pi dt=4.*Tp/N DO i=0,N-1 f(i)=sin(t) if (f(i).LT.0.) f(i)=0. t=t+dt END DO CALL DFT(f,N,re,im,omega) OPEN (UNIT=1,FILE='Prob1911.dat',STATUS='unknown') DO i=0,N-1 WRITE(1,*) i,f(i),re(i),im(i) END DO CLOSE(1) END SUBROUTINE DFT(f,N,re,im,omega) IMPLICIT NONE INTEGER k,nn,N REAL f(0:127),re(0:127),im(0:127),angle,omega DO k=0,N-1 DO nn=0,N-1 angle=k*omega*nn re(k)=re(k)+f(nn)*cos(angle)/N im(k)=im(k)-f(nn)*sin(angle)/N END DO END DO END 19.11 The results for the n = 32 case are displayed below: index 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 f(t) 0 0.7071 1 0.7071 0 0 0 0 0 0.7071 1 0.7071 0 0 0 0 0 0.7071 1 0.7071 0 real imaginary 0.3018 0 0 0 0 0 0 0 0 -0.25 0 0 0 0 0 0 -0.125 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.0518 0 0 0 0 0 0 0 0 0 21 22 23 24 25 26 27 28 29 30 31 0 0 0 0 0.7071 1 0.7071 0 0 0 0 0 0 0 -0.125 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.25 0 0 0 The runs for N = 32, 64 and 128 were performed with the following results obtained. (Note that even though we used a slow PC, we had to call the function numerous times to obtain measurable times. These times were then divided by the number of function calls to determine the time per call shown below) N 32 64 128 time (s) 0.09 0.37 1.48 A power (log-log) model was fit (see plot below) to this data to yield log(time) = −4.08 + 2.02 log(N). Thus, the result verifies that the execution time ∝ N2. 10 time 1 10 0.1 100 1000 N 0.01 19.12 Here’s a Fortran 90 code that implements the FFT. It is set up to solve Prob. 19.13. PROGRAM FFourier IMPLICIT NONE INTEGER i,N REAL f(0:127),re(0:127),im(0:127),omega,pi,t,Tp,dt pi=4.*ATAN(1.) N=32 t=0. Tp=2.*pi dt=4.*Tp/N DO i=0,N-1 re(i)=sin(t) if (re(i).LT.0.) re(i)=0. f(i)=re(i) t=t+dt END DO CALL FFT(N,re,im) DO i=0,N-1 PRINT *, i,f(i),re(i),im(i) END DO CLOSE(1) END SUBROUTINE FFT (N, x, y) IMPLICIT NONE INTEGER :: i,j,N,m,N2,N1,k,l REAL :: f(0:127),re(0:127),im(0:127),omega,pi,t,Tp,dt,xN,angle REAL :: arg,c,s,xt,x(0:n),y(0:n),yt xN=N m = INT(LOG(xN) / LOG(2.)) pi = 4. * ATAN(1.) N2 = N DO k = 1, m N1 = N2 N2 = N2 / 2 angle = 0. arg = 2 * pi / N1 DO j = 0, N2 - 1 c = COS(angle) s = -SIN(angle) DO i = j, N - 1, N1 l = i + N2 xt = x(i) - x(l) x(i) = x(i) + x(l) yt = y(i) - y(l) y(i) = y(i) + y(l) x(l) = xt * c - yt * s y(l) = yt * c + xt * s END DO angle = (j + 1) * arg END DO END DO j = 0 DO i = 0, N - 2 IF (i.LT.j) THEN xt = x(j) x(j) = x(i) x(i) = xt yt = y(j) y(j) = y(i) y(i) = yt END IF k = N / 2 DO IF (k.GE.j+1) EXIT j = j - k k = k / 2 END DO j = j + k END DO DO i = 0, N - 1 x(i) = x(i) / N y(i) = y(i) / N END DO END 19.13 Note that the results for the n = 32 case should be the same as for the DFT as in the first part of the solution of Prob. 19.11 as shown above. The runs for N = 32, 64 and 128 were performed with the following results obtained. (Note that even though we used a slow PC, we had to call the function numerous times to obtain measurable times. These times were then divided by the number of function calls to determine the time per call shown below) N time (s) 32 64 128 0.0135 0.031 0.068 A plot of time versus N log2N yielded a straight line (see plot below). Thus, the result verifies that the execution time ∝ N log2N. time 0.05 0 0 500 1000 N log2N 19.14 Using a similar approach to that described in Example 19.3, the Excel Chart Wizard and the Trendline tool can be used to create the following fit: 8 y = 10.051x-0.5304 R2 = 0.9515 4 0 0 10 20 19.15 Using a similar approach to Example 19.4, the following spreadsheet can be set up: T 0 8 16 24 32 40 T^2 0 64 256 576 1024 1600 T^3 0 512 4096 13824 32768 64000 T^4 0 4096 65536 331776 1048576 2560000 o 14.621 11.843 9.87 8.418 7.305 6.413 The Data Analysis Toolpack can then be used to generate SUMMARY OUTPUT Regression Statistics Multiple R 0.99999994 R Square Adjusted R Square Standard Error Observations 0.99999988 0.99999939 0.00239377 6 ANOVA df 4 1 5 SS 47.0093523 5.7302E-06 47.009358 MS 11.75234 5.73E-06 Coefficients 14.6208492 -0.4113267 0.0090115 -0.0001302 8.4432E-07 Std Error 0.00238902 0.0011012 0.00013149 5.1867E-06 6.4426E-08 t Stat 6120.018 -373.527 68.53234 -25.1078 13.10526 Regression Residual Total Intercept X Variable 1 X Variable 2 X Variable 3 X Variable 4 F Significance F 2050962 0.0005237 P-value Lower 95% Upper 95% 0.000104 14.59049395 14.6512 0.001704 -0.425318685 -0.39733 0.009289 0.007340736 0.010682 0.025342 -0.000196129 -6.4E-05 0.048483 2.57132E-08 1.66E-06 The polynomial along with the data can be plotted as 15 10 5 0 10 20 30 40 19.16 Linear regression can be implemented with the Data Analysis Toolpack in a fashion similar to Example 19.4. After setting the data ranges, the confidence interval box should be checked an set to 90% in order to generate 90% confidence intervals for the coefficients. The result is SUMMARY OUTPUT Regression Statistics Multiple R 0.98465 R Square 0.969535 Adjusted R 0.963442 Square Standard 1.625489 Error Observations 7 ANOVA Regression df 1 SS 420.4375 MS 420.4375 F 159.1232 Significance F 5.56E-05 Residual Total Intercept X Variable 1 5 6 13.21107 433.6486 2.642214 Coefficients Standard Error 1.373789 0.153594 t Stat 0.714286 1.9375 0.519938 12.6144 P-value Lower 95% Upper 95% 0.625298 5.56E-05 -2.81715 1.542674 4.245717 2.332326 Lower 90.0% -2.05397 1.628 Upper 90.0% 3.482538 2.247 The 90% confidence interval for the intercept is from -2.05 to 3.48, which encompasses zero. The regression can be performed again, but with the “Constant is Zero” box checked on the regression dialogue box. The result is SUMMARY OUTPUT Regression Statistics Multiple R 0.983813 R Square 0.967888 Adjusted R 0.801221 Square Standard 1.523448 Error Observations 7 ANOVA Regression Residual Total Intercept X Variable 1 df 1 6 7 SS 419.7232 13.92536 433.6486 MS 419.7232 2.320893 Coefficients Standard Error #N/A 0.064377 t Stat 0 2.008929 #N/A 31.20549 F 180.8456 Significance F 4.07E-05 P-value Lower 95% Upper 95% #N/A 7.19E-08 #N/A 1.851403 The data along with both fits is shown below: 30 20 10 0 0 19.17 Using MATLAB: >> x=[0 2 4 7 10 12]; >> y=[20 20 12 7 6 5.6]; 8 16 #N/A 2.166455 Lower 90.0% #N/A 1.883832 Upper 90.0% #N/A 2.134026 >> xi=0:.25:12; >> yi=spline(x,y,xi); >> plot(x,y,'o',xi,yi) >> spline(x,y,3) ans = 16.0669 19.18 Using Mathcad i 0.. 63 .2 .π .i sin .2 .π .i rnd x cos 3 10 () 1 .5 i 63 63 x i 0 0 5 i 2. π . 63 c j fft( x) 0.. 20 4 c j 2 0 0 10 20 j 19.19 As in Example 19.5, the data can be entered as >> x=[0 2 4 7 10 12]; >> y=[20 20 12 7 6 5.6]; Then, a set of x values can be generated and the interp1 function used to generate the linear interpolation >> xi=0:.25:12; >> yi=interp1(x,y,xi); These points can then be plotted with >> plot(x,y,'o',xi,yi) 20 15 10 5 0 2 4 6 8 10 12 The 5th-order interpolating polynomial and plot can be generated with >> p=polyfit(x,y,5) p = 0.0021 -0.0712 >> yi=polyval(p,xi); >> plot(x,y,'o',xi,yi) 0.8909 -4.5982 6.1695 20.0000 24 22 20 18 16 14 12 10 8 6 4 0 2 4 6 8 10 12 The cubic spline and plot can be generated with >> yi=spline(x,y,xi); >> plot(x,y,'o',xi,yi) 22 20 18 16 14 12 10 8 6 4 0 2 4 6 8 10 12 19.20 The following MATLAB session develops the fft along with a plot of the power spectral density versus frequency. >> >> >> >> >> >> t=0:63; y=cos(3*2*pi*t/63)+sin(10*2*pi*t/63)+randn(size(t)); Y=fft(y,64); Pyy=Y.*conj(Y)/64; f=1000*(0:31)/64; plot(f,Pyy(1:32)) 25 20 15 10 5 0 0 19.21 100 200 300 400 500 PROGRAM Fitpoly Use IMSL Implicit NONE Integer::ndeg,nobs,i,j Parameter (ndeg=4, nobs=6) Real:: b (ndeg + 1), sspoly(ndeg + 1), stat(10), X(nobs), y(nobs), ycalc (nobs) Data x/0,8,16,24,32,40/ Data y/14.621,11.843,9.870,8.418,7.305,6.413/ Call Rcurv(nobs, X, y, ndeg, b, sspoly, stat) Print *, 'Fitted polynomial is' Do i = 1,ndeg+1 Print 10, i - 1, b(i) End Do Print * Print 20, stat(5) Print * Print *, ' No. X Y YCALC' Do i = 1,nobs ycalc = 0 Do j = 1,ndeg+1 ycalc(i) = ycalc(i) + b(j)*x(i)**(j-1) End Do Print 30, i, X(i), y(i), ycalc(i) End Do 10 Format(1X, 'X^',I1,' TERM: ',F8.4) 20 Format(1X,'R^2: ',F8.3,'%') 30 Format(1X,I8,3(5X,F8.4)) End Output: Fitted polynomial is X^0 TERM: 14.6208 X^1 TERM: -0.4113 X^2 TERM: 0.0090 X^3 TERM: -0.0001 X^4 TERM: 0.0000 R^2: 100.000% No. 1 2 3 4 5 X 0.0000 8.0000 16.0000 24.0000 32.0000 Y 14.6210 11.8430 9.8700 8.4180 7.3050 YCALC 14.6208 11.8438 9.8685 8.4195 7.3042 6 40.0000 6.4130 6.4132 19.22 Using Excel, plot the data and use the trend line function to fit a polynomial of specific order. Obtain the R – squared value to determine the goodness of fit. Dye Concentraion vs. Time 5 y = 0.0008x3 - 0.0642x2 + 1.1609x - 2.6606 R2 = 0.8947 4 Dye Concentration 3 2 1 0 0 5 10 15 20 25 30 -1 Seconds after injection Dye Concentraion vs. Time 4.5 y = 0.0003x4 - 0.0151x3 + 0.2092x2 - 0.5741x + 0.3917 4 R2 = 0.9838 3.5 Dye Concentration 3 2.5 2 1.5 1 0.5 0 0 5 10 15 20 -0.5 Seconds after injection Use the 4th order polynomial: C = 0.0003t 4 − 0.0151t 3 + 0.2092t 2 − 0.5741t + 0.3917 Integrate to find the area under the curve: 25 30 24 ∫ 0 .0003t 4 − 0 .0151t 3 + 0.2092t 2 − 0.5741t + 0.3917 dt = 33. 225 2 Area under curve: 33.225 mg sec/L Cardiac output = 5 mg = 0.15049 L / sec = 9 L / min 33.225 mg sec/ L Cardiac output ≅ 9 L/min 19.23 Plug in Ao = 1 and T = 1 ⇒ 4 f (t ) = ∞   ∑  ( 2n −o1)π  sin  n =1 4A Make table and plot in Excel ⇒ Shown on the following pages 2π ( 2 n − 1)t   T  0 0 0 0 0 0 0 0 0.01 0.02 0.03 0.11 0.12 0.13 0.14 0.15 0.317 0.613 0.872 1.075 1.211 1.271 1.251 1.152 0.981 0.748 0.469 0.291 0.424 0.327 0.053 -0.249 -0.417 -0.358 -0.106 0.204 0.404 0.384 0.242 0.150 -0.150 -0.242 0.000 0.242 0.150 -0.150 -0.242 0.000 0.242 0.179 -0.067 -0.154 0.125 0.107 -0.165 -0.045 0.182 -0.023 -0.173 0.088 0.109 -0.139 0.068 0.052 -0.135 0.119 -0.018 -0.097 0.141 -0.083 -0.035 0.043 -0.079 0.105 -0.116 0.110 -0.089 0.056 -0.015 -0.029 0.068 -0.098 1.180 0.901 1.068 0.947 1.044 0.962 1.035 0.967 1.033 0.964 1.050 0.160 0.156 0.150 0.140 0.128 0.114 0.847 -0.160 -0.156 -0.150 -0.140 -0.128 -0.114 -0.847 -0.469 -0.384 -0.242 -0.088 0.035 0.098 -1.050 -0.748 -0.404 0.000 0.173 0.083 -0.068 -0.964 0.04 0.05 0.06 0.07 0.08 0.09 0.1 0.16 0.17 1.5 1 0.5 0 -0.5 -1 -1.5 0 0.1 0.2 0.3 0.4 0.5 Time, t 0.6 0.7 0.8 0.18 0.19 0.2 -0.981 -1.152 -1.251 -1.271 -1.211 -0.204 0.106 0.358 0.417 0.249 0.242 0.150 -0.150 -0.242 0.000 0.023 -0.182 0.045 0.165 -0.107 -0.141 0.097 0.018 -0.119 0.135 0.029 0.015 -0.056 0.089 -0.110 -1.033 -0.967 -1.035 -0.962 -1.044 Sum of the First Six Terms of the Fourier Series f(t) time-> n 1 2 3 4 5 6 Sum 0.9 1