DE.03 Laplace Transform and Fourier Analysis
TUTORIALS

\(Needs["Graphics`Colors`"]; 

    Graphics`Colors`GosiaGreen = RGBColor[0, \ 0.392187, \ 0]; 

    Needs["Graphics`FilledPlot`"]; 


    If[MemberQ[{"Global`", "Graphics`Common`GraphicsCommon`", 
          "Vector3D`"}, \ Context[ScaleFactor]], 
      Remove["*`ScaleFactor"]]; 
<< "Graphics`Legend`"; 
    
<< "Graphics`Common`GraphicsCommon`"; 

    Remove["*`ScaleFactor"]; 

Off[General::"spell"]; 

    \(Off[General::"spell1"]; \)\)Off[ParametricPlot3D::"ppcom"]; Off[ParametricPlot::"ppcom"]; \(Off[Plot::"plnr"]; 
Needs["Calculus`Master`"]; 

    Off[Unique::"usym"]; 
SetOptions[Limit, Analytic -> True]; 
    \)

CMView = {2.7, 1.6, 1.2}; \(Clear[FastFourierfit, Fourierfitters, F, Fvalues, n, k, jump, num, 
      numtab, coeffs, t, L]; 
jump[n_] := \(jump[n] = N[1\/\(2\ n\)]\); 

    \(Fvalues[F_, L_, n_] := N[Table[F[L\ t], {t, 0, 1 - jump[n], jump[n]}]]; 
    \)\)numtab[n_] := numtab[n] = Table[k, {k, 1, n}]; Fourierfitters[L_, n_, t_] := 
      Table[(E^(((2 π I k t)/L))), {k, -n + 1, n - 1}]; coeffs[n_, list_] := 
      (1/N)[(Sqrt[Length])[list]] 
        (Join[Reverse[
              Fourier[list][[numtab[n]
                  ]]], 
            InverseFourier[list][[Drop[numtab[n], 1]
                ]]])FastFourierfit[F_, L_, n_, t_] := 
      Chop[Fourierfitters[L, n, t] . coeffs[n, Fvalues[F, L, n]]]; \(L := 
      Expand[a\ ∂\_{#2, 2}#1 + b\ ∂\_#2 #1 + c\ #1]&\)

Off[Graphics::"gpt"]

New Vector Drawers

BeginPackage["Vector3D`", "Graphics`Colors`"];

\(Axes3D::"usage" = 
      "Axes3D[a, b] creates a Graphics3D object of cartesian axes with x, y, and z running from -a/3 to a, and with axes labels b units beyond the tips of the axes.  Axes3D[a] is Axes3D[a, a/8]."; 

    \(Perpend::"usage" = 
      "Perpend[a] returns a unit vector which is perpendicular to the vector a."; \)\)Arrow::"usage" = 
      "Arrow[a] produces a 2 or 3 dimensional vector from the origin to a. \
Arrow[a, Tail -> tail] gives a vector from tail to a."; Vector::"usage" = 
      "Vector[a, b] gives a 2D or 3D vector from a to b."; ArrowHead::"usage" = 
      "ArrowHead[a, vec] produces an arrowhead with its tip placed at the \
point a, pointing in the direction of the vector vec."; Tail::"usage" = 
      "Tail -> point puts the tail of the vector at point."; Norm::"usage" = 
      "Norm[vec] returns the Euclidean length of vec."; Normalize::"usage" = 
      "Normalize[vec] returns a unit vector in the direction of vec, \
provided vec is nonzero."; Aperture::"usage" = 
      "Aperture is the ratio of the radius of the base to the length of the \
head of a vector."; TipSize::"usage" = 
      "TipSize specifies an absolute size for the inner tip length of the \
head of a vector."; TipRatio::"usage" = 
      "TipRatio is the ratio of the inner length of the arrowhead to the \
length of the shaft of a vector."; HeadRatio::"usage" = 
      "HeadRatio is the ratio of the outer length of the arrowhead to the \
length of the shaft of a vector."; HeadSize::"usage" = 
      "HeadSize specifies an absolute size for the head of a vector."; 
    ScaleFactor::"usage" = 
      "ScaleFactor specifies the amount to scale a vector in length and may \
be either a number or a pure function."; ZeroVectorPointSize::"usage" = 
      "ZeroVectorPointSize is the size of the point used to represent a \
zero vector."; Options[Vector] = 
      Options[Arrow] = 
        Options[ArrowHead] = {HeadRatio -> 0.18, 
            TipRatio -> 0.14, Aperture -> 0.3, 
            ScaleFactor -> 1, HeadType -> Polygon, 
            ShaftQ -> True, EdgesQ -> True, 
            VectorColor -> RGBColor[0, 0, 1], ShaftWidth -> 0.005, 
            ZeroVectorPointSize -> 0.01}; Begin["`Private`"]; SetOptions[ParametricPlot, AspectRatio -> Automatic]; SetOptions[Plot, AspectRatio -> Automatic]; SetOptions[Graphics, AspectRatio -> Automatic]; Axes3D[u_, v_] := 
      Graphics3D[{{Blue, Line[{{-(u/3), 0, 0}, {u, 0, 0}}]}, 
          Text["x", {u + v, 0, 0}], {Blue, 
            Line[{{0, -(u/3), 0}, {0, u, 0}}]}, 
          Text["y", {0, u + v, 0}], {Blue, 
            Line[{{0, 0, -(u/3)}, {0, 0, u}}]}, 
          Text["z", {0, 0, u + v}]}]; Axes3D[u_] := Axes3D[u, (u/8)]; Norm[a_?VectorQ] := N[(Sqrt[a . a])]; Normalize[a_?VectorQ] := If[a . a == 0, a, (a/Norm)[a]]; Perpend[a_?VectorQ] := 
      Normalize[{a[[2]], 
            -a[[1]]}] /; 
        Length[a] == 2Perpend[a_?VectorQ] := 
      Normalize[
          If[a[[1]] == 0, {1, 0, 0}, {0, 
              a[[3]], 
              -a[[2]]}]] /; 
        Length[a] == 3Base3D = 
      Table[{Cos[((2 k π)/8.) ], Sin[((2 k π)/8.) ]}, {k, 0, 
          8}]; Vector[from_?VectorQ, to_?VectorQ, opts___] := 
      Module[{sq, sw, ht, hs, hr, tr, ap, sf, vc, zvps, diff = to - from, 
            tip}, {sq, sw, ht, hs, hr, ts, tr, ap, sf, vc, zvps} = 
            {ShaftQ, ShaftWidth, HeadType, HeadSize, HeadRatio, TipSize, 
                  TipRatio, Aperture, ScaleFactor, VectorColor, 
                  ZeroVectorPointSize} /. {opts} /. Options[Vector]; 
          If[from == to, 
            Return[Graphics[{PointSize[zvps], vc, Point[from]}]]]; 
          diff = If[NumberQ[N[sf]], sf diff, sf[diff]]; tip = from + diff; 
          If[NumberQ[hs], hr = (hs/Norm)[diff]; tr = 0.75 hr]; 
          If[NumberQ[ts], tr = (ts/Norm)[diff]]; 
          Graphics[{
              If[sq, {vc, Thickness[sw], Line[{from, tip - tr diff}]}, {}], 
              With[{trans = ap hr Norm[diff] Perpend[diff]}, {vc, 
                  Thickness[sw], 
                  ht[{tip, tip - hr diff + trans, tip - tr diff, 
                      tip - hr diff - trans, tip}]}]}]] /; 
        Length[from] == Length[to] == 2Vector[from_?VectorQ, to_?VectorQ, opts___] := 
      Module[{sq, eq, sw, hs, hr, tr, ap, sf, vc, zvps, diff = to - from, 
            tip}, {sq, eq, sw, hs, hr, ap, sf, vc, zvps} = 
            {ShaftQ, EdgesQ, ShaftWidth, HeadSize, HeadRatio, Aperture, 
                  ScaleFactor, VectorColor, ZeroVectorPointSize} /. {opts}
                   /. Options[Vector]; 
          If[from == to, 
            Return[Graphics3D[{PointSize[zvps], vc, Point[from]}]]]; 
          diff = If[NumberQ[N[sf]], sf diff, sf[diff]]; tip = from + diff; 
          If[NumberQ[hs], hr = (hs/Norm)[diff]]; tr = (hr/2); 
          Graphics3D[{
              If[sq, {vc, Thickness[sw], Line[{from, tip - tr diff}]}, {}], 
              If[eq, {}, EdgeForm[]], {SurfaceColor[vc], 
                With[{temp = Perpend[diff]}, 
                  (Polygon[Append[#1, tip]]&)/@
                    Partition[
                      (tip - hr diff + #1&)/@
                        ((ap hr Norm[diff] Base3D) . {temp, 
                              Normalize[temp\[Cross]diff]}), 2, 1]]}}]] /; 
        Length[from] == Length[to] == 3Arrow[a_, Tail -> b_, opts___] := Vector[b, a + b, opts]Arrow[a_, opts___] := Vector[Table[0, {Length[a]}], a, opts]\(ArrowHead[a_, b_, opts___] := 
      Vector[a - b, a, ShaftQ -> False, opts]; 

End[]; 

    \(EndPackage[]; \)\)

VectorField and Phase Lines

BeginPackage["PhaseStuff`", "Vector3D`", 
      "Utilities`FilterOptions`"]; PhaseLine::"usage" = 
      "PhaseLine[fcn,{var,a,b},tint,pos] makes a vertical phase line for \
the function, fcn, of the variable, var, for the variable going from a to b, \
with  a color given by tint,  at the x-coordinate pos."; PhaseLine3D::"usage" = 
      "PhaseLIne3D[fcn,{var,a,b},tint,loc] makes a vertical phase line for \
the function, fcn, of the variable, var, for the variable going from a to b, \
with  a color given by tint, and base point at loc = {x,y}"; VectorField::"usage" = 
      "VectorField[{function,{var1,a,b,space1},{var2,c,d,space2}] will make \
a vector field of normalized vectors, colored by magnitude, for the function \
of the variables on the range {r1,r2} with spacing between points given by  \
{space1, space2}.  If a single function of a single variable is given, the \
vector field {1,f} will be plotted."; StikPikXZ::"usage" = 
      "StikPikXZ[picture,a:_0] takes a 2 dimensional Mathematica picture \
and embeds it on the XZ plane at y = a.   Some optional characteristics of \
the graphs are lost in the  process.  StikPikYZ and StikPikXZ are \
similar."; StikPikYZ::"usage" = 
      "StikPikYZ[picture,a:_0] takes a 2 dimensional Mathematica picture \
and embeds it on the XZ plane at y = a.   Some optional characteristics of \
the graphs are lost in the  process.  StikPikYZ and StikPikXZ are \
similar."; StikPikXY::"usage" = 
      "StikPikXY[picture,a:_0] takes a 2 dimensional Mathematica picture \
and embeds it on the XZ plane at y = a.   Some optional characteristics of \
the graphs are lost in the  process.  StikPikYZ and StikPikXZ are \
similar."; Begin["`Private`"]; PhaseLine[fcn_, {var_, a_, b_}, tint_, pos_:  0] := 
      Module[{len, len1, k, zset, zero, zero1, signs, stick, dots, pointers}, 
        zset = Solve[fcn == 0, var]; 
        zero = If[Length[zset] == 0, {}, var /. zset]; zero = Sort[zero]; 
        zero = Select[zero, #1 ≥ a && #1 ≤ b&]; 
        zero1 = Sort[zero ⋃ {a, b}]; len = Length[zero]; 
        len1 = Length[zero1]; 
        signs = Table[
            Sign[fcn /. 
                var -> 
                  (1/2) (
                      zero1[[k + 1]] + 
                        zero1[[k]])], {
              k, 1, len1 - 1}]; 
        stick = {Thickness[0.007], Line[{{pos, a}, {pos, b}}]}; 
        dots = (Point[{pos, #1}]&)/@zero; 
        dots = Join[{PointSize[0.05], tint}, dots]; 
        heads = Table[
            MakeHead[zero1[[k]], 
              zero1[[k + 1]], pos, 
              signs[[k]], b - a, tint], {
              k, 1, len1 - 1}]; 
        Return[{Graphics[Join[{tint}, stick, dots]], heads}]]MakeHead[a_, b_, c_, s_, slen_, tint_] := 
      Module[{u, v, w}, u = Min[2 (b - a), ((3 slen)/4)]; 
        Return[ArrowHead[{c, a ( .5 - (s/5)) + b ( .5 + (s/5))}, {0, 
              s u}, VectorColor -> tint]]]PhaseLine3D[fcn_, {var_, a_, b_}, tint_, pos_: {0, 0}] := 
      Module[{len, len1, k, zset, zero, zero1, signs, stick, dots, pointers}, 
        zset = Solve[fcn == 0, var]; 
        zero = If[Length[zset] == 0, {}, var /. zset]; zero = Sort[zero]; 
        zero = Select[zero, #1 ≥ a && #1 ≤ b&]; 
        zero1 = Sort[zero ⋃ {a, b}]; len = Length[zero]; 
        len1 = Length[zero1]; 
        signs = Table[
            Sign[fcn /. 
                var -> 
                  (1/2) (
                      zero1[[k + 1]] + 
                        zero1[[k]])], {
              k, 1, len1 - 1}]; 
        stick = {Thickness[0.007], 
            Line[{{pos[[1]], 
                  pos[[2]], a}, {
                  pos[[1]], 
                  pos[[2]], b}}]}; 
        dots = 
          (Point[{pos[[1]], 
                    pos[[2]], #1}]&)/@
            zero; dots = Join[{PointSize[0.03], tint}, dots]; 
        heads = Table[
            MakeHead3D[zero1[[k]], 
              zero1[[k + 1]], pos, 
              signs[[k]], b - a, tint], {
              k, 1, len1 - 1}]; 
        Return[{Graphics3D[Join[{tint}, stick, dots]], heads}]]MakeHead3D[a_, b_, c_, s_, slen_, tint_] := 
      Module[{u, v, w}, u = Min[2 (b - a), (slen/2)]; 
        Return[ArrowHead[{c[[1]], 
              c[[2]], 
              a ( .5 - (s/5)) + b ( .5 + (s/5))}, {0, 0, s u}, 
            VectorColor -> tint]]]VectorField[
        function_List, {var1_, a_, b_, sp1_:  1}, {var2_, c_, d_, 
          sp2_:  1}] := 
      Module[{i, j, U, V, UU, VV, ulow, uhigh, vlow, vhigh, maxval, minval, 
          lengths, directions, locations, mlen}, mlen = Min[sp1, sp2]; 
        normf[uu_, vv_] = 
          (Sqrt[function[[1]]^2 + 
                function[[2]]^2]) /. {
              var1 -> uu, var2 -> vv}; 
        lengths = Table[N[normf[U, V]], {U, a, b, sp1}, {V, c, d, sp2}]; 
        maxval = Max[lengths]; minval = Min[lengths]; 
        If[minval == maxval, maxval += 1]; 
        locations = Table[{uu, vv}, {uu, a, b, sp1}, {vv, c, d, sp2}]; 
        directions = 
          Table[If[N[normf[U, V]] == 0, {0, 0}, 
              ((mlen (function /. {var1 -> U, var2 -> V})
                    )/(1.3 N[normf[U, V]]))], {U, a, b, sp1}, {V, c, d, 
              sp2}]; {UU, VV} = 
          Dimensions[directions][[{1, 2}
              ]]; 
        Return[Table[
            Arrow[directions[[i, j]], 
              Tail -> 
                locations[[i, j]], 
              Aperture -> 1, HeadRatio ->  .4, 
              ShaftWidth -> 0.008, 
              VectorColor -> 
                RGBColor[
                  ((lengths[[i, j]] - 
                      minval)/(maxval - minval)), 0, 
                  ((maxval - 
                      lengths[[i, j
                        ]])/(maxval - minval))]], {i, 
              UU}, {j, VV}]]]VectorField[
        function_, {var1_, a_, b_, sp1_:  1}, {var2_, c_, d_, sp2_:  1}] := 
      VectorField[{1, function}, {var1, a, b, sp1}, {var2, c, d, sp2}]Needs["Utilities`FilterOptions`"]FixOptions[a_] := 
      Module[{picture = a, opt3}, 
        If[Options[picture] ≠ {}, 
          opt3 = (FilterOptions[Graphics3D, #1]&)/@Options[picture]; 
          opt3 = opt3 //. {
                (AxesLabel -> _) -> AxesLabel -> None, 
                (PlotRange -> _) -> 
                  PlotRange -> Automatic, 
                (PlotRegion -> _) -> 
                  PlotRegion -> Automatic, 
                (Ticks -> _) -> Ticks -> Automatic}; 
          picture = List@@picture; 
          picture[[-1]] = opt3; 
          picture = Graphics3D@@picture]; Return[picture]]StikPikXZ[picture_, a_:  0] := 
      Module[{PIK3 = picture}, 
        PIK3 = PIK3 //. {
              Graphics -> Graphics3D, {u_?NumberQ, v_?NumberQ} 
                :> {u, a, v}}; PIK3 = FixOptions[PIK3]; 
        Return[PIK3]]StikPikYZ[picture_, a_:  0] := 
      Module[{PIK3 = picture}, 
        PIK3 = PIK3 //. {
              Graphics -> Graphics3D, {u_?NumberQ, v_?NumberQ} 
                :> {a, u, v}}; PIK3 = FixOptions[PIK3]; 
        Return[PIK3]]StikPikXY[picture_, a_:  0] := 
      Module[{PIK3 = picture}, 
        PIK3 = PIK3 //. {
              Graphics -> Graphics3D, {u_?NumberQ, v_?NumberQ} 
                :> {u, v, a}}; PIK3 = FixOptions[PIK3]; 
        Return[PIK3]]End[]; EndPackage[];

ThreeAxes[u,v]

ThreeAxes[u_, v_] := 
      Graphics3D[{{Blue, Line[{{-u, 0, 0}, {u, 0, 0}}]}, 
          Text["x", {u + v, 0, 0}], {Blue, 
            Line[{{0, -u, 0}, {0, u, 0}}]}, 
          Text["y", {0, u + v, 0}], {Blue, Line[{{0, 0, 0}, {0, 0, u}}]}, 
          Text["z", {0, 0, u + v}]}]; ThreeAxes[u_] := ThreeAxes[u, (u/8)]; ThreeAxes::"usage" = 
      "ThreeAxes[a,b] makes a standard cartesian axis graphics object with \
x, y, and z running from -a to a, and with axis labels b units beyond the \
tips of the axes.  ThreeAxes[a] is ThreeAxes[a,a/8].";

Crisscross

\(Clear[x]; 

    crisscross = 
      Plot[{Cos[x], Sin[x]}, {x, π\/8, \(3\ π\)\/8}, 
        PlotStyle -> {{RGBColor[1, 0, 0], Thickness[0.02]}}, 
        Axes -> None, AspectRatio -> 1\/3, 
        DisplayFunction -> Identity]; \)

Pendulum Stuff

modticksX[gr_Graphics, m_] := 
      Module[{f, x}, f[{a_, "", st___}] := {a, ""}; 
        f[{a_, b_, st___}] := {a, b m}; 
        Show[gr, AxesOrigin -> Reverse[FullOptions[gr, AxesOrigin]], 
          Ticks -> {
              f/@FullOptions[gr, Ticks][[1
                  ]], Automatic}]]

modticksY[gr_Graphics, m_] := 
      Module[{f, x}, f[{a_, "", st___}] := {a, ""}; 
        f[{a_, b_, st___}] := {a, b m}; 
        Show[gr, AxesOrigin -> Reverse[FullOptions[gr, AxesOrigin]], 
          Ticks -> {Automatic, 
              f/@FullOptions[gr, Ticks][[2
                  ]]}]]

Clear[center, angle, length, a, l, pendulum, calibratedpendulum];

pendulum[length_, angle_] := 
      Module[{center := {length Sin[π - angle], 
              length Cos[π - angle]}}, 
        Show[Graphics[{Blue, Thickness[ .01], 
              Line[{{0, 0}, {0, -(length + 3)}}]}], 
          Graphics[{Blue, 
              Polygon[{{-5, -(length + 3)}, {-5, 
                    -(length + 3.5)}, {-5.3, 
                    -(length + 3.5)}, {-5.3, 
                    -(length + 4)}, {5.3, -(length + 4)}, {
                    5.3, -(length + 3.5)}, {5, 
                    -(length + 3.5)}, {5, -(length + 3)}, {
                    -5, -(length + 3)}}]}], 
          Graphics[{Black, 
              Disk[{0, 0}, 
                length, {-(π/2) + angle - (π/84), 
                  -(π/2) + angle + (π/84)}]}], 
          Graphics[{Black, 
              Disk[center, 
                length, {(π/2) + angle - (π/500), 
                  (π/2) + angle + (π/500)}]}], 
          Graphics[{Red, PointSize[ .05], Point[{0, 0}]}], 
          Graphics[{Red, Disk[center, (length/5), {0, 360 Degree}]}], 
          Graphics[{GoldenrodLight, 
              Disk[center, (length/6), {330 Degree, 350 Degree}]}], 
          Graphics[{Red, Disk[center, (length/7), {0, 360 Degree}]}], 
          Graphics[{Black, PointSize[ .02], Point[center]}], 
          AspectRatio -> Automatic, 
          PlotRange 
            -> {{length + 5, -(length + 3)}, {length + 4, 
                -(length + 4)}}, AxesOrigin -> {0, 0}, 
          DisplayFunction -> Identity]]

calibratedpendulum[length_, angle_]:=Module[{center := {length Sin[π - angle], length Cos[π - angle]}}
            ,Show[Graphics[{Blue, Thickness[ .01], 
                  Line[{{0, 0}, {0, -(length + 3)}}]}],Graphics[{Blue, 
                  Polygon[{{-5, -(length + 3)}, {-5, 
                        -(length + 3.5)}, {-5.3, 
                        -(length + 3.5)}, {-5.3, 
                        -(length + 4)}, {5.3, 
                        -(length + 4)}, {5.3, 
                        -(length + 3.5)}, {5, 
                        -(length + 3.5)}, {5, 
                        -(length + 3)}, {-5, 
                        -(length + 3)}}]}],Graphics[{Black, 
                  Disk[{0, 0}, 
                    length, {-(π/2) + angle - (π/84), 
                      -(π/2) + angle + (π/84)}]}],Graphics[{Black, 
                  Disk[center, 
                    length, {(π/2) + angle - (π/500), 
                      (π/2) + angle + (π/500)}]}],Graphics[{Red, PointSize[ .05], Point[{0, 0}]}],Graphics[{Red, Disk[center, (length/5), {0, 360 Degree}]}],Graphics[{GoldenrodLight, 
                  Disk[center, (length/6), {330 Degree, 350 Degree}]}],Graphics[{Red, 
                  Disk[center, (length/7), {0, 360 Degree}]}],Graphics[{Black, PointSize[ .02], Point[center]}],Graphics[Text["y[t] =  \!\(angle\)",{(8/8) length Sin[((15 π)/12)], 
                      (8/8) length Cos[((5 π)/6)]}]],AspectRatio -> Automatic,PlotRange 
                -> {{length + 5, -(length + 3)}, {length + 4, 
                    -(length + 4)}},AxesOrigin -> {0, 0},DisplayFunction -> Identity]]

T.1) Using Fourier Fit to try to detect periodic forcing functions that force some undamped oscillators to beat or go into near or true resonance

T.1.a.i)

Activate this code.

\(Clear[FastFourierfit, Fourierfitters, F, Fvalues, n, k, 
\ \ \ \ \ \ 
      jump, num, numtab, coeffs, t, L]; 

    jump[n_]\  := \ \(jump[n]\  = \ N[1\/\(2\ n\)]\); 

    Fvalues[F_, L_, n_]\  := \ 
\ \ 
      N[Table[F[L\ t], {t, 0, 1\  - \ jump[n], jump[n]}]]; 
\ \ 


    numtab[n_]\  := \ \(numtab[n]\  = \ Table[k, {k, 1, n}]\); 


    Fourierfitters[L_, n_, t_]\  := \ 
      Table[E\^\(\(2\ π\ I\ k\ t\)\/L\), 

        \ \ \ \ \ \ \ \ \ \ {k, \(-n\)\  + \ 1, n\  - \ 1}]; 

    coeffs[n_, list_]\  := \ 
      Join[Reverse[Part[Fourier[list], numtab[n]]], 
	
          Part[InverseFourier[list], Drop[numtab[n], 1]]]/
		
        N[Sqrt[Length[list]]]
\)\(\(FastFourierfit[F_, L_, n_, t_]\  := \ 
\ \ \ \ \ 
      Chop[Fourierfitters[L, n, t] . coeffs[n, Fvalues[F, L, n]]]; \)\ \ \ \ 
    \)

Look at this periodic function f[t]:

\(Clear[f, t]; 
f[t_] = \ 14.8\ \ \((0.35\ t - Round[0.35\ t])\)\^2\ \ ; 
    
thigh\  = \ 12; 

    fplot = Plot[f[t], {t, 0, thigh}, 
        PlotStyle -> {{Thickness[0.02], Blue}}, PlotRange -> All, 
        AxesLabel -> {"t", "f[t]"}, AspectRatio -> 1\/2]; 
    \)

Now look at this  fast Fourier fit of f[t]:

L = (1/0.35); n = 7; \(Clear[approxf]; 

    approxf[t_] = Chop[ComplexExpand[FastFourierfit[f, L, n, t]]]\)fitplot = 
      Plot[{f[t], approxf[t]}, {t, 0, 20}, AspectRatio -> (1/2), 
        PlotRange -> All, 
        PlotStyle 
          -> {{Thickness[0.015], Blue}, {Thickness[0.01], 
              CadmiumOrange}}];

Use what you see to come up with values of c so that the undamped forced oscillator 
           y''[t]+c y[t]=f[t]
are fried by going into resonance or near resonance (beating) .

Take another look at that good Fourier fit:

approxf[t]

That big coefficient on the Cos[2.19911 t] term tells you that when you put c = (2.19911^2), you can count on  the undamped forced oscillator 
           y''[t]+c y[t]=f[t]
to go into resonance or near resonance (beating) .
See it happen:

Clear[f, t, approxY, approxy, s];
b =  0;
c  =  (2.19911^2);
endtime  =  120;
f[t_]  =  approxf[t];
ystarter  =  Random[Real, {-1, 1}];
yprimestarter  =  Random[Real, {-1, 1}];
approxY[s_]  =  
        (1/((s^2) + b s + c)) 
          (LaplaceTransform[f[t], t, s] +  b ystarter +  s ystarter  + 
               yprimestarter);
approxy[t_]  =  
        Chop[Expand[InverseLaplaceTransform[approxY[s], s, t]]];

\(approxplot\  = 
        Plot[approxy[t], {t, 0, endtime}, 
          PlotStyle -> {{Thickness[0.01], Red}}, PlotPoints -> 2\ endtime, 
          AspectRatio -> 1/GoldenRatio, \ \ 
		\ \ \ \ \ \ \ \ \ 
          PlotRange -> All, AxesLabel -> {"t", "y[t]"}, 
          PlotLabel -> c\ ]\);

Fried.
In fact, if you take any c near (2.19911^2),
then you can still count on resonance or near resonance (beating):

Clear[f, t, approxY, approxy, s];
b =  0;
c  =  (2.19911^2)  + Random[Real, {-0.3, 0.3}];
endtime  = 120;
f[t_]  =  approxf[t];
ystarter  =  Random[Real, {-1, 1}];
yprimestarter  =  Random[Real, {-1, 1}];
approxY[s_]  =  
        (1/((s^2) + b s + c)) 
          (LaplaceTransform[f[t], t, s] +  b ystarter +  s ystarter  + 
               yprimestarter);
approxy[t_]  =  
        Chop[Expand[InverseLaplaceTransform[approxY[s], s, t]]];

\(approxplot\  = 
        Plot[approxy[t], {t, 0, endtime}, 
          PlotStyle -> {{Thickness[0.01], Red}}, PlotPoints -> 2\ endtime, 

          		AspectRatio -> 1/GoldenRatio, \ \ 
		\ \ \ \ \ \ \ \ \ 
          PlotRange -> All, AxesLabel -> {"t", "y[t]"}, 
          PlotLabel -> c]\);

To find even more resonating c's, look again at the good Fourier fit:

approxf[t]

That fairly big coefficient on the Cos[4.39823 t] term tells you that when you put c = (4.39823^2), you can count on  the undamped forced oscillator 
           y''[t]+c y[t]=f[t]
to go into resonance or near resonance (beating) .
See it happen:

\(Clear[f, t, approxY, approxy, s]; 
b = \ 0; 
c\  = \ 4.39823\^2; 

    endtime\  = \ 80; 
f[t_]\  = \ approxf[t]; 

    ystarter\  = \ Random[Real, {\(-1\), 1}]; 

    yprimestarter\  = \ Random[Real, {\(-1\), 1}]; 

    approxY[s_]\  = \ 
      \(1\/\(s\^2 + b\ s + c\)\) 
        \((LaplaceTransform[f[t], t, s] + \ b\ ystarter + \ s\ ystarter\  + \ 
            yprimestarter)\); 

    approxy[t_]\  = \ 
      Chop[Expand[InverseLaplaceTransform[approxY[s], s, t]]]; 


    approxplot\  = 
      Plot[approxy[t], {t, 0, endtime}, 
        PlotStyle -> {{Thickness[0.01], Red}}, PlotPoints -> 2\ endtime, 
        AspectRatio -> 1/GoldenRatio, \ \ 
		\ \ \ \ \ \ \ \ \ 
        PlotRange -> All, AxesLabel -> {"t", "y[t]"}, 
        PlotLabel -> c\ ]; \)

To find even more resonating c's, look again at the good Fourier fit:

approxf[t]

That fairly big coefficient on the Cos[6.59734457253856465` t] term tells you that when you put \(c = 6.59734457253856465`\^2\), you can count on  the undamped forced oscillator 
           y''[t]+c y[t]=f[t]
to go into resonance or near resonance (beating) .
See it happen:

Clear[f, t, approxY, approxy, s];
b =  0;
c  =  (6.59734^2);
endtime  =  80;
f[t_]  =  approxf[t];
ystarter  =  Random[Real, {-1, 1}];
yprimestarter  =  Random[Real, {-1, 1}];
approxY[s_]  =  
        (1/((s^2) + b s + c))(
          LaplaceTransform[f[t], t, s] +  b ystarter +  s ystarter  +  
            yprimestarter);
approxy[t_]  =  
        Chop[Expand[InverseLaplaceTransform[approxY[s], s, t]]];

\(approxplot\  = 
        Plot[approxy[t], {t, 0, endtime}, 
          PlotStyle -> {{Thickness[0.01], Red}}, PlotPoints -> 2\ endtime, 
          AspectRatio -> 1/GoldenRatio, \ \ 
		\ \ \ \ \ \ \ \ \ 
          PlotRange -> All, AxesLabel -> {"t", "y[t]"}, 
          PlotLabel -> c\ ]\);

\(TextForm\`\(You\ get\ the\ \(idea . \)
		\)\)

T.1.a.ii)

Why did that work?

When you go with y''[t]+(c^2) y[t]=0 ( with c > 0).
 the natural output consists of combinations of 
          Cos[c t] and Sin[c  t].
 When you force y''[t]+(c^2) y[t]=0 with these combinations you get
          y''[t]+(c^2) y[t]= A Cos[c t]  +  B Sin[c t]
which is guaranteed to be in true resonance.
So when you see an approximate forcing function such as
          2.0 Sin[2 t] - 0.4  Cos[4  t] + 0.5 Sin[6 t],
 then you can be sure that each of the following undamped forceds oscillators is in resonance or near resonance:
          y''[t] + 4 y[t]  = 2.0 Sin[2 t] - 0.4 Cos[4 t] + 0.5 Sin[6 t],
          y''[t] + 16 y[t]  = 2.0 Sin[2 t] - 0.4 Cos[4 t] + 0.5 Sin[6 t], and
          y''[t] + 36 y[t]  = 2.0 Sin[2 t] - 0.4 Cos[4 t] + 0.5 Sin[6 t].

T.2)  The Laplace transform of yzeroinput[t] is the same as the product of Laplace transform of f[t] and the Laplace transform of yunitimpulse[t]

T.2.a.i)

Given a forced oscillator diffeq
		y''[t]+b y'[t]+c y[t]=f[t],
you know  yunitimpulset[t]  solves
             y''[t]+b y'[t]+c y[t]=0
              with y[0] = 0 and y'[0]=1.
The Laplace transform,Yunitimpulset[t]  of yunitimpulset[t]  is

Clear[f, t, b, c, Yunitimpulse , y, ystarter, yprimestarter, s];
f[t_]  =  0;
ystarter  =  0;

yprimestarter =  1;

Yunitimpulse[s_] = (1/((s^2) + b s + c))(
          LaplaceTransform[f[t], t, s] +  b ystarter +  s ystarter  +  
            yprimestarter)

Short and sweet.

Given a forced oscillator diffeq
		y''[t]+b y'[t]+c y[t]=f[t],
you know yzeroinput[t]  solves
             y''[t]+b y'[t]+c y[t]=0
              with y[0] = 0 and y'[0]=1.
The Laplace transform, Yzeroinput[t]  of yzeroinput[t]  is

\(Clear[f, t, b, c, Yzerinput, y, ystarter, yprimestarter, s]; 

    ystarter\  = \ 0; 
yprimestarter = \ 0; 


    Yzeroinput[s_]\  = \ 
      \(1\/\(s\^2 + b\ s + c\)\((
        LaplaceTransform[f[t], t, s] + \ b\ ystarter + \ s\ ystarter\  + \ 
          yprimestarter)\)\)\)

Two benefits:
-> This allows you to write down the Laplace transform of yzeroinput[t] directly in terms of the Laplace transform of f[t] and the coefficients b and c with no extra calculation. 
 
-> And this reveals that the Laplace transform of yzeroinput[t] is the same as the product of Laplace transform of f[t] and the Laplace transform of yunitimpulse[t]:

LaplaceTransform[f[t], t, s]  Yunitimpulse[s]

Is the fact that the Laplace transform of yzeroinput[t] is the same as the product of Laplace transform of f[t] and the Laplace transform of yunitimpulse[t] an accident?

Get ready.

T.2.a.ii)

Look at this calculation of the convolution of two cleared functions f[t] and g[t]:

\(Clear[f, g, h, t, s, x, y]; 

    h[t_] = ∫\_0\%t\( g[t - x]\ f[x]\) \[DifferentialD]x\)

Calculate the Laplace transform of h[t]:

LaplaceTransform[h[t], t, s]

This tells you that the Laplace transform of the convolution
             h[t] = (∫_0^t) f[t - x] g[x] dx
is the Laplace transform of f[t] times the Laplace transform of g[t].

Use this fact to explain why the Laplace transform of  yzeroinput[t] is the same as the product of Laplace transform of f[t] and the Laplace transform of yunitimpulse[t].

Remember that 
       \(TextForm
      \`yzeroinput[t]\  = 
        ∫\_0\%t yunitimpulse[t - x]\ f[x] 
            \(dx . \)\)
This is the same as 
       \(TextForm
      \`h[t] = ∫\_0\%t\ f[x]\ g[t - x] dx\)
with g[t - x] = yunitimpulse[t - x]  and h[t] = yzeroinput[t] .

According to what was said above,the Laplace transform of  
       yzeroinput[t] = h[t] 
is the product of the  Laplace transform of yunitimpulse[t] = g[t]  and the Laplace transform of f[t].

T.3) A good way of screwing up your Fourier fits

T.3.a)

Activate this code.

\(Clear[FastFourierfit, Fourierfitters, F, Fvalues, n, k, 
\ \ \ \ \ \ 
      jump, num, numtab, coeffs, t, L]; 

    jump[n_]\  := \ \(jump[n]\  = \ N[1\/\(2\ n\)]\); 

    Fvalues[F_, L_, n_]\  := \ 
\ \ 
      N[Table[F[L\ t], {t, 0, 1\  - \ jump[n], jump[n]}]]; 
\ \ 


    numtab[n_]\  := \ \(numtab[n]\  = \ Table[k, {k, 1, n}]\); 


    Fourierfitters[L_, n_, t_]\  := \ 
      Table[E\^\(\(2\ π\ I\ k\ t\)\/L\), 

        \ \ \ \ \ \ \ \ \ \ {k, \(-n\)\  + \ 1, n\  - \ 1}]; 

    coeffs[n_, list_]\  := \ 
      Join[Reverse[Part[Fourier[list], numtab[n]]], 
	
          Part[InverseFourier[list], Drop[numtab[n], 1]]]/
		
        N[Sqrt[Length[list]]]
\)\(\(FastFourierfit[F_, L_, n_, t_]\  := \ 
\ \ \ \ \ 
      Chop[Fourierfitters[L, n, t] . coeffs[n, Fvalues[F, L, n]]]; \)\ \ \ \ 
    \)

Look at this plot of a periodic function:

\(Clear[f, t]; 
f[t_] = \ Cos[1.3  t]\ Sign[Sin[1.3\ \ t]]\  + \ 1\ \ ; 
    
thigh\  = \ 12; 


    fplot = Plot[f[t], {t, 0, thigh}, 
        PlotStyle -> {{Thickness[0.02], Blue}}, PlotRange -> All, 
        AxesLabel -> {"t", "f[t]"}, AspectRatio -> 1\/2]; 
    \)

Now look at this attempt at a fast Fourier fit of f[t]:

\(L = π; 
n = 15; 
Clear[approxf]; 


    approxf[t_] = Chop[ComplexExpand[FastFourierfit[f, L, n, t]]]\)fitplot = 
      Plot[{f[t], approxf[t]}, {t, 0, thigh}, AspectRatio -> (1/2), 
        PlotRange -> All, 
        PlotStyle 
          -> {{Thickness[0.015], Blue}, {Thickness[0.01], 
              CadmiumOrange}}];

What went wrong?
How do you fix it?

Take another look:

Show[fitplot];

This plot resulted from feeding L = π  into the FourierFit instruction.
The trouble is that L = π  is not a period of f[t]. 
The formula for f[t] is:

f[t]

Calculate the period of f[t]:

Solve[1.3 t == 2 π, t]

Set L = 4.83322 and run everything again:

\(L = 4.83322; 
n = 10; 
Clear[approxf]; 

    approxf[t_] = Chop[ComplexExpand[FastFourierfit[f, L, n, t]]]\)fitplot = 
      Plot[{f[t], approxf[t]}, {t, 0, thigh}, AspectRatio -> (1/2), 
        PlotRange -> All, 
        PlotStyle 
          -> {{Thickness[0.015], Blue}, {Thickness[0.01], 
              CadmiumOrange}}];

Lookin' real nice.
To get a good Fourier fit, you've got to set the period L of the Fourier fit to be equal to the period of the function you are fitting.

T.4) Inserting an extra t when the characteristic equation doesn't give two different solutions

T.4.a.i)

Here's an oscillator differential equation with starter data:

b = 6; c = 9; Clear[t, y];
oscdiffeq=y''[t]+b y'[t]+c y[t]==0

Here's what happens when you use the charcteristic equation to go after formulas for the solution of this oscillator differential equation with sample starter data 
      y[0] = 1.2  and  y'[0]=3.4.

\(Clear[z]; 
charequation = z\^2 + b\ z + c == 0; 

    zsols = Solve[charequation, z]\)

\(Clear[combo, C1, C2]; 

    combo[t_] = 
      C1\ E\^\(zsols[[1, 1, 2]]\ t\) + 
        C2\ E\^\(zsols[[2, 1, 2]]\ t
              \)\)

ystarteq = combo[0] == 1.2

yprimestarteq=combo'[0]==3.4

Csols = Solve[{ystarteq, yprimestarteq}, {C1, C2}]

What went wrong?

Look at combo[t]:

combo[t]

At first glance there appear to be two degrees of freedom here, but because the charactertistic equation has only one solution, you have only one degree of freedom. This is not enough to allow you to match up with the starter data on y[0] and y'[0].

T.4.a.i) Inserting a t

What can you do to get around this problem?

Old DiffEq hands know how to tweak combo[t]. 
Here's how they do it:
Look at combo[t]:

combo[t]

Insert a t onto one of the terms:

Clear[bettercombo];
bettercombo[t_]=combo[t]/.C2 (E^(-3 t))->t C2 (E^(-3 t))

Try it out:

\(ExpandAll[

      oscdiffeq /. {y[t] -> bettercombo[t], 
\ \ \ \ \ \ \ \ \ \ \ \ 
          \(y'\)[t] -> \(bettercombo'\)[t], 
\ \ \ \ \ \ \ \ \ \ \ \ 
          \(\(y'\)'\)[t] -> \(\(bettercombo'\)'\)[t]}]\)

The upshot:
bettercombo[t] solves the differential equation and gives you the two genuine degrees of freedom you need to match up with the starter data on y[0] and y'[0].

ystarteq = bettercombo[0] == 1.2; yprimestarteq=bettercombo'[0]==3.4;Csols = Solve[{ystarteq, yprimestarteq}, {C1, C2}]

This gives you the exact formula:

\(Clear[yformula]; 

    yformula[t] = 
      bettercombo[t] /. Csols[[1]]\)

Inserting the extra  t did the trick.

T.4.a.ii)

When do you insert an extra t as in part i) above?

If there are two different solutions of the characteristic equation, then you don't insert any extra t. In fact if you do, you will have a guaranteed screw-up.
If the characteristic equation has only one solution (double root), then you need to insert the t as in part i) above.

T.4.a.iii) How important is the issue of the extra t?

How important is the issue of the extra t?

If you want an exact formula, it's important.
If you are willing to settle for a very good approximate formula then it's not an issue at all.
Watch this:

b = 3; c = ((((b/2)))^2); Clear[t, y];
oscdiffeq=y''[t]+b y'[t]+c y[t]==0

Here's what happens when you use the characteristic equation to go after formulas for the solution of this oscillator differential equation with starter data 
      y[0] = 2  and  y'[0]=1.

\(Clear[z]; 
charequation = z\^2 + b\ z + c == 0; 

    zsols = Solve[charequation, z]\)

Only one solution of the characteristic equation. If you want an exact formula, then you need to insert the t:

Clear[exacty];
exactsol=DSolve[{oscdiffeq,y[0] == 2,y'[0]==1},y[t],t];
exacty[t_] = 
        y[t] /. exactsol[[1]]

There's that extra t. 
If you are willing to settle for a good approximate formula, you look at the original differential equation:

b = 3; c = ((((3/2)))^2); Clear[t, y];
diffeq=y''[t]+b y'[t]+c y[t]==0

Tweak the damping term b a little bit to get a nearby oscillator:

tweakedb = 1.0001 b; Clear[t, y];
nearbydiffeq=y''[t]+tweakedb y'[t]+c y[t]==0

Solve the nearby oscillator diffeq via the characteristic equation:

\(Clear[z]; 
charequation = z\^2 + tweakedb\ z + c == 0; 

    zsols = Solve[charequation, z]\)

\(Clear[combo, C1, C2]; 

    combo[t_] = 
      C1\ E\^\(zsols[[1, 1, 2]]\ t\) + 
        C2\ E\^\(zsols[[2, 1, 2]]\ t
              \)\)

ystarteq = combo[0] == 2

yprimestarteq=combo'[0]==1

Csols = Solve[{ystarteq, yprimestarteq}, {C1, C2}]

\(Clear[approxy]; 

    approxy[t_] = 
      combo[t] /. Csols[[1]]\)

No extra t.
Compare the two solutions:

endtime = 12; Plot[{exacty[t], approxy[t]}, {t, 0, endtime}, 
      PlotStyle -> {{Blue, Thickness[0.02]}, {Red, Thickness[0.01]}}, 
      PlotRange -> All, AxesLabel -> {"t", ""}, 
      AspectRatio -> (1/2)];

Sharing ink all the way.
You can't tell the difference without a scorecard.

The upshot: The matter of inserting the t is important only when theoretical exactness is paramount.
In practical situations, you can almost always tweak the damping term to get good, practical results.

T.4.b.i)

Use the Laplace transform to explain where the extra t comes from.

Here's another oscillator differential equation which needs the extra t:

\(\(\ b\  = \ 8; 
\ c\  = \ 16; 
Clear[t, y]; 

    oscdiffeq\  = \ 
      \((\(\(y'\)'\)[t]\  + \ b\ \(y'\)[t]\  + \ c\ y[t]\  == \ 0)\)\)\)

Here's what happens when you use the characteristic equation to go after formulas for the solution of this oscillator differential equation with starter data 
      y[0] = 2  and  y'[0]=3.

\(Clear[z]; 
charequation = z\^2 + b\ z + c == 0; 

    zsols = Solve[charequation, z]\)

Uh-oh.
Only one solution of the characteristic equation. To see where the extra t comes from, go to the Laplace Transform:

\(Clear[t, Y, approxy, s]; 
f[t_] = 0; ystarter = 2; yprimestarter = 3; 
    Y[s_] = \(LaplaceTransform[f[t], t, s] + b\ ystarter + s\ ystarter + 
          yprimestarter\)\/\(s\^2 + b\ s + c\); 
    y[t_] = Chop[Expand[InverseLaplaceTransform[Y[s], s, t]]]\)

Now you can see where the extra t comes from. 
Look at the Laplace transform of the solution y[t]:

Apart[Y[s]]

The inverse Laplace transform of 
          (2/(s + 4))   is   2 (E^((-4 t)))
and the inverse Laplace transform of 
           (11/(((4 + s))^2)) = -11 D[((1/(s + 4))), s]  is   11 t (E^((-4 t))):

2 InverseLaplaceTransform[(1/(s + 4)), s, t]

11 InverseLaplaceTransform[-D[(1/(s + 4)), s], s, t]

That's where the t comes from.

T.4.b.ii)

In the last part, knowing that the inverse Laplace transform of
          (2/(s + 4))   is   2 (E^((-4 t))),
 you recognized that the inverse Laplace transform of 
         (11/(((4 + s))^2)) = -11 D[(1/(s + 4)), s]   is  - t (E^((-4 t))).
(That's where the extra t really comes from.)
Is this a special case of a general rule?

You bet your sweet pocket calculator!
The rule is this:
If 
         Y[s] is the Laplace transform of y[t], 
then the inverse Laplace transform of 
         Y'[s] is -t y[t].
It's even programmed into Mathematica:

Clear[s, t, y, Y];
Y[s_] = LaplaceTransform[y[t], t, s];
InverseLaplaceTransform[Y'[s],s,t]

You can explain this sweet little rule as follows:
       Y[s] = 
        (∫_0^∞) 
          (E^(-s t)) y[t] dt
       Y'[s] = 
        (∂_s)(
              (∫_0^∞) 
                (E^(-s t)) y[t] dt) = 
          (∫_0^∞) 
            (∂_s)((E^(-s t)) y[t]) 
              dt
       Y'[s] = 
        (∫_0^∞) 
          (E^(-s t)) y[t] (-t) dt
       Y'[s] = 
        (∫_0^∞) 
          (E^(-s t)) (- t y[t]) dt

This tells you that Y'[s] is the Laplace transform of  -t y[t].
Nice, simple calculus.

T.4.b.iii)

Is there a way you can see where the extra t comes from without resorting to the black box of the Laplace transform?

Yes there is!
The prototype oscillator differential equation whose characteristic equation has only one solution is:

\(Clear[oscdiffeq, r, y, t]; 
b = 2\ r; 
\(c = r\^2; \)\)oscdiffeq=y''[t]+b y'[t]+c y[t]==0

Check:

\(Clear[z]; 
charequation = z\^2 + b\ z + c == 0; 

    zsols = Solve[charequation, z]\)

Good. 
Now tweak this oscillator this way:

Clear[nearbyoscdiffeq, y, h, t];
nearbyoscdiffeq[h_]=y''[t]+(2 r + h) y'[t]+r (r + h) y[t]==0

When h is very small, then the two differential equations are almost the same.

ColumnForm[{oscdiffeq, nearbyoscdiffeq[0.00001]}]

Now look at formulas for solutions of both:

Clear[yexact, p, q];
exactsol=DSolve[{oscdiffeq,y[0] == p,y'[0]==q},y[t],t];
yexact[t_] = 
        y[t] /. exactsol[[1]]

\(Clear[ynearby, p, q]; 

    nearbysol\  = \ 
\ \ \ 
      DSolve[{nearbyoscdiffeq[h], y[0]\  == \ p, \(y'\)[0]\  == \ q}, y[t], 
        t]; 
ynearby[t_, h_]\  = \ y[t] /. nearbysol[\([1]\)]\)

Put:

\(Clear[g, x]; 
g[x_] = \(-E\^\(\(-x\)\ t\)\)\)

And notice that ynearby[t, h] is given by:

\(Clear[altynearby]; 

    altynearby[t_, h_] = 
      \(\((q + p\ r)\)\ \((g[r + h] - g[r])\)\)\/h + p\/E\^\(r\ t\)\)

Check:

Simplify[ynearby[t, h] - altynearby[t, h]]

Good.
As h closes in on 0,
      nearbyoscdiffeq[h] closes in on oscdiffeq;
so 
      ynearby[t, h] closes in on yexact[t].
But 
      ynearby[t, h]=(((q + p r) (g[r + h] - g[r]))/h) + 
            (p/(E^((r t)))), 
and as h closes in on 0, this closes in on:

(q + p r) g'[r]+(p/(E^(r t)))

This is the same as:

yexact[t]

So the extra t comes from taking a derivative.

T.5)  The Gibbs phenomenon

Activate this code.

\(Clear[FastFourierfit, Fourierfitters, F, Fvalues, n, k, 
\ \ \ \ \ \ 
      jump, num, numtab, coeffs, t, L]; 

    jump[n_]\  := \ \(jump[n]\  = \ N[1\/\(2\ n\)]\); 

    Fvalues[F_, L_, n_]\  := \ 
\ \ 
      N[Table[F[L\ t], {t, 0, 1\  - \ jump[n], jump[n]}]]; 
\ \ 


    numtab[n_]\  := \ \(numtab[n]\  = \ Table[k, {k, 1, n}]\); 


    Fourierfitters[L_, n_, t_]\  := \ 
      Table[E\^\(\(2\ π\ I\ k\ t\)\/L\), 

        \ \ \ \ \ \ \ \ \ \ {k, \(-n\)\  + \ 1, n\  - \ 1}]; 

    coeffs[n_, list_]\  := \ 
      Join[Reverse[Part[Fourier[list], numtab[n]]], 
	
          Part[InverseFourier[list], Drop[numtab[n], 1]]]/
		
        N[Sqrt[Length[list]]]
\)\(\(FastFourierfit[F_, L_, n_, t_]\  := \ 
\ \ \ \ \ 
      Chop[Fourierfitters[L, n, t] . coeffs[n, Fvalues[F, L, n]]]; \)

    \ \ \ \ \)

Look at this plot:

\(Clear[f, t]; 
f[t_] = Sign[Sin[π\ t]]; 
\(L = 2; \)\)\(cycles = 3; 


    fplot = Plot[f[t], {t, 0, cycles\ L}, 
        PlotStyle -> {{Thickness[0.02], Blue}}, 
        AxesLabel -> {"t", "f[t]"}, 
        PlotLabel -> "cycles"\ cycles, 
        Epilog -> {{Red, Thickness[0.02], Line[{{0, 0}, {L, 0}}]}, {
              Text["One Period", {L\/2, 0.1}]}}]; \)

Now look at this fast Fourier fit of f[t] shown with f[t]:

n = 8; \(Clear[realfit]; 

    realfit[t_] = Chop[ComplexExpand[FastFourierfit[f, L, n, t]]]; 

    \(fitplot = 
      Plot[{f[t], realfit[t]}, {t, 0, L}, 
        PlotStyle -> {{Thickness[0.01], Red}, {Blue}}]; \)\)

T.5.a)

Late last century, the American mathematician Gibbs noticed theoretically (with no machine help) the same phenomenom you are seeing here. Lots of folks call this by the name "Gibbs phenomenon."
What is the Gibbs phenomenon?

Look again:

Show[fitplot];

Increase the quality of the Fourier fit:

\(\(n\  = \ 12; 
Clear[realfit]; 

    realfit[t_]\  = \ 
Chop[ComplexExpand[FastFourierfit[f, L, n, t]]]; 


    newfitplot\  = \ 
      Plot[{f[t], realfit[t]}, {t, 0, L}, 
\ \ \ \ \ 
        PlotStyle -> {{Thickness[0.01], Red}, {Blue}}]; 
\ \ \ \ \ \)\)

Look at the jump f[t] takes in the middle.  Notice that the fast Fourier fit overshoots on one side of the jump with a compensating undershoot on the the other side of the jump.  
Gibbs proved theoretically that when you go with a high quality fast Fourier fit (big n) of any periodic function f[t], then you will always get this phenomenon at a jump of f[t].