Surfing in Bombay

> restart;

Preliminary definitions

> FwdDiff := proc(Y)
local k;
[seq(Y[k]-Y[k-1], k=2..nops(Y))]
end:#Define the FwdDiff function.

> SSE := proc(f,Xd,Yd) local k;
sum((Yd[k]-f(Xd[k]))^2, k=1..nops(Xd));
end:# Define the Sum of Squared Errors.

> SinReg2 := proc(w,Xd,Yd) local eq;
eq := stats[fit,leastsquare[[x,y],
y = a*cos(w*x)+b*sin(w*x)+d,
{a,b,d}]]([Xd,Yd]);
unapply(rhs(eq),x);
end:

> sdCount := proc(y) local fd,sd,ct,sg,k;
fd := FwdDiff(y);
sd := FwdDiff(fd);
ct := 0;
sg := sd[1];
for k from 2 to nops(sd) do
if sg*sd[k] < 0 then
sg := sd[k];
ct := ct+1;
fi;
od;
ct;
end:

> SinReg1 := proc(wlo,whi,Xd,Yd)
local k,w1,w2,w3,w4,rs1,rs2,rs3,rs4,
SS1,SS2,SS3,SS4,ww,rs,SS,q,a,b,c,ans,gf;
userinfo(3,SinReg,`wlo,whi:`,wlo,whi);
w1 := wlo;
rs1 := SinReg2(w1,Xd,Yd);
SS1 := SSE(rs1,Xd,Yd);
w4 := whi;
rs4 := SinReg2(w4,Xd,Yd);
SS4 := SSE(rs4,Xd,Yd);
gf := evalf((sqrt(5)-1)/2);
w3 := w1 + gf*(w4-w1);
rs3 := SinReg2(w3,Xd,Yd);
SS3 := SSE(rs3,Xd,Yd);
w2 := w1 + gf*gf*(w4-w1);
rs2 := SinReg2(w2,Xd,Yd);
SS2 := SSE(rs2,Xd,Yd);
userinfo(3,SinReg,"w1,w2,w3,w4:",w1,w2,w3,w4);
userinfo(3,SinReg,"SSE:",SS1,SS2,SS3,SS4);
q := x -> a*x^2 + b*x + c;
for k from 1 to 25 do
if ((SS2 < SS1) and (SS2 < SS3) and (SS3 < SS4)) then
w4 := w3;
rs4 := rs3;
SS4 := SS3;
ans := solve({q(w1)=SS1, q(w2)=SS2, q(w3)=SS3},{a,b,c});
ww := -subs(ans,b)/(2*subs(ans,a));
if ((w1 < ww) and (ww < w2)) then
w3 := w2;
rs3 := rs2;
SS3 := SS2;
w2 := ww;
rs2 := SinReg2(w2,Xd,Yd);
SS2 := SSE(rs2,Xd,Yd);
elif ((w2 < ww) and (ww < w3)) then
w3 := ww;
rs3 := SinReg2(w3,Xd,Yd);
SS3 := SSE(rs3,Xd,Yd);s;
else
w3 := w2;
rs3 := rs2;
SS3 := SS2;
w2 := w1 + gf*(w3-w1);
rs2 := SinReg2(w2,Xd,Yd);
SS2 := SSE(rs2,Xd,Yd);
fi;
elif (SS2 < SS1) and (SS3 < SS2) and (SS3 < SS4) then
w1 := w2;
rs1 := rs2;
SS1 := SS2;
ans := solve({q(w2)=SS2, q(w3)=SS3, q(w4)=SS4},{a,b,c});
ww := -subs(ans,b)/(2*subs(ans,a));
if ((w2 < ww) and (ww < w3)) then
w2 := ww;
rs2 := SinReg2(w2,Xd,Yd);
SS2 := SSE(rs2,Xd,Yd);
elif ((w3 < ww) and (ww < w4)) then
w2 := w3;
rs2 := rs3;
SS2 := SS3;
w3 := ww;
rs3 := SinReg2(w3,Xd,Yd);
SS3 := SSE(rs3,Xd,Yd);
else
w2 := w3;
rs2 := rs3;
SS2 := SS3;
w3 := w2 + gf*gf*(w4-w2);
rs3 := SinReg2(w3,Xd,Yd);
SS3 := SSE(rs3,Xd,Yd);
fi;
elif SS2 < SS3 then
w4 := w3;
rs4 := rs3;
SS4 := SS3;
w3 := w2;
rs3 := rs2;
SS3 := SS2;
w2 := w1 + gf*gf*(w4-w1);
rs2 := SinReg2(w2,Xd,Yd);
SS2 := SSE(rs2,Xd,Yd);
else
w1 := w2;
rs1 := rs2;
SS1 := SS2;
w2 := w3;
rs2 := rs3;
SS2 := SS3;
w3 := w1 + gf*(w4-w1);
rs3 := SinReg2(w3,Xd,Yd);
SS3 := SSE(rs3,Xd,Yd);
fi;
userinfo(3,SinReg,"w1,w2,w3,w4:",w1,w2,w3,w4);
userinfo(3,SinReg,"SSE:",SS1,SS2,SS3,SS4);
if (abs(w4-w1) < 0.001*w4) then
break;
fi;
od;
SS := min(SS1,SS2,SS3,SS4);
if (SS1 = SS) then
rs1;
elif (SS2 = SS) then
rs2;
elif (SS3 = SS) then
rs3;
else
rs4;
fi;
end:# Define the Sinusoidal Regression function

> SinReg := proc(Xd,Yd)
local k,cnt,w0,Sk,SS,kk,rs,A,B,a,b,c,d;
cnt := sdCount(Yd);
userinfo(3,SinReg,`Count:`,cnt);
if cnt < 2 then
cnt := 2;
fi;
# We assume that the data covers at least a quarter of the period.
w0 := evalf(0.5*Pi/abs(Xd[nops(Xd)]-Xd[1]));
kk := 1;
rs := SinReg2(kk*w0,Xd,Yd);
SS := SSE(rs,Xd,Yd);
userinfo(3,SinReg,`w, SS:`,w0, SS);
for k from 2 to 2*cnt+2 do
rs := SinReg2(k*w0,Xd,Yd);
Sk := SSE(rs,Xd,Yd);
userinfo(3,SinReg,`w, SS:`,k*w0, Sk);
if Sk < SS then
kk := k;
SS := Sk;
fi;
od;
if kk = 1 then
rs := SinReg1(w0,3*w0,Xd,Yd);
elif kk = 2*cnt+2 then
rs := SinReg1(2*cnt*w0,(2*cnt+2)*w0,Xd,Yd);
else
rs := SinReg1((kk-1)*w0,(kk+1)*w0,Xd,Yd);
fi;
userinfo(3,SinReg,`rs`,rs(x));
A := op(1,op(1,rs(x)));
B := op(1,op(2,rs(x)));
a := sqrt(A^2+B^2);
b := op(1,op(1,op(2,op(1,rs(x)))));
c := arctan(A,B);
d := op(3,rs(x));
unapply(a*sin(b*x+c)+d,x);
end:# Define the Sinusoidal Regression function

The Tide dat
a

> Hour := [ 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20];
Depth := [2.5, 2.7, 4.5, 7.5, 11.1, 14.5, 17.0, 18.0, 17.2, 15.0, 11.7, 8.0, 4.9, 2.9, 2.5];

> pts := [seq([Hour[k],Depth[k]],k=1..nops(Hour))]:
plot(pts,style=point);

> Bombay := SinReg(Hour,Depth);

> P1 := plot(pts,style=point,color=blue):
P2 := plot(Bombay(x),x=6..20,2..20):
plots[display](P1,P2);
SS5 := SSE(Bombay,Hour,Depth);

Zooming In

> pt := [[11,Bombay(11)]];
a := 8;
b := 14;
P1 := plot(pt,view=[8..14,6..18],style=point):
P2 := plot(Bombay(x),x=a..b):
plots[display](P1,P2);
(Bombay(b)-Bombay(a))/(b-a);

> pt := [[11,Bombay(11)]];
a := 10;
b := 12;
P1 := plot(pt,view=[10..12,11..17],style=point):
P2 := plot(Bombay(x),x=a..b):
plots[display](P1,P2);
(Bombay(b)-Bombay(a))/(b-a);

> pt := [[11,Bombay(11)]];
a := 10.9;
b := 11.1;
P1 := plot(pt,view=[10.9..11.1,14.2..14.8],style=point):
P2 := plot(Bombay(x),x=a..b):
plots[display](P1,P2);
(Bombay(b)-Bombay(a))/(b-a);

> pt := [[11,Bombay(11)]];
a := 10.99;
b := 11.01;
P1 := plot(pt,view=[10.99..11.01,14.49...14.56],style=point):
P2 := plot(Bombay(x),x=a..b):
plots[display](P1,P2);
(Bombay(b)-Bombay(a))/(b-a);