Fitting a Sine Curve to DataSteven R. DunbarDepartment of MathematicsUniversity of Nebraska-LincolnLincoln, NE, USA 68588sdunbar@math.unl.edu http://www.math.unl.edu/~sdunbar Copyright: 1999, 2005
<Text-field layout="Heading 2" style="Heading 2"><Font family="Times New Roman"> Physics Content</Font></Text-field>Making a best sine curve fit to a set of sparse data from observations of the star 51 Pegasi.Making a best sine curve fit to a set of sparse data from observation of the tides in the Bay of Fundy.
<Text-field layout="Heading 2" style="Heading 2"><Font family="Times New Roman">Mathematics Content</Font></Text-field>Curve fitting, and simple statistics, least-squares optimization
<Text-field layout="Heading 2" style="Heading 2"><Font family="Times New Roman">Keywords</Font></Text-field>best fit curve, least-squares fitting, sine curve fit, sparse data
<Text-field layout="Heading 2" style="Heading 2"><Font family="Times New Roman">Objective</Font></Text-field>Illustrate nonlinear curve fitting with Maple, using both elementary commands and sophisticated tools.
<Text-field layout="Heading 2" style="Heading 2"><Font family="Times New Roman">Audience</Font></Text-field>Anyone interested in using Maple to do simple curve curve fitting.
<Text-field layout="Heading 2" style="Heading 2"><Font family="Times New Roman">Prerequisite</Font></Text-field>Understanding of curve fitting problem.Understanding of mean of a set of data.Understanding of sine and arcsine functions.Understanding vertical, horizontal shifting of a graph, understanding change of scale of a graph.Understanding of least-squares line fitting.
<Text-field layout="Heading 1" style="Heading 1"><Font family="Times New Roman">Plan</Font></Text-field>We want to find the parameters A, f, phi, K, so that the curvey = A*sin(f*(x - phi)) + Kfits the set of data points { [x_i, y_i]} well.The mean of the curve is clearly K, so we can take the mean of the data to be K.The vertical range of the curve is 2*A, so we can take the vertical range of the vertically rescaled data to be 2*A, and solve for A.Then(y - K)/A = sin(f*(x - phi))and(1/f)* arcsin( (y-K)/A ) + phi = xso the line(1/f)*z + phi = xshould be a good fit to the set of data points { [arcsin( (y_i -K)/A, x_i)] }
restart;
<Text-field layout="Heading 1" style="Heading 1"><Font family="Times New Roman">Fitting a Sine Curve to Sparse Data from Observations of the Star 51 Pegasi</Font></Text-field>
<Text-field layout="Heading 2" style="Heading 2"><Font family="Times New Roman">Data</Font></Text-field>Velocity of the star 51 Pegasi measured in m/sVelocity := [2.65, 2.80, 2.96, 3.80, 3.90, 4.60, 4.80, 4.90, 5.65, 5.92];Days, measured in Julian DaysDays := [-45.5, -48.8, -54.0, -13.5, -7.0, 42.0, 50.5, 54.0, 36.0, 14.5];Pairs := zip((x,y) -> [x,y], Velocity, Days);plot( Pairs, symbol = DIAMOND, style = POINT);
<Text-field layout="Heading 2" style="Heading 2"><Font family="Times New Roman">Mean of the Data, the vertical translation of the data, range</Font></Text-field>with(stats);K := describe[mean](Days);RescaledDays := map( y -> y - K, Days);ScaledDayRange := describe[range](RescaledDays);A := max( abs(op(1,ScaledDayRange)), abs(op(2, ScaledDayRange)));
<Text-field layout="Heading 2" style="Heading 2"><Font family="Times New Roman">Transform the Data</Font></Text-field>ArcsinScaledDays := map( y -> arcsin( (y-K)/A), Days);
<Text-field layout="Heading 2" style="Heading 2"><Font family="Times New Roman">Fitting Transformed Data</Font></Text-field>FittedLine := fit[leastsquare[[z,x], x=RecipFreq*z+phi]]([ArcsinScaledDays,Velocity]);f := 1/coeff(rhs(FittedLine), z, 1); phi := coeff(rhs(FittedLine), z, 0);DataPlot := plot( Pairs, symbol = DIAMOND, style = POINT):CurvePlot := plot( A*sin(f*(x-phi)) + K, x=Velocity..Velocity[nops(Velocity)], color = blue):plots[display](DataPlot, CurvePlot);RMSError := sqrt( sum( ( A*sin(f*(Velocity[i]-phi))+K - Days[i])^2, i = 1..nops(Velocity) ) );
<Text-field layout="Heading 2" style="Heading 2"><Font family="Times New Roman">Some Additional Fitting by Eye</Font></Text-field>f := 1.5; CurvePlot := plot( A*sin(1.5*(x-phi)) + K, x=Velocity..Velocity[nops(Velocity)], color = blue):plots[display](DataPlot, CurvePlot);RMSError := sqrt( sum( ( A*sin(f*(Velocity[i]-phi))+K - Days[i])^2, i = 1..nops(Velocity) ) ); phi := 4.0; CurvePlot := plot( A*sin(f*(x-phi)) + K, x=Velocity..Velocity[nops(Velocity)], color = blue):plots[display](DataPlot, CurvePlot);RMSError := sqrt( sum( ( A*sin(f*(Velocity[i]-phi))+K - Days[i])^2, i = 1..nops(Velocity) ) );
<Text-field layout="Heading 2" style="Heading 2"><Font family="Times New Roman">Curve Fitting by Least-Squares (Nonlinear) Optimization</Font></Text-field>with(Optimization);NonlinearResiduals := [ seq( A0*sin(f0*(Velocity[i]-phi0))+K0 - Days[i], i = 1..nops(Velocity)) ]:lssoln := LSSolve( NonlinearResiduals, initialpoint={A0=A, f0=f, phi0=phi, K0=K});CurvePlot := plot( subs( lssoln, A0*sin(f0*(x-phi0)) + K0), x=Velocity..Velocity[nops(Velocity)], color = blue):plots[display]( DataPlot, CurvePlot);RMSError := sqrt( sum( ( subs( lssoln, A0*sin(f0*(Velocity[i]-phi0))+K0) - Days[i])^2, i = 1..nops(Velocity) ) );The RMSError and the objective function measured by the Least-Squares minimization are related by Least-Squares = (1/2)*RMSError^2: lssoln, (1/2)*RMSError^2;
<Text-field layout="Heading 1" style="Heading 1"><Font family="Times New Roman">Fitting a Sine Curve to Sparse Data from observations of the Tides in the Bay of Fundy</Font></Text-field>restart;The Bay of Fundy in Nova Scotia has some of the largest tides in the world: the water level can rise and fall 50 feet in one day. Nova Scotia bends when the tide comes in! As 14 billion tonnes (14 cubic kilometers) of sea water flow into Minas Basin twice daily, the Nova Scotia countryside actually tilts slightly under the immense load. In this project I will analyze some data concerning those tides.The table below gives times and heights for low and high tides at Cape Blomidon. The heights are measured in meters and represent the depth of the water compared to a fixed "standard" waterline.From the data, find a function that models the water depth at Cape Blomidon as a function of time. The model should cover the entire range of time covered by the data. Make clear the units and the meaning of the independent variable, and the function values. NiNJJURhdGVHNiI=NiNJJERheUc2Ig==NiMqJkklSGlnaEc2IiIiIkklVGltZUdGJUYmNiMqJkklSGlnaEc2IiIiIkknSGVpZ2h0R0YlRiY=NiMqJkkkTG93RzYiIiIiSSVUaW1lR0YlRiY=NiMqJkkkTG93RzYiIiIiSSdIZWlnaHRHRiVGJg==NiMqJkklSGlnaEc2IiIiIkklVGltZUdGJUYmNiMqJkklSGlnaEc2IiIiIkknSGVpZ2h0R0YlRiY=NiMqJkkkTG93RzYiIiIiSSVUaW1lR0YlRiY=NiMqJkkkTG93RzYiIiIiSSdIZWlnaHRHRiVGJg==NiMsJEkkU2VwRzYiIiNHNiNJJFR1ZUc2Ig==NiMiJUM+NiMkIiMpKSEiIw==NiMsJEkkU2VwRzYiIiNINiNJJFdlZEc2Ig==NiMiJFEiNiMkIiVyNyEiIw==NiMiJFooNiMkIiMpKSEiIw==NiMiJSs5NiMkIiVqNyEiIw==NiMiJTI/NiMkIiN5ISIjNiMsJEkkU2VwRzYiIiNJNiNJJFRodUc2Ig==NiMiJEEjNiMkIiVmNyEiIw==NiMiJEcpNiMkIiQxIiEiIw==NiMiJVQ5NiMkIiVkNyEiIw==NiMiJVw/NiMkIiMpKSEiIw==NiNJJE9jdEc2Ig==NiNJJEZyaUc2Ig==NiMiJC4kNiMkIiVLNyEiIw==NiMiJDMqNiMkIiRTIiEiIw==NiMiJT86NiMkIiVONyEiIw==NiMiJUlANiMkIiQ6IiEiIw==NiMsJEkkT2N0RzYiIiIjNiNJJFNhdEc2Ig==NiMiJFgkNiMkIiUlPiIhIiM=NiMiJFsqNiMkIiQlPSEiIw==NiMiJSw7NiMkIiUtNyEiIw==
<Text-field layout="Heading 2" style="Heading 2"><Font family="Times New Roman">Data</Font></Text-field>TidalHeights := [0.88, 12.71,0.88,12.63,0.78,12.59,1.06,12.57,0.88,12.32, 1.40, 12.35, 1.15,11.94, 1.84, 12.02];Times, measured from midnight, 28 SeptemberTidalTimes := map( evalf, [19+24/60, 24+1+38/60, 24+7+47/60, 24+14, 24+20+7/60, 48+2+22/60, 48+8+28/60, 48+14+41/60, 48+20+49/60, 72+3+3/60, 72+9+8/60, 72+15+20/60, 72+21+30/60, 96+3+45/60, 96+9+48/60, 96+16+1/60]); Pairs := zip((x,y) -> [x,y], TidalTimes, TidalHeights);plot( Pairs, symbol = DIAMOND, style = POINT);
<Text-field layout="Heading 2" style="Heading 2"><Font executable="false" family="Times New Roman">Mean of the Data, the vertical translation of the data, range</Font></Text-field>with(stats);K := describe[mean](TidalHeights);RescaledTidalHeights := map( y -> y - K, TidalHeights);ScaledTidalHeightsRange := describe[range](RescaledTidalHeights);A := max( abs( op(1,ScaledTidalHeightsRange)), abs( op(2, ScaledTidalHeightsRange)));
<Text-field layout="Heading 2" style="Heading 2"><Font family="Times New Roman">Transform the Data</Font></Text-field>ArcsinScaledTidalHeights := map( y -> arcsin( (y-K)/A), TidalHeights);
<Text-field layout="Heading 2" style="Heading 2"><Font family="Times New Roman">Fitting Transformed Data</Font></Text-field>FittedLine := fit[leastsquare[[z,x], x=RecipFreq*z+phi]]([ArcsinScaledTidalHeights,TidalTimes]);f := 1/coeff(rhs(FittedLine), z, 1); phi := coeff(rhs(FittedLine), z, 0);DataPlot := plot( Pairs, symbol = DIAMOND, style = POINT):CurvePlot := plot( A*sin(f*(x-phi)) + K, x=TidalTimes..TidalTimes[nops(TidalTimes)], color = blue):plots[display](DataPlot, CurvePlot);RMSError := sqrt( sum( ( A*sin(f*(TidalTimes[i]-phi))+K - TidalHeights[i])^2, i = 1..nops(TidalTimes) ) );
<Text-field layout="Heading 2" style="Heading 2"><Font executable="false" family="Times New Roman">Some Additional Fitting by Eye</Font></Text-field>AverageHigh := describe[mean]([ seq(TidalHeights[2*i], i = 1..nops(TidalHeights)/2 )]); AverageLow := describe[mean]([ seq(TidalHeights[2*i-1], i = 1..nops(TidalHeights)/2 )]); A := (AverageHigh - AverageLow)/2; K := (AverageHigh + AverageLow)/2;period := describe[mean]( [ seq( TidalTimes[i+2]-TidalTimes[i], i=1..nops(TidalTimes)-2) ]);f := evalf(2*Pi/period);phi := solve( f*(TidalTimes-x) = -Pi/2, x);CurvePlot := plot( A*sin(f*(x-phi)) + K, x=TidalTimes..TidalTimes[nops(TidalTimes)], color = blue):plots[display](DataPlot, CurvePlot);RMSError := sqrt( sum( ( A*sin(f*(TidalTimes[i]-phi))+K - TidalHeights[i])^2, i = 1..nops(TidalTimes) ) );
<Text-field layout="Heading 2" style="Heading 2"><Font executable="false" family="Times New Roman">Curve Fitting by Least-Squares (Nonlinear) Optimization</Font></Text-field>with(Optimization);A;f;phi;K;NonlinearResiduals := [ seq( A0*sin(f0*(TidalTimes[i]-phi0))+K0 - TidalHeights[i], i = 1..nops(TidalTimes)) ]:lssoln := LSSolve( NonlinearResiduals, initialpoint={A0=A, f0=f, phi0=phi, K0=K});CurvePlot := plot( subs( lssoln, A0*sin(f0*(x-phi0)) + K0), x=TidalTimes..TidalTimes[nops(TidalTimes)], color = blue):plots[display](DataPlot,CurvePlot);RMSError := sqrt( sum( ( subs( lssoln, A0*sin(f0*(TidalTimes[i]-phi0))+K0 - TidalHeights[i])^2), i = 1..nops(TidalTimes) ) );The RMSError and the objective function measured by the Least-Squares minimization are related by Least-Squares = (1/2)*RMSError^2: lssoln, (1/2)*RMSError^2;