DE.03 Laplace Transform and Fourier Analysis
 BASICS

\(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_] := 
      Join[Reverse[
            Fourier[list][[numtab[n]
                ]]], 
          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]]]; \(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]]

B.1) The Laplace transform Y[s] of y[t] is 
                          Y[s] = 
        (∫_0^(infinity ))(E^(-s t)) y[t] 
            dt. 
The Laplace transform  of y'[t]  is   s Y[s] - y[0].
The Laplace transform  of y''[t]  is   (s^2) Y[s]-s y[0]-y'[0]

Here's the inside scoop on the Laplace transform:
The functions legal for the Laplace transform are all functions f[t] with the property that
        (E^(-s t)) f[t] -> 0 as t -> ∞
for all large positive s.
      Sin[p t], Cos[p t] , (E^(k t)),  Log[t]
and any quotient of polynomials are all legal.
The Laplace transform of a given function f[t] is another function, F[s], given by
      F[s] = 
        (∫_0^(infinity ))(E^(-s t)) f[t] 
            dt.

For instance, here's the Laplace transform of
       f[t] = t:

\(Clear[f, F, s, t]; 
f[t_] = t; 

    F[s_] = LaplaceTransform[f[t], t, s]\)

Check:

(∫_0^∞) (E^(-s t)) f[t] 
        \[DifferentialD]t

This is Mathematica's way of telling you that if s > 0, then 
     (∫_0^infinity) (E^(-s t)) f[t] dt = 
        (1/(s^2)).
If Mathematica is given the Laplace transform of f[t], then Mathematica can often recover the formula for f[t]:

InverseLaplaceTransform[F[s], s, t]

Some more play:

Experiment 1:

\(Clear[f, F, s, t]; 
f[t_] = Sin[t]; 

    F[s_] = LaplaceTransform[f[t], t, s]\)

InverseLaplaceTransform[F[s], s, t]

InverseLaplaceTransform[F[s], s, t] == f[t]

Experiment 2:

\(Clear[f, F, s, t]; 
f[t_] = E\^\(-t\); 

    F[s_] = LaplaceTransform[f[t], t, s]\)

InverseLaplaceTransform[F[s], s, t]

InverseLaplaceTransform[F[s], s, t] == f[t]

Experiment 3:

\(Clear[f, F, s, t]; 
F[s_] = 2\/\(s\^2 + 4\)\)

f[t_] = InverseLaplaceTransform[F[s], s, t]

LaplaceTransform[f[t], t, s]

Experiment 4:

\(Clear[f, Lf, s, x]; 
F[s_] = 5.2\/\(s + 2\)\)

f[t_] = InverseLaplaceTransform[F[s], s, t]

LaplaceTransform[f[t], t, s]

Different functions have different Laplace transforms.

B.1.a.i) The Laplace transform of the first derivative

Look at this:

Clear[y, t, s];
LaplaceTransform[y'[t],t,s]

One reason folks fall for the charms of the Laplace transform is that if 
      g[t] = y'[t], 
then you are guaranteed that the Laplace transforms line up like this:
      G[s] = s Y[s] - y[0]
where 
      G[s] is the Laplace transform of g[t]=y'[t]
and 
      Y[s] is the Laplace transform of y[t].

This is a good fringe benefit because it tells you how to write down the Laplace transform of g[t]=y'[t] once you've got your hands on y[t] and its Laplace transform. 
Play with it:

Experiment 1:

Clear[y, Y, s, t];
y[t_] = (E^(-0.4 t));
LaplaceTransform[y'[t],t,s]

Together[s LaplaceTransform[y[t], t, s] - y[0]]

Experiment 2:

Clear[y, s, t];
y[t_] = (E^(-t)) Cos[π t];
Together[LaplaceTransform[y'[t],t,s]]

Together[s LaplaceTransform[y[t], t, s] - y[0]]

Even when Mathematica can't handle the calculations, you can be certain that if 
        g[t]=y'[t], 
then you are guaranteed that
        G[s] = s Y[s] - y[0]
because you can use integration by parts to explain this formula.  Do it.

Remember,
        Y[s] = 
        (∫_0^∞) 
          (E^(-s t)) y[t] dt, 
and because g[t]=y'[t]
        G[s] = 
        (∫_(0 )^(∞ )) 
          (E^(-s t)) y'[t] dt, .

The upshot:
Explaining why 
        G[s] = s Y[s] - y[0]
is the same as explaining why
       (∫_0^∞)(E^(-s t)) y'[t]dt=s (∫_0^∞) 
              (E^(-s t)) y[t] dt - y[0].
To get a grasp on the left side, remember
       (∫_0^∞)(E^(-s t)) y'[t]dt  =  limit of (∫_0^(n ))(E^(-s t)) y'[t]dt as n -> ∞

Now work on
       (∫_0^(n ))(E^(-s t)) y'[t]dt 
using integration by parts:
The integration by parts formula is
       (∫_0^n)u[t] v'[t]dt   =    u[t] v[t] \(|\_0\%n\) - (∫_0^n)v[t] u'[t]dt.
Make the assignments:
        u[t] = (E^(-s t)) and  v'[t]=y'[t].
This gives 
        u'[t]=-s (E^(-s t)) and v[t] = y[t]
Plug into the integration by parts formula to get
       (∫_0^(n ))(E^(-s t)) y'[t]dt  
    =  (E^(-s t)) y[t]\(|\_0\%n\)  -  (∫_0^n)-s (E^(-s t)) y[t] 
          dt
    =  (E^(-s t)) y[t]\(|\_0\%n\)  +  
      s (∫_0^n) 
          (E^(-s t)) y[t] dt
    =  (E^(-s n)) y[n]- y[0] +  
      s (∫_0^n) 
          (E^(-s t)) y[t] dt.
As n -> ∞,
       (∫_0^(n ))(E^(-s t)) y'[t]dt   ->   G[s],
        (E^(-s n)) y[n]  ->   0,
and 
        s 
        (∫_0^n) (E^(-s t)) y[t] dt  ->   s Y[s].
So
        G[s] = -y[0] + s Y[s]  =   s Y[s] - y[0]
as advertised.

If someone ever asks you whether you have seen the Laplace transform, now you can say that you've seen a little bit about it.  And you can tell them that you saw it first in DiffEq&Mathematica.

B.1.a.ii)

Look at this:

Clear[y, t, s];
LaplaceTransform[y''[t],t,s]

You already know the Laplace transform of y'[t] is
        s Y[s] - y[0].
How does this fact guarantee that the Laplace transform of y''[t] is
        (s^2) Y[s]-s y[0]-y'[0]?

Put 
       h[t]=y''[t] and g[t]=y'[t]
so that 
       h[t]=g'[t].

According to part a.i) above, you are guaranteed that
        H[s] = s G[s] - g[0]:

Clear[G, H, g, y, Derivative, s];
H[s_]=s G[s]-y'[0]

The formula in part a.i) also says that
        G[s] = s Y[s] - y[0]:

Expand[H[s] /. G[s] -> s Y[s] - y[0]]

So
         H[s]=(s^2) Y[s]-s y[0]-y'[0],
as advertised.

Try this formula out:

Experiment 1:

Clear[y, s, x, t];
y[t_] = Sin[5 t];
LaplaceTransform[y''[t],t,s]

Together[(s^2) LaplaceTransform[y[t], t, s]-s y[0]-y'[0]]

Experiment 2:

Clear[y, s, t];
y[t_] = Sin[t] + (1/2) Sin[2 t];
Together[LaplaceTransform[y''[t],t,s]]

Together[(s^2) LaplaceTransform[y[t], t, s]-s y[0]-y'[0]]

It works in theory even in cases Mathematica can't handle.

B.2) Using the magic of the Laplace Transform to  come up with formulas for forced exponential diffeqs and forced oscillator diffeqs

B.2.a) Using the Laplace transform to solve a forced exponential diffeq

Showcase the Laplace transform by using it to come up with the exact formula for the solution of
     y'[t]+1.2 y[t]=4.75 (E^(-0.09 t)) Cos[3.5  t]
     with  y[0] = 5.2.

With the Laplace transform, this is duck soup. You look at the diffeq:
         y'[t]+1.2 y[t]=4.75 (E^(-0.09 t)) Cos[3.5  t]
with  y[0] = 5.2.
and then you enter

\(Clear[f, t, Y, y, s]; 
r = 1.2; 

    f[t_] = 4.75\ E\^\(\(-0.09\)\ t\)\ Cos[3.5\ t]; 
starter = 5.2; 
    Y[s_] = \(LaplaceTransform[f[t], t, s] + starter\)\/\(r + s\)\)

This is the Laplace transform of the solution. 
To get the formula for the solution y[t] you calculate the inverse Laplace transform of Y[s].
The formula for the solution y[t] is:

y[t_] = Chop[Expand[InverseLaplaceTransform[Y[s], s, t]]]

That's all there is to it.
Check whether this y[t] solves   y'[t]+1.2 y[t]=4.75 (E^(-0.09 t)) Cos[3.5  t]:

Chop[Expand[y'[t]+1.2 y[t]]]

Yepper.
Check whether this y[0] = 5.2 for this y[t]:

y[0]

Everything checks.

B.2.b)  Using the Laplace transform to solve a forced oscillator diffeq

Showcase the Laplace transform by using it to come up with the exact formula for the forced, damped oscillator coming from 
       y''[t]+1.3 y'[t]+6.5 y[t]=9.1 (E^(-0.8 t))
     with  y[0] = 1.6  and  y'[0]=3.7.

With the Laplace transform, this is also duck soup. 
You look at the diffeq:
         y''[t]+1.3 y'[t]+6.5 y[t]=9.1 (E^(-0.8 t))
     with  y[0] = 1.6  and  y'[0]=3.7 . 
and then you enter

\(Clear[f, t, Y, y, s]; 
b = 1.3; c = 6.5; 

    f[t_] = 9.1\ E\^\(\(-0.8\)\ t\); 
ystarter = 5.2; yprimestarter = 3.7; 
    Y[s_] = \(LaplaceTransform[f[t], t, s] + b\ ystarter + s\ ystarter + 
          yprimestarter\)\/\(s\^2 + b\ s + c\)\)

This is the Laplace transform of the solution. 
To get the formula for the solution y[t] you calculate the inverse Laplace transform of Y[s].
The formula for the solution y[t] is:

y[t_] = Chop[Expand[InverseLaplaceTransform[Y[s], s, t]]]

That's all there is to it.
Check whether this y[t] solves y''[t]+1.3 y'[t]+6.5 y[t]=9.1 (E^(-0.8 t)):

Chop[Expand[((y'))'[t]+1.3 y'[t]+6.5 y[t]]]

Yepper.
Check whether  y[0] = 1.6  and  y'[0]=3.7 .

y[0]

y'[0]

Everything checks. 
Amazing thing this Laplace transform.

B.2.c)  Using the Laplace transform to solve a forced oscillator diffeq

You are using Mathematica to handle the Laplace transform.
How was it done in your grandparents' diffeq class?

In your grandparents' time, they would look at a diffeq such as 

         y''[t]+2.0y'[t]+6.5 y[t]=9.1 (E^(-1.9 t))
     with  y[0] = 1.5  and  y'[0]=-4.7 . 
and write down
Y[s]  =  
        ((LaplaceTransform[f[t], t, s] +  b ystarter +  s ystarter  +  
            yprimestarter)/((s^2) + b s + c))
  with b = 2.0, c = 6.5,  f[t]=9.1 (E^(-1.9 t)), ystarter = 1.5 and yprimestarter = -4.7.
Then your grandparents would look in a printed table to find which function y[t] has Y[s] for its Laplace transform.
  
This is similar to what you are doing, but you're using the computer to find which function y[t] has Y[s] for its Laplace transform.

B.2.d)

Part of the job of mathematics is to explain why things are guaranteed to work out. 
Explain why the methods used in parts a) and b) above worked.

Explanation for the forced exponential diffeq 
         y'[t] + r y[t] = f[t]:

It's just a little algebra.
Enter the differential equation with everything cleared:

\(Clear[y, t, f, Y, r, s]; 

    diffeq\  = \ \(y'\)[t]\  + \ r\ y[t]\  == \ f[t]\)

Replace all the functions involved by their Laplace transforms:

\(laplaced\  = \ 
      diffeq /. 

        \ \ \ \ \ \ \ \ \ {y[t] -> \ LaplaceTransform[y[t], t, s], 

          \ \ \ \ \ \ \ \ \ \ 
          \(y'\)[t] -> \ LaplaceTransform[\(y'\)[t], t, s], 

          \ \ \ \ \ \ \ \ \ \ f[t] -> LaplaceTransform[f[t], t, s]}\)

Solve this for the Laplace Transform of y[t]:

sol = Simplify[Solve[laplaced, LaplaceTransform[y[t], t, s]]]

There is is.
This tells you that the Laplace transform Y[s] of the solution y[t] is
                      Y[s] = ((LaplaceTransform[f[t], t, s] +  y[0])/(r + s))
which is exactly the formula used in part a).

Explanation for the forced oscillator diffeq 
         y''[t] + b y'[t]  + c y[t] = f[t]:

Again, it's just a little algebra.
Enter the differential equation with everything cleared:

\(Clear[y, t, f, b, c, s, Y]; 

    diffeq\  = \ 
      \(\(y'\)'\)[t]\  + \ b\ \(y'\)[t]\  + \ c\ y[t]\  == \ f[t]\)

Replace all the functions involved by their Laplace transforms:

\(laplaced\  = \ 
      diffeq /. 

        \ \ \ \ \ \ \ \ \ {y[t] -> \ LaplaceTransform[y[t], t, s], 

          \ \ \ \ \ \ \ \ \ \ 
          \(y'\)[t] -> \ LaplaceTransform[\(y'\)[t], t, s], 

          \ \ \ \ \ \ \ \ \ \ 
          \(\(y'\)'\)[t] -> \ LaplaceTransform[\(\(y'\)'\)[t], t, s], 

          \ \ \ \ \ \ \ \ \ \ f[t] -> LaplaceTransform[f[t], t, s]}\)

Solve this for the Laplace Transform of y[t]:

sol = Solve[laplaced, LaplaceTransform[y[t], t, s]]

There is is.
This tells you that the Laplace transform Y[s] of the solution y[t] is
                      Y[s]=((LaplaceTransform[f[t], t, s] + b y[0]+ s y[0]+ y'[0])/(c +  b s + (s^2)))
which is exactly the formula used in part b).

B.3) Fourier fit of periodic functions

B.3.a)

Look at these:

Function 1:

\(Clear[f, t]; 
f[t_] = Cos[2\ π\ t]; 
\(L = 1; \)\)cycles = 2; 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}]}}];

Function 2:

\(Clear[f, t]; 
f[t_] = t\/2 - Floor[t\/2]; 
\(L = 2; \)\)\(cycles = 4; 


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

Function 3:

\(Clear[f, t]; 
f[t_] = \(-0.4\) + E\^\(\(-1.5\)\ \((t - Floor[t])\)\); 
    
\(L = 1; \)\)\(cycles = 4; 


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

Function 4:

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


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

Folks call all of these functions periodic.
How do you recognize periodic functions?

Most folks, like you, just eyeball a plot looking for cycles.  
Once you spot a cycle, then you can say that its length L represents one period of the function.  This gives you 
          f[t + L] = f[t] for all t's.

B.3.b.i) Fast Fourier fit of periodic functions

At the December, 1807 meeting of the French Academy of Sciences, mathematician and engineer Joseph Fourier announced that you can approximate an arbitrary function by sums of sine and cosine waves. Those in attendance thought he was slightly off his rocker, but eventually they were all forced to change their minds because time has confirmed that Fourier's idea was one of the great scientific break-throughs of all time.  In this problem, you get the opportunity to get your hands on Fourier's idea in action.
Go for it.

\(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]]]; \)\ \ \ \ 
    \)

To get an idea of what these new instructions can do for you, look at this plot of a function f[t] whose period is L = 2, plotted for three cycles:

\(Clear[f, t]; 
f[t_] = Cos[1.5\ \((t\/2 - Floor[t\/2])\)]; 
\(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.03], Line[{{0, 0}, {L, 0}}]}, {
              Text["One Period", {(L/2), 0.1}]}}, AspectRatio -> 1];

Now stay with the same function f[t]:

f[t]

Stay with L = 2 and go with n = 4 and look at:

n = 4; FastFourierfit[f, L, n, t]

That ran really fast.
Now look at this plot of f[t] and FastFourierfit[f, L, n, t] with 
          n = 4 and L = 2:

\(Clear[realfit]; 

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

See how realfit[t] plots out:

fitplot = 
      Plot[{f[t], realfit[t]}, {t, 0, cycles L}, 
        PlotStyle -> {{Thickness[0.02], Blue}, {Thickness[0.01], Red}}, 
        AxesLabel -> {"t", ""}, 
        PlotLabel -> "cycles" cycles, 
        AspectRatio -> (1/GoldenRatio), 
        Epilog -> {{Red, Thickness[0.02], Line[{{0, 0}, {L, 0}}]}, {
              Text["One Period", {(L/2), 0.1}]}}];

What's happening here?

Look at this embellished plot of the first cycle still going with L = 2 and n = 4:

first = 
      Plot[{f[t], realfit[t]}, {t, 0, L}, 
        PlotStyle -> {{Thickness[0.01], Blue}, {Thickness[0.01], Red}}, 
        AxesLabel -> {"t", ""}, 
        DisplayFunction -> Identity]; fdata = 
      Table[N[{t, f[t]}], {t, 0, L - (L/(2 n - 1)), (L/(2 n - 1))}]; 
    fdataplot = 
      ListPlot[fdata, PlotStyle -> PointSize[0.04], 
        DisplayFunction -> Identity]; Show[first, fdataplot, DisplayFunction -> $DisplayFunction];

Evidently, Fourierfit[f, n, t] picks 2 n - 1 equally spaced data points off the plot of f[t] between t = 0 and t = L and tries to fit these points with a combination of complex exponentials.

When you kick n up, you get a really good fit:

n = 8; FastFourierfit[f, L, n, t]

\(Clear[realfit]; 

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

cycles = 3; fitplot = 
      Plot[{f[t], realfit[t]}, {t, 0, cycles L}, 
        PlotStyle -> {{Thickness[0.02], Blue}, {Thickness[0.01], Red}}, 
        AxesLabel -> {"t", ""}, 
        PlotLabel -> "cycles" cycles, 
        AspectRatio -> (1/GoldenRatio), 
        Epilog -> {{Red, Thickness[0.02], Line[{{0, 0}, {L, 0}}]}, {
              Text["One Period", {(L/2), 0.1}]}}];

Hot plot.

And you get an extra benefit:
When you put the period L of f[t] into the FastFourierFit[f, L, t, n] instruction, then the periods of the complex exponentials used by the fit match the period of the forcing function f[t], and you continue to get that good fit all the way:

cycles = 8; fitplot = 
      Plot[{f[t], realfit[t]}, {t, 0, cycles L}, 
        PlotStyle -> {{Thickness[0.02], Blue}, {Thickness[0.01], Red}}, 
        AxesLabel -> {"t", ""}, 
        PlotLabel -> "cycles" cycles, 
        AspectRatio -> (1/GoldenRatio), 
        Epilog -> {{Red, Thickness[0.02], Line[{{0, 0}, {L, 0}}]}, {
              Text["One Period", {(L/2), 0.1}]}}];

You've got to agree that Fourier knew what he was talking about.

B.3.b.ii)

You could use the Mathematica Fit instruction to get the same effect as the FastFourierfit instruction.
Why would anyone want this new FastFourierfit instruction?

When you try to fit using many data points, the Fit instruction can bog down badly.
On the other hand, the FastFourierfit instruction really smokes.
Try it out on another function whose period is 2 π.

\(Clear[f, n, t]; 
f[t_] = E\^Sin[t]; 

    FastFourierfit[f, 2\ π, 4, t]\)

FastFourierfit[f, 2 π, 8, t]

FastFourierfit[f, 2 π, 12, t]

Almost instant responses.

The demands of modern science made it essential to be able to do Fourier fits as fast as possible. In fact, fast Fourier fits continue to play the starring role in most everything electronic including hot fields such as signal analysis, data compression and high definition television.
The ideas behind the workings of a fast Fourier fit instruction have attracted the attention of some of the best mathematical scientists of all time.  None other than Gauss came up with a fast Fourier fit algorithm in 1866.  The modern version of the workings behind the fast Fourier fit sprang out of the 1960's with the work by John Tukey and his cohorts at IBM and Princeton.

Mathematica's workings behind the FastFourierfit instruction depend heavily on Tukey's idea.

B.4) Combining Fourier fit and the Laplace transform to come up with good approximate formulas for periodically forced oscillators

B.4.a)

You are given this periodic forcing function f[t]:

endtime = 15.0; \(Clear[f, t]; 
f[t_] = 5.3\ Sin[3\ \((t\/2 - Floor[t\/2])\)]; 

    \(fplot = 
      Plot[f[t], {t, 0, endtime}, 
        PlotStyle -> {{Thickness[0.02], Blue}}, 
        AxesLabel -> {"t", "f[t]"}]; \)\)

And you are asked to come up with a formula for the solution of
      y''[t]+1.8 y'[t]+8.3 y[t]=f[t]
           with y[0] = 4.7 and y'[0]=-0.9.
First you enter everything:

\(b\  = \ 1.8; \ 
c\  = \ 8.3; \ 
startery\  = \ 4.7; \ 

    starteryprime\  = \ \(-0.9\); \ 
Clear[t, \ y]; 


    oscdiffeq\  = 
      \ {\(\(y'\)'\)[t]\  + \ b\ \(y'\)[t]\  + \ c\ y[t]\  == \ f[t], \ 

        \ \ \ \ y[0]\  == \ startery, \ \(y'\)[0]\  == \ starteryprime}; \ 

    ColumnForm[oscdiffeq]\)

Then you see how the solution plots out:

\(ndsol = 
      NDSolve[oscdiffeq, y[t], {t, 0, endtime}, MaxSteps -> 2000]; 

    Clear[ysol]; 


    ysol[t_] = y[t] /. ndsol[[1]]; 
    ndsplot = 
      Plot[ysol[t], {t, 0, endtime}, 
        PlotStyle -> {{Thickness[0.02], Blue}}, PlotRange -> All, 
        AxesLabel -> {"t", "y[t]"}]; \)

About what you expected.
Now the job is to come up with a formula for this solution.
You try the convolution integral method and Mathematica bogs down. The integrals are just too complicated.
Then you think of the LapalceTransform and experiment with the Laplace transform of f[t]:

LaplaceTransform[f[t], t, s]

Mathematica again balks. 
The integrals involved are just to complicated.
What do you do?

Instead of solving
            y''[t]+1.8 y'[t]+8.3 y[t]=f[t]
           with y[0] = 4.7 and y'[0]=-0.9,
you use Fourier fit to come up with a good sine and cosine wave approximation, approxf[t], of f[t] and then you use the Laplace Transform to solve the nearby oscillator diffeq
            y''[t]+1.8 y'[t]+8.3 y[t]=approxf[t]
           with y[0] = 4.7 and y'[0]=-0.9.
Here you go:
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

Show[fplot];

The period of this function is L = 2.
Get a good Fourier fit of f[t]:

L = 2; n = 3; \(Clear[approxf]; 

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

Nice fit:
Now move in with the Laplace transform to get a formula for the solution of

            y''[t]+1.8 y'[t]+8.3 y[t]=approxf[t]
           with y[0] = 4.7 and y'[0]=-0.9.
Here you go:

\(Clear[f, t, approxY, approxy, s]; 
b = 1.8; c = 8.3; 

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

This is the Laplace transform of the the approximate solution. 
To get the formula for the approximate solution, approxy[t]. you calculate the inverse Laplace transform of approxY[s].
The formula for the solution approxy[t] is:

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

There it is!
A simple approximate formula for the solution.
See it:

approxplot = 
      Plot[approxy[t], {t, 0, endtime}, 
        PlotStyle -> {{Thickness[0.01], Red}}, PlotRange -> All, 
        AxesLabel -> {"t", "y[t]"}];

Compare with what you got from NDSolve:

Show[ndsplot, approxplot];

Both plots are fakes but they are both very good fakes.

B.5)  Fourier integral fit of periodic functions

Here's how you can use integrals to go after a Fourier fit on an interval [0, L]:
Start with the function and its plot on [0, L]:

\(Clear[f, t]; 
f[t_] = Sin[π\ t]\ \((t - 4)\); 
\(L = 4; \)\)Plot[f[t], {t, 0, L}, PlotStyle -> {{Thickness[0.02], Blue}}, 
      AxesLabel -> {"t", "f[t]"}];

To try to fit f[t] on [0, L], you fit with combinations of
        (E^(((I k 2 π t)/L)))
this way:

Clear[A, t];
A[k_]:=N[(((∫_0^L)f[t] (E^(-((I k (2 π) t)/L)))\[DifferentialD]t)/L)];
Clear[k, complexfitter];
m = 5;
complexfitter[t_]=Chop[(∑_(k = -m)^m)A[k] (E^(((I k (2 π) t)/L)))]

\(Clear[realfitter, t]; 

    realfitter[t_] = Chop[ComplexExpand[complexfitter[t]]]\)

Check out the fit with a plot:

Plot[{f[t], realfitter[t]}, {t, 0, L}, 
      PlotStyle -> {{Thickness[0.02], Blue}, {Thickness[0.01], Red}}, 
      AxesLabel -> {"t", ""}];

The periodic nature of the Fourier fit ruins the fit at the ends, but inside [0, L], the fit is not all that bad.

B.5.a)

How do you try to go after a better fit?

The same way you go after a better fast Fourier fit: You increase m.  
Try it:

\(Clear[A, t]; 

    A[k_] := \(A[k] = 
        NIntegrate[
            f[t]\ E\^\(-\(\(I\ k\ \((2\ π)\)\ t\)\/L\)\), {t, 0, L}]\/L
          \); 
Clear[k, complexfitter]; 
m = 13; 

    complexfitter[t_] = 
      Chop[∑\+\(k = \(-m\)\)\%m A[k]\ 
            E\^\(\(I\ k\ \((2\ π)\)\ t\)\/L\)]\)

\(Clear[realfitter, t]; 

    realfitter[t_] = Chop[ComplexExpand[complexfitter[t]]]\)

Check out the fit with a plot:

Plot[{f[t], realfitter[t]}, {t, 0, L}, 
      PlotStyle -> {{Thickness[0.02], Blue}, {Thickness[0.01], Red}}, 
      AxesLabel -> {"t", ""}];

That's almost as good as it can get. If f[0] had been the same as f[L], then the trouble near the endpoints would have disappeared.

B.5.b.i) The formula 
            A[k] = 
        (1/L) 
          (∫_(0 )^L)
               f[t] (E^(-((I k (2 π) t)/L))) 
              dt

What's the idea behind the formula 
          A[k] = 
        (1/L) 
          (∫_0^L) 
            f[t] (E^(-((I k (2 π) t)/L))) 
              dt
as used above?

The idea:
You go for the whole ball of wax in one gulp.
Given a function f[t] with period L, you go ahead and assume that f[t] can be written in a form like this:

m = 3; \(Clear[k, f, complexfitter, A, L]; 

    complexfitter[t_] = 
      ∑\+\(k = \(-m\)\)\%m A[k]\ E\^\(\(I\ k\ 2\ π\ t\)\/L\)\)

Once you go with the assumption that f[t] = complexfitter[t], you gotta agree that
     
            
      (∫_0^L)
           complexfitter[t] (E^(-((I p (2 π) t)/L))) 
          dt    =    
      (∫_0^L) f[t] (E^(-((I p (2 π) t)/L))) 
          dt,
         
no matter what p you go with.
Peek at what this tells you:

\(Clear[p, f]; 

    ColumnForm[
      Table[{Expand[
            ∫\_0\%L
                \( complexfitter[t]\ 
                  E\^\(-\(\(I\ p\ \((2\ π)\)\ t\)\/L\)\)\) 
                \[DifferentialD]t], " must be", 
          ∫\_0\%L
              \( f[t]\ E\^\(-\(\(I\ p\ \((2\ π)\)\ t\)\/L\)\)\) 
              \[DifferentialD]t}, {p, \(-m\), m}]]\)

Your eyes lead you to the formula
           A[k] =  
        (1/L) 
          (∫_(0  )^L)
               f[t] (E^(-((I k (2 π) t)/L))) 
              dt.
This is the formula that was used in part i) so successfully.
Explanation over.

B.5.b.ii)

Not so fast.
Wasn't this supposed to be a math course?
You know very well that a false assumption like 
        f[t] = complexfitter[t]
can be used to explain anything.
How do you justify this false assumption?

Good question.  
The reason it works is that it is almost true.
Reason: For most functions f[t] with period L, only a bureaucratic bean-counter armed with a good magnifying glass could ever tell the difference between f[t] and the complexfitter you get with a very high m.
Try it on
          f[t] = (Sqrt[1 - Cos[π t]])
Because the period of Cos[t] is 2 π, the period of f[t] is L = 2.

\(Clear[f, t, x]; 
f[t_] = \@\(1 - Cos[π\ t]\); 
\(L = 2; \)\)m = 6; \(Clear[k, complexfitter, A]; 

    A[k_] := N[1\/L]\ 
        NIntegrate[
          f[x]\ E\^\(-\(\(I\ k\ \((2\ π)\)\ x\)\/L\)\), {x, 0, L}, 
          AccuracyGoal -> 3]; 

    complexfitter[t_] = 
      Chop[∑\+\(k = \(-m\)\)\%m A[k]\ 
            E\^\(\(I\ k\ \((2\ π)\)\ t\)\/L\)]\)

\(Clear[realfitter, t]; 

    realfitter[t_] = Chop[ComplexExpand[complexfitter[t]]]\)

cycles = 4; Plot[{f[t], realfitter[t]}, {t, 0, L}, 
      PlotStyle -> {{Thickness[0.02], Blue}, {Thickness[0.01], Red}}, 
      AxesLabel -> {"t", ""}, PlotRange -> All];

Tell that miserable bean-counter to get out of your life.

B.5.c.i) The connection between fast Fourier Fit and Fourier integral fit

Go with f[t] = t (3 - t) on [0, L] with L = 3 and compare what you get from Fourier integral fit and fast Fourier fit.

\(Clear[f, t]; 
f[t_] = t\ \((3 - t)\); 
\(L = 3; \)\)\(Clear[A, t]; 

    A[k_] := \(A[k] = 
        N[\(∫\_0\%L
                \( f[t]\ E\^\(-\(\(I\ k\ \((2\ π)\)\ t\)\/L\)\)\) 
                \[DifferentialD]t\)\/L]\); 
Clear[k, complexintegralfitter]; 
    
m = 5; 

    complexintegralfitter[t_] = 
      Chop[∑\+\(k = \(-m\)\)\%m A[k]\ 
            E\^\(\(I\ k\ \((2\ π)\)\ t\)\/L\)]\)

n = m + 1; FastFourierfit[f, L, n, t]

They are almost the same.
Subtract them:

Expand[complexintegralfitter[t] - FastFourierfit[f, L, n, t]]

Look at those puny coefficients reflecting how close the fast Fourier fit is to the Fourier integral fit.

Try it again with a bigger m:

\(Clear[k, complexintegralfitter]; 
m = 10; 

    \(complexintegralfitter[t_] = 
      Chop[∑\+\(k = \(-m\)\)\%m A[k]\ 
            E\^\(\(I\ k\ \((2\ π)\)\ t\)\/L\)]; \)\)n = m + 1; Expand[complexintegralfitter[t] - FastFourierfit[f, L, n, t]]

How do you explain this striking similarity of the two Fourier fits?

The fast Fourier fit tries to fit at equally spaced points on the plot of f[t] on [0, L]. 

When you go with a high n in FastFourierfit[f, L, n, t], you are fitting so many densely packed points that FastFourierfit[f, L, n, t] is almost fitting all the points on the plot of f[t] on [0, L]. 

This nearly replicates the Fourier integral fit of f[t] on [0, L] because the Fourier integral fit tries to fit all the points on the plot of f[t] on [0, L].

B.5.c.ii) Practicalities

Explain the following statement:
"Fast Fourier fit is almost always fast and practical; Fourier integral fit is sometimes slow and sometimes impractical and even impossible."

Answer:

The Achilles heel of the Fourier integral fit is the integral 
          (1/L) 
        (∫_0^L)
             f[t] (E^(-((I k (2 π) t)/L))) 
            dt
that must be done to calculate the coefficients. For some functions f[t], the only practical possibility is to go with NIntegrate. But the highly oscillatory nature of 
        f[t] (E^(-((I k (2 π) t)/L)))
sometimes makes the integrals impractical even for for NIntegrate.

Fast Fourier fit carries none of this disappointing baggage. It runs quickly and accurately, provided no impulses are around.

B.5.c.iii)

Why do folks want to do Fourier integral fit?

If the professors demand that the students work by hand, then the student has little choice between Fourier integral fit and fast Fourier fit.
Reason:
Doing fast Fourier fit by hand is completely out of the question.
Doing integrals by hand is not.

In theoretical situations, most good folks want to use the Fourier integral fit because it gives them a specific formula to work with.  But after the theory is developed and calculations begin, then the same folks usually go with fast Fourier fit.

B.5.d.i) The advantage of fast Fourier fit over Fourier integral fit

Fourier integral fit has been around for more than 200 years. 
Fast Fourier fit has been around less than 50 years and really came into its own with the computer revolution of which you are a part.
Why are many compuer-literate folks so quick to de-emphasize Fourier integral fit in favor of fast Fourier fit?

The Achilles heel of the Fourier integral fit is the integral 
          (1/L) 
        (∫_0^L)
             f[t] (E^(-((I k (2 π) t)/L))) 
            dt
which must be done to calculate the coefficients of the approximation. Even with Mathematica, these integrals can either be very slow or impossible. Numerical integration (via the NIntegrate instruction) is theoretically possible, but the highly oscillatory nature of 
          f[t] (E^(-((I k (2 π) t)/L)))
sometimes makes the integrals impractical even for for numerical integration.

Fast Fourier fit does not care how the function is presented and it doesn't care about doing any integrals. Fast Fourier works rapidly and well for any function - other than impulse hits.
It is robust and fast.
Compare them on the square wave:

\(Clear[f, t]; 
f[t_] = Sign[Sin[π\ t]] + 1; 
\(L = 2; \)\)cycles = 4; fplot = 
      Plot[f[t], {t, 0, cycles L}, 
        PlotStyle -> {{Thickness[0.02], Blue}}, 
        AxesLabel -> {"t", "f[t]"}, 
        PlotLabel -> "cycles" cycles, 
        AspectRatio -> (1/GoldenRatio), 
        Epilog -> {{Red, Thickness[0.02], Line[{{0, 0}, {L, 0}}]}, {
              Text["One Period", {(L/2), 0.1}]}}];

Go for a fast Fourier fit:

\(n = 6; 
Clear[fastapproxf]; 

    fastapproxf[t_] = FastFourierfit[f, L, n, t]\)

Check it out:

\(Clear[realfit]; 
realfit[t_] = Chop[ComplexExpand[fastapproxf[t]]]; 

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

Now do the same thing with the integral Fourier fit:

Try numerical integration via NIntegrate.

\(Clear[A, t]; 

    A[k_] := \(A[k] = 
        N[NIntegrate[
              f[t]\ E\^\(-\(\(I\ k\ \((2\ π)\)\ t\)\/L\)\), {t, 0, L}]\/L]
          \); 
Clear[k, integralfitter]; 
m = 5; 

    integralfitter[t_] = 
      Chop[∑\+\(k = \(-m\)\)\%m A[k]\ 
            E\^\(\(I\ k\ \((2\ π)\)\ t\)\/L\)]\)

This gives an answer, but Mathematica is telling you that it doesn't quite trust the answer.
Check it out:

\(Clear[realfit]; 
realfit[t_] = Chop[ComplexExpand[integralfitter[t]]]; 
    
\(integralfitplot = 
      Plot[{f[t], realfit[t]}, {t, 0, cycles\ L}, 
        PlotStyle -> {{Thickness[0.02], Blue}, {Thickness[0.01], Red}}, 
        AxesLabel -> {"t", ""}, 
        PlotLabel -> "cycles"\ cycles, 
        AspectRatio -> 1\/GoldenRatio, 
        Epilog -> {{Red, Thickness[0.02], Line[{{0, 0}, {L, 0}}]}, {
              Text["One Period", {L\/2, 0.1}]}}]; \)\)

The calculation of the integral fit took a lot longer than the calculation of the fast fit.
Fast Fourier fit went through this like a hot knife in butter. 
Fourier integral fit went through this like a butter knife cutting steak.

B.5.d..ii)

If you are working with hand calculations only, why is Fourier integral fit your only real choice?

Fast Fourier fit involves solving many linear equations. Solving these equations by hand is just out of the question. 
This leaves you no choice but Fourier integral fit. That's why you often see the Fourier integral fit in the old-fashioned courses that emphasize hand calculation.