MGM.11 Function Spaces and Root-Mean-Square Approximation
TUTORIALS

Needs["LinearAlgebra`MatrixManipulation`"];
Needs["LinearAlgebra`Orthogonalization`"];
Needs["Calculus`Master`"];
Off[Unique::"usym"];
Needs["Graphics`Colors`"];
Off[General::"spell"];
Off[General::"spell1"];
Off[ParametricPlot3D::"ppcom"];
Off[ParametricPlot::"ppcom"];
Off[Plot::"plnr"];
SetOptions[Limit, Analytic -> True];
$Pre = Chop;

CMView = {2.7, 1.6, 1.2};

Vector Drawer

\(\(\( (*\ \(\(:\)\(Date : \ \ Copyright\ 2001\)\), \ 
      Math\ Everywhere, \ \(\(Inc\)\(.\)\)\ 
*) \)\(
\)\( (*
Original\ Conception\ to\ create\ a\ vector\ drawer\ is\ due\ to\ Bill\ Davis . 
Current\ code\ is\ by\ Bruce\ \(\(Carpenter\)\(.\)\)
*) \)\(
\)\( (*\ \ \(\(\(:\)\(Drawbacks : \ \ Since\ the\ vectors\ are\ drawn\ out\ of\)\)\ 

      \ \ \ \ \ context\), \ 
      the\ shapes\ of\ the\ arrowheads\ are\ only\ 
\ \ \ \ \ correct\ with\ AspectRatio\ set\ to\ \(\(Automatic\)\(.\)\)
*) \)\)\)

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]==2;

Perpend[a_?VectorQ]:=
    Normalize[If[a[[1]]==0,{1,0,0},{a[[2]],-a[[1]],0}]]/;
      Length[a]==3;
      
Base3D = Table[{Cos[2 k Pi/8.], Sin[2 k Pi/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]==2

Vector[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[#,tip]]& /@ Partition[(tip - hr diff + #)&/@((ap hr Norm[diff] Base3D) . {temp,Normalize[Cross[temp,diff]]}),2,1]]}}]
] /; Length[from]==Length[to]==3


Arrow[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[];
Null

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

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].";

Needs["Statistics`ContinuousDistributions`"];

Clear[matrixplucker];
matrixplucker[expr_, 
        vars_] := Simplify[
          Outer[D, {expr}, 
            vars]][[1]]

T.1) Scaling for other intervals

T.1.a.i) Polynomial approximation on intervals [-L,L] with L ≠ 1

Here is a function f[t] plotted on [-L,L] with L = 2.7

L =  2.7;
Clear[f, t];
f[t_]  =  Sin[t]  +  (2/(1.3  +  t^2));
fplot = Plot[f[t],{t, -L,  L},PlotStyle -> {{Thickness[0.01], Blue}},AxesLabel -> {"t", "f[t]"}];

Come up with a pretty good root-mean square approximation of f[t] on [-L.L] by polynomials.

If the given interval had been [-1,1], this would have been just a matter of copying, pasting and editing from the Basics. 
To take advantage of this fact, you just scale f[t] on [-L,L] this way:

Clear[scaledf];
scaledf[t_] = f[t  L]
scaledfplot = Plot[scaledf[t],{t, -1, 1},PlotStyle -> {{Thickness[0.01], Blue}},AxesLabel -> {"t", "scaledf[t]"}];

The formula 
          scaledf[t] = f[t L]
tells you that as t advances from -1 to 1, scaledf[t] does the same thing f[t] does as t advances from -L to L.

Show[GraphicsArray[{fplot, scaledfplot}]];

Now approximate scaledf[t] with polynomials on [-1,1]:

khigh  =  4;
Clear[s, k, t];
(Subscript[s, k])_[t_]  :=  LegendreP[k - 1,  t];
a  =  -1;
\(\(b\  = \ 1;\)
\)
Clear[fouriercoeff, Sclosest, k];
fouriercoeff[k_]  :=  
        NIntegrate[scaledf[t] (Subscript[s, k])[t], {t, a, b}]/NIntegrate[(Subscript[s, k])[t] (Subscript[s, k])[t], {t, a, b}];
Sclosest[
        t_]   =  (∑_(k = 1)^khigh) fouriercoeff[k] (Subscript[s, k])[
            t]

Check the quality of the approximation:

scaledfitplot = Plot[{scaledf[t], Sclosest[t]},{t, -1, 1},PlotStyle -> {{Thickness[0.02], Blue}, {Thickness[0.01], 
                  Red}},AxesLabel -> {"t", ""},Epilog -> {{Blue, 
                  Text["scaledf[t]", {a  +  0.5, 
                      f[a  +  0.51]}]}, {Red, 
                  Text["Sclosest[t]", {b  -  0.5, 
                      Sclosest[b - 0.5]}]}}];

Not great.
Use more Legendre Polynomials to go after a better approximation:

khigh  =  12;
Clear[s, k, t];
(Subscript[s, k])_[t_]  :=  
        LegendreP[k - 1,  t];
a  =  -1;
b  =  1;

Clear[fouriercoeff, Sclosest, 
        k];
fouriercoeff[k_]  :=  
        NIntegrate[scaledf[t] (Subscript[s, k])[t], {t, a, b}]/NIntegrate[(Subscript[s, k])[t] (Subscript[s, k])[t], {t, a, b}];
Sclosest[
        t_]   =  (∑_(k = 1)^khigh) fouriercoeff[k] (Subscript[s, k])[
            t]

Check the quality of the approximation:

scaledfitplot = Plot[{scaledf[t], Sclosest[t]},{t, -1, 1},PlotStyle -> {{Thickness[0.02], Blue}, {Thickness[0.01], 
                  Red}},AxesLabel -> {"t", ""},Epilog -> {{Blue, 
                  Text["scaledf[t]", {a  +  0.5, 
                      f[a  +  0.5]}]}, {Red, 
                  Text["Sclosest[t]", {b  -  0.5, 
                      Sclosest[b - 0.5]}]}}];

Lookin' pretty good.
Now unscale Sclosest[t] via 
            unscaledSclosest[t]=Sclosest[((t )/L)]:

Clear[unscaledSclosest];
unscaledSclosest[t_]=Sclosest[((t )/L)]

See the fit:

fitplot = Plot[{f[t], unscaledSclosest[t]},{t, -L, L},PlotStyle -> {{Thickness[0.015], 
                  Blue}, {Thickness[0.01], Red}},AxesLabel -> {"t", ""}];

Just fine!.

T.1.a.ii) Sine wave approximation on intervals [0,L] with L ≠ π

Here is a function f[t] plotted on [0,L] with L = 4:

L = 4
Clear[f, t];
f[t_] = t ((6.2/(1.3 + (E^t))) - 0.4);
fplot = 
        Plot[f[t], {t, 0, L}, PlotStyle -> {{Thickness[0.01], Blue}}, 
          AxesLabel -> {"t", "f[t]"}];

Come up with a pretty good root-mean square approximation of f[t] on [0,L] by sine waves.

If the given interval had been [0,π], this would have been just a matter of copying, pasting and editing from the Basics. 
To take advantage of this fact, you just scale f[t] on [0,L] this way:

Clear[scaledf];
scaledf[t_]=f[t(L/π)];
scaledfplot = 
        Plot[scaledf[t], {t, 0, π}, 
          PlotStyle -> {{Thickness[0.01], Blue}}, 
          AxesLabel -> {"t", "scaledf[t]"}];

The formula 
          scaledf[t]=f[t(L/π)]
tells you that as t advances from 0 to π, scaledf[t] does the same thing as f[t] does as t advances from 0 to L.

Show[GraphicsArray[{fplot, scaledfplot}]];

Now approximate scaledf[t] with Sine waves on [0,π]:

khigh  =  15;
Clear[s, k, t];
(Subscript[s, k])_[
        t_]  :=  Sin[k t]
a  =  0;
b =π;

Clear[fouriercoeff, Sclosest, 
        k];
fouriercoeff[k_]  :=  
        NIntegrate[scaledf[t] (Subscript[s, k])[t], {t, a, b}]/NIntegrate[(Subscript[s, k])[t] (Subscript[s, k])[t], {t, a, b}];
Sclosest[
        t_]   =  (∑_(k = 1)^khigh) fouriercoeff[k] (Subscript[s, k])[
            t]

Check the quality of the approximation:

scaledfitplot  =  
        Plot[{scaledf[t], Sclosest[t]}, {t, a, b}, 
          PlotStyle -> {{Thickness[0.01], Blue}, {Thickness[0.01], Red}}, 
          AxesLabel -> {"t", ""}, 
          Epilog -> {{Blue, 
                Text["scaledf[t]", {a  +  1, f[a  +  1]}]}, {Red, 
                Text["Sclosest[t]", {b  -  1, 
                    Sclosest[b - 1]}]}}];

Not bad - except at the right endpoint.
Now unscale Sclosest[t] via 
            unscaledSclosest[t]=Sclosest[((t π)/L)]:

Clear[unscaledSclosest];
unscaledSclosest[t_] = Sclosest[t  (π/L)]

See the fit:

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

Not bad except at the right endpoint.

T.2) Weighted root-mean -square -distance between two functions
       Weighted dot product of two functions
       Families orthogonal with respect to a weight function
       Chebyshev polynomials
       Best weighted root-mean-square -distance  approximation

T.2.a.i) Examples of weight functions

A weight function on an interval [a,b] is any function weight[t] with 
weight[t]  ≥  0 for all t's with a ≤ t ≤ b.
Here's an example of a weight function on [-2,2]

a  =  -2; 
b  =  2;
Clear[weight, t];
\(\(weight[t_]\  = \ 5\ t\^2\ E\^\(-t\^2\);\)
\)
Plot[weight[t], {t, a, b}, PlotStyle -> {{Thickness[0.01], Red}}, 
        AxesLabel -> {"t", "weight[t]"}, 
        AspectRatio -> 1/GoldenRatio];

This function puts no weight on t = 0 and puts heavier weight on what's going on around t = 1 and t = -1.

Here's an example of a weight function on [-1,1]

a  =  -1; 
b  =  1;
Clear[weight, t];
\(\(weight[t_]\  = \ \@\(1 - t\^2\);\)
\)
Plot[weight[t], {t, a, b}, PlotStyle -> {{Thickness[0.01], Red}}, 
        AxesLabel -> {"t", "weight[t]"}, 
        AspectRatio -> 1/GoldenRatio];

This function puts no weight on t = -1 and no weight on t = 1 but puts heavier weight on what's going on around t =0.

T.2.a.ii) Weighted root-mean-square distance between two functions on an interval.

.Here's  a weight function on [a,b] with  a = -1 and b = 1:

a = -1; 
b = 1; 
Clear[weight, t]; 
weight[t_] = (Sqrt[1 - (t^2)]); 
Plot[weight[t], {t, a, b}, PlotStyle -> {{Thickness[0.01], Red}}, 
      AxesLabel -> {"t", "weight[t]"}, 
      AspectRatio -> (1/GoldenRatio)];

Here are two functions f[t] and g[t]  on [a,b] with  a = -1 and b = 1:

Clear[f, g, t];
f[t_]  = Cos[3 t]^2;
g[t_]  =  Sin[2 t]^2 - 0.3;
Plot[{f[t], g[t]}, {t, a, b}, 
        PlotStyle -> {{Thickness[0.01], CadmiumOrange}, {Thickness[0.01], 
              Black}}, AxesLabel -> {"t", ""}];

Calculate the weighted root-mean-square distance from f[t] to g[t] with respect to the given weight function on [a,b].

The weighted root-mean-square distance from f[t] to g[t] with respect to this weight function on [a,b] is  
       (Sqrt[(∫_a^b)(((f[t] - g[t]))^2)weight[t] dt]):

(Sqrt[NIntegrate[(((f[t] - g[t]))^2) weight[t],{t, a, b}]])

See what happens when you go with new functions f[t] and g[t]:

Clear[f, g, t];
m  =  Random[Real, {1, 3}];
n  =  Random[Real, {3, 6}];
f[t_]  =  Cos[m t]^2;
\(\(g[t_]\  = \ Cos[n\ t]\^2;\)
\)
Plot[{f[t], g[t]}, {t, a, b}, 
        PlotStyle -> {{Thickness[0.01], CadmiumOrange}, {Thickness[0.01], 
              Black}}, AxesLabel -> {"t", ""}];

The weighted root-mean-square distance from f[t] to g[t] with respect to this weight function on [a,b] is  
       (Sqrt[(∫_a^b)(((f[t] - g[t]))^2) weight[t] 
            dt]):

(Sqrt[NIntegrate])[(((f[t] - g[t]))^2) weight[t], {t, a, b}]

T.2.a.iii) Weighted dot product of  two functions on an interval

.Here's  a weight function on [a,b] with  a = -1 and b = 1:

a  =  -1; 
b  =  1;
Clear[weight, t];
\(\(weight[t_]\  = \ Sqrt[1\  + \ t];\)
\)
Plot[weight[t], {t, a, b}, PlotStyle -> {{Thickness[0.01], Red}}, 
        AxesLabel -> {"t", "weight[t]"}, 
        AspectRatio -> 1/GoldenRatio];

Here are two functions f[t] and g[t]  on [a,b] with  a = -1 and b = 1:

Clear[f, g, t];
f[t_] = (E^(-0.4 t)) - 1.0;
\(\(g[t_] = E\^\(\(-2.7\)\ t\^2\);\)
\)
Plot[{f[t], g[t]}, {t, a, b}, 
        PlotStyle -> {{Thickness[0.01], Orange}, {Thickness[0.01], Magenta}}, 
        AxesLabel -> {"t", ""}];

Calculate weighted dot product  f•g   of f[t] and g[t] on [a,b] with respect to the given weight function.

Just calculate:
              f•g =(∫_a^b)f[t] g[t]weight[t]dt:

NIntegrate[f[t] g[t] weight[t] ,{t,  a, b}]

They are almost perpendicular with respect to the weighted dot product.

T.2.b.i) Weighted root-mean-square approximation via Chebyshev polynomials system

The set of functions   
        {(Subscript[s, 1])[t], (Subscript[s, 2])[t], (Subscript[s, 3])[t], . . ., 
        (Subscript[s, k])[t],  . . .}
for
         (Subscript[s, k])[t]  =  ChebyshevT[k  - 1,  t] 
is an orthogonal set with respect to the weight function 
          weight[t] = 1/(Sqrt[1 - (t^2)])
on the interval a ≤ t ≤ b with a = -1 and b = 1.
Here are the formulas for the first seven:

Clear[k];
ColumnForm[Table[ChebyshevT[k - 1,  t], {k, 1, 7}]]

Here are plots of the first six:

Clear[s, g, k, t];
\(\(s\_k_[t_] := ChebyshevT[k - 1, t];\)
\)
a = -1;
\(\(b = 1;\)
\)
Plot[{(Subscript[s, 1])[t], (Subscript[s, 2])[t], (Subscript[s, 3])[t], (Subscript[s, 4])[t], (Subscript[s, 5])[t], (Subscript[s, 6])[t]}, {t, a, 
          b}, PlotStyle -> {{Black}, {Blue}, {Magenta}, {Red}, {Green}, {CadmiumOrange}}, 
        PlotLabel -> Chebyshev Polynomials on[a, b]];

To illustrate the fact that this set of functions   
        {(Subscript[s, 1])[t], (Subscript[s, 2])[t], (Subscript[s, 3])[t], . . ., 
        (Subscript[s, k])[t],  . . .}
is an orthogonal set with respect to the weight function 
          weight[t] = (1/(Sqrt[1 - (t^2)]))
on the interval a ≤ t ≤ b with a = -1 and b = 1,
look at this calculation of
              (∫_a^b) (Subscript[s, k])[t] (Subscript[s, p])[t] weight[t] dt
for random positive integers p and k with p ≠k

Clear[weight];
weight[t_] = 
        (1/(Sqrt[1 - (t^2)]));
k = 
        Random[Integer, {1, 8}];
p = 
        Random[Integer, {k + 1, 2 k}];
(∫_a^b)(Subscript[s, k])[t] (Subscript[s, p])[t] weight[t]\[DifferentialD]t

Here is a function f[t] plotted on [a,b] = [-1,1]:

Clear[f, t];
a  =  1;
b  =  -1;
\(\(f[t_]\  = \ Sin[π\ t]\  + \ 0.4\ Cos[π\ t];\)
\)
fplot  =  
        Plot[f[t], {t, a, b}, PlotStyle -> {{Thickness[0.01], Indigo}}, 
          AxesLabel -> {"t", "f[t]"}];

Go with the function space S[a,b] spanned by:

khigh  =  5;
Table[ (Subscript[s, k])[t], {k, 1, khigh}]

Come up with the function Sclosest[t] in  S[a,b]  so that
the weighted root-mean-square distance between f[t] and Sclosest is as small as possible

It's not bad at all.
For each k you take the weighted component of f[t] in the direction of (Subscript[s, k])[t]:
           f•((Subscript[s, k])/((Subscript[s, k]) •(Subscript[s, k]))) (Subscript[s, k])[t]

(This uses the weighted dot productf•g = 
        ∫(∫_a^b) f[t] g[t] weight[t] 
              dt.)

and then you add 'em up to get
              Sclosest[t]=(∑_(k = 1)^khigh) f•((Subscript[s, k])/((Subscript[s, k]) •(Subscript[s, k])))(Subscript[s, k])[t]
See the formula:

Clear[weightedfouriercoeff, Sclosest, k];
weightedfouriercoeff[k_] := ((NIntegrate[f[t] (Subscript[s, k])[t] weight[t],{t, a, b}])/(NIntegrate[(Subscript[s, k])[t] (Subscript[s, k])[t] weight[t],{t, a, b}]));
Sclosest[
        t_]   =  (∑_(k = 1)^khigh) weightedfouriercoeff[k] (Subscript[s, k])[
            t]

See the quality of the approximation of  f[t] by Sclosest[t]:

fitplot  =  
        Plot[{f[t], Sclosest[t]}, {t, a, b}, 
          PlotStyle -> {{Thickness[0.015], Indigo}, {Thickness[0.01], Gold}}, 
          AxesLabel -> {"t", ""}, 
          Epilog -> {{Indigo, 
                Text["f[t]", {a  +  0.5, f[a  +  0.5]}]}, {Gold, 
                Text["Sclosest[t]", {b  -  0.5, 
                    Sclosest[b - 0.5]}]}}];

T.2.b.ii) Increasing the quality of the approximation

Keep everything the same as in part i), but go with the function space S[a,b] spanned by more Chebyshev polynomials:

khigh  =  8;
Table[ (Subscript[s, k])[t], {k, 1, khigh}]

Come up with the function Sclosest[t] in  S[a,b]  so that
the weighted root-mean-square distance between f[t] and Sclosest is as small as possible

Clear[weightedfouriercoeff, Sclosest, k];
weightedfouriercoeff[k_] := ((NIntegrate[f[t] (Subscript[s, k])[t] weight[t],{t, a, b}])/(NIntegrate[(Subscript[s, k])[t] (Subscript[s, k])[t] weight[t],{t, a, b}]));
Sclosest[
        t_]   =  (∑_(k = 1)^khigh) weightedfouriercoeff[k] (Subscript[s, k])[
            t]

See the quality of the approximation of  f[t] by Sclosest[t]:

fitplot  =  
        Plot[{f[t], Sclosest[t]}, {t, a, b}, 
          PlotStyle -> {{Thickness[0.02], Indigo}, {Thickness[0.01], Gold}}, 
          AxesLabel -> {"t", ""}, 
          Epilog -> {{Indigo, 
                Text["f[t]", {a  +  0.5, f[a  +  0.5]}]}, {Magenta, 
                Text["Sclosest[t]", {b  -  0.5, 
                    Sclosest[b - 0.5]}]}}];

Sharing ink all the way..
Great!
The weighted root-mean-square distance between f[t] and Sclosest[t] is:

(Sqrt[NIntegrate])[(((f[t] - Sclosest[t]))^2) weight[t], {t, a, b}]

Love those Chebyshev polynomials.

T.3)  Fourier Sine approximation and the heat equation
    Fourier Sine approximation and the wave equation

T.3.a.i) Fourier sine approximation and the heat equation

Start with a heated wire L units long with the temperature allowed to vary from position to position on the wire. Think of the wire as the interval 0 ≤ x ≤ L.  At the start of the experiment,  you instantly cool the ends at x = 0 and x = L  and maintain these ends at temperature 0, and you take pains to guarantee that the rest of the wire is perfectly insulated.
At the start of this particular experiment, the temperature of the wire at position x (for 0 ≤ x ≤ L = π) is given by the following function startertemp[x]:

L = π;
Clear[startertemp, x];
\(\(startertemp[x_] = 
        x\ \((x - π\/4)\)\ \((x - 2)\)\ \((π - x)\);\)
\)
starterplot = 
        Plot[startertemp[x], {x, 0, L}, 
          PlotStyle -> {{Thickness[0.015], Blue}}, 
          AxesLabel -> {"x", "starting temperature"}];

To fully understand this plot, look at this:

xvalue = 2.5;
Clear[pointer];
\(\(pointer[x_] := 
        Arrow[{0, startertemp[x]}, Tail -> {x, 0}, 
          VectorColor -> Red];\)
\)
Show[starterplot, pointer[xvalue]];

Think of the interval [0,L] = [0 π] as the wire.  
The tip of the pointer tells you the starting temperature (at time t = 0) at the tail of the pointer.
Your problem here is to use Fourier Sine approximation on [0,π] of startertemp[x] come up with a function temp[x,t] that estimates the  temperature of the wire at position x at time t after the experiment begins. 
Do it.

First do a Sine approximation of f[x] = startertemp[x] on [0,Pi]:

khigh  =  3;
Clear[s, g, f, k, 
        x];
f[x_]  =  startertemp[x];
(Subscript[s, k])_[x_] := Sin[ k x];
a  = 0;
b  = π;

Clear[fouriercoeff, Sclosest, 
        k];
fouriercoeff[k_]  :=  fouriercoeff[k]  =  
          NIntegrate[f[t] (Subscript[s, k])[t], {t, a, b}]/NIntegrate[(Subscript[s, k])[t] (Subscript[s, k])[t], {t, a, b}];
Sclosest[x_]  = (∑_(k = 1)^khigh)fouriercoeff[k] (Subscript[s, k])[x]

Check out the quality of the fit:

fitplot = 
        Plot[{startertemp[x], Sclosest[x]}, {x, a, b}, 
          PlotStyle -> {{Thickness[0.02], Blue}, {Thickness[0.01], 
                Gold}}, AxesLabel -> {"x", ""}];

That's not the greatest fit. 
Make khigh bigger:

khigh  =  8;
Clear[s, g, f, k, 
        x];
f[x_]  =  startertemp[x];
(Subscript[s, k])_[x_] := Sin[ k x];
a  = 0;
b  = π;

Clear[fouriercoeff, Sclosest, 
        k];
fouriercoeff[k_]  :=  fouriercoeff[k]  =  
          NIntegrate[f[t] (Subscript[s, k])[t], {t, a, b}]/NIntegrate[(Subscript[s, k])[t] (Subscript[s, k])[t], {t, a, b}];
Sclosest[x_]  = (∑_(k = 1)^khigh)fouriercoeff[k] (Subscript[s, k])[x]

Check out the quality of the fit:

fitplot = 
        Plot[{startertemp[x], Sclosest[x]}, {x, a, b}, 
          PlotStyle -> {{Thickness[0.02], Blue}, {Thickness[0.01], 
                Gold}}, AxesLabel -> {"x", ""}];

Beautiful fit - go with this one.
To get temp[x,t], you look at:

Sclosest[x]

Pick off the coefficients of the Sin[k x] terms:

Clear[A, k];
\(\(A[k_] := Coefficient[Sclosest[x], Sin[k\ x]];\)
\)
Table[A[k], {k, 1, khigh}]

Now you're done because you can write down temp[x,t].
It's just:

Clear[temp, x, t];
temp[x_, t_]=(∑_(k = 1)^khigh)A[k] (E^(-(k^2) t))Sin[k x]

See the wire cool down:

Clear[tempplotter];
\(\(tempplotter[t_] := 
        Plot[temp[x, t], {x, 0, π}, 
          PlotStyle -> {{Thickness[0.01], Blue}}, 
          PlotRange -> {\(-1\), 2}, 
          AxesLabel -> {"x", " temperature"}, 
          PlotLabel -> N[t]\ "= t", 
          AspectRatio -> 1\/2];\)
\)
\(\(timejump = 0.2;\)
\)
Table[tempplotter[t], {t, 0, 1.4, timejump}]

It takes about 1.2 time units for the wire to cool all the way down to 0.

T.3.a.ii) The heat equation explains why you want the exponentials

You get temp[x,t] by taking a rigged sine fit of starterposition[x]:

Sclosest[x]

And you get temp[x,t] by inserting Cosines into each term of the rigged sine fit:

temp[x, t]

Why did that trick of inserting the exponentials work?

Engineering studies have shown that after the appropriate unit adjustments are made, the function temp[x,t] satisfies the  partial differential equation known as the heat equation
        \(∂\_{x, 2}temp[x, t] = ∂\_t temp[x, t]\)

with 
->  temp[x, 0] = Sclosest[x]
 and
->  temp[0, t]=0  and temp[L, t] = 0 for all t's.
Inserting the exponentials made tempt[x,t] so all this:

Check
         \(∂\_{x, 2}temp[x, t] = ∂\_t temp[x, t]\):

\(∂\_{x, 2}temp[x, t]\)

(∂_t) temp[x, t]

Good.
Check the rest:

temp[x, 0]

Sclosest[x]

temp[0, t]

temp[π, t]

Done!

T.3.b.i) Fourier Sine approximation and the wave equation

The ends of a guitar string are anchored at 0 and L on the x-axis and the string is pulled to an initial position and then allowed to vibrate on its own starting with initial velocity 0.

At the start of this particular experiment, the position of the guitar string at position x (for 0 ≤ x ≤  L = Pi) is given by the following function starterposition[x]:

L = π;
Clear[starterposition, x];
\(\(starterposition[x_] = 
        x\ \((x - π\/4)\)\ \((x - 2)\)\ \((π - x)\);\)
\)
starterplot = 
        Plot[starterposition[x], {x, 0, L}, 
          PlotStyle -> {{Thickness[0.01], Blue}}, PlotRange -> All, 
          AspectRatio -> (1/4), 
          AxesLabel -> {"x", "starting position"}];

Think of the curve as the starting position of the guitar string.
To fully understand this plot, look at this:

xvalue = 2.5;
Clear[pointer];
\(\(pointer[x_] := 
        Arrow[{0, starterposition[x]}, Tail -> {x, 0}, 
          VectorColor -> Red];\)
\)
Show[starterplot, pointer[xvalue]];

The tip of the pointer tells you the starting position (at time t = 0) at the tail of the pointer.
Your problem here is to  use Fourier Sine approximation on [0,π] of starterposition[x] come up come up with a function position[x,t] that estimates the position of the guitar string at position x on the x-axis at time t after the experiment begins. 
Do it.

First do a Sine approximation of f[x] = starterposition[x] on [0,Pi]:

khigh  =  8;
Clear[s, g, f, k, x];
f[x_] = starterposition[x];
(Subscript[s, k])_[x_] := Sin[ k x];
a  = 0;
b  =  Pi;

Clear[fouriercoeff, Sclosest, 
        k];
fouriercoeff[k_]  :=  fouriercoeff[k]  =  
          NIntegrate[f[t] (Subscript[s, k])[t], {t, a, b}]/NIntegrate[(Subscript[s, k])[t] (Subscript[s, k])[t], {t, a, b}];
Sclosest[x_]  = (∑_(k = 1)^khigh)fouriercoeff[k] (Subscript[s, k])[x]

Check out the quality of the fit:

fitplot = 
        Plot[{starterposition[x], Sclosest[x]}, {x, a, b}, 
          PlotStyle -> {{Thickness[0.02], Blue}, {Thickness[0.01], 
                Gold}}, AxesLabel -> {"x", ""}];

Beautiful fit - go with this one.
To get position[x,t], you look at:

Sclosest[x]

Pick off the coefficients of the Sin[k x] terms:

Clear[A, k];
\(\(A[k_] := Coefficient[Sclosest[x], Sin[k\ x]];\)
\)
Table[A[k], {k, 1, khigh}]

Now you're done because you can write down position[x,t].
It's just:

Clear[position, x, t];
position[x_, t_]=(∑_(k = 1)^khigh)A[k] Cos[k t]Sin[k x]

See the string vibrate::

Clear[tempplotter];
\(\(tempplotter[t_] := 
        Plot[position[x, t], {x, 0, π}, 
          PlotStyle -> {{Thickness[0.01], Blue}}, 
          PlotRange -> {\(-2\), 2}, 
          AxesLabel -> {"x", " position"}, 
          PlotLabel -> N[t]\ "= t", 
          AspectRatio -> 1\/2];\)
\)
\(\(timejump = π\/8;\)
\)
Table[tempplotter[t], {t, 0, 2 π, timejump}]

It takes  2 π time units for the string to return to its original position..

T.3.b.ii) A couple of calculations

Keep everything the same as above and look at these calculations:

(∂_t) position[x, t] /. t -> 3

(∂_t) position[x, t] /. x -> 2

What do these calculations measure?

Take them one at a time:
The velocity of the string at x at time t = 3 is:

(∂_t) position[x, t] /. t -> 3

See it :

Plot[(∂_t) position[x, t] /. t -> 3, {x, 0, Pi}, 
      PlotStyle -> {{Red, Thickness[0.01]}}, 
      AxesLabel -> {"x", "velocity when t = 3"}];

The velocity of the string at x = 2 at time t is:

(∂_t) position[x, t] /. x -> 2

See it :

Plot[Evaluate[(∂_t) position[x, t] /. x -> 2],{t,0,4  Pi},PlotStyle -> {{Red, Thickness[0.01]}},AxesLabel -> {"t", "velocity at x = 2"}];

T.3.b.iii) The wave equation explains why you want the Cosines

You get position[x,t] by taking a rigged sine fit of starterposition[x]:

Sclosest[x]

And you get position[x,t] by inserting Cosines into each term of the rigged sine fit:

position[x, t]

Why do you insert those Cosines?

Engineering studies have shown that after the appropriate unit adjustments are made, the function position[x,t] satisfies the wave equation
                  \(∂\_{x, 2}position[x, t] = 
        ∂\_{t, 2}position[x, t]\).

with 
-> position[x, 0] = Sclosest[x]
-> position[0, t] = 0 andposition[π, t] = 0 for all t's

and with 
-> (∂_t) position[x, t]=0 for t = 0 .

Inserting the Cosines made position[x,t] do all this:

Check \(∂\_{x, 2}position[x, t] = 
        ∂\_{t, 2}position[x, t]\):

\(∂\_{x, 2}position[x, t]\)

\(∂\_{t, 2}position[x, t]\)

This checks.
Check: position[x, 0] = Sclosest[x]

position[x, 0]

Sclosest[x]

This checks.
Check:  position[0, t] = 0 and position[Pi, t] = 0 for all t's

position[0, t]

position[Pi, t]

This checks.
Check:  (∂_t) position[x, t] = 0 for t = 0 .

(∂_t) position[x, t] /. t -> 0

T.4) Using Fourier methods to try to bring the Dirac Delta function to life

T.4.a)

For a lot of folks the Dirac Delta function DiracDelta[t] is a total mystery.
If t ≠ 0, then DiracDelta[t] = 0:

t  =  Random[Real, {-4, 4}];
DiracDelta[t]

If t = 0, then DiracDelta[t] is not defined as a real number:

t  =  0;
DiracDelta[t]

Yet,
               (∫_a^b) DiracDelta[t] dt = 1
 provided a < 0 < b.

a  =  Random[Real, {-2, -0.0001}];
b  =  Random[Real, {0.0001, 2}];
Clear[t];
(∫_a^b) DiracDelta[t] \[DifferentialD]t

But   
          (∫_a^b) DiracDelta[t] dt = 0 
for 0 < a < b:

a  =  Random[Real, {0.00001, 3}];
b  =  Random[Real, {a, 3}];
(∫_a^b) DiracDelta[t] \[DifferentialD]t

And  
        (∫_a^b) DiracDelta[t] dt = 0 
for a < b < 0:

a  =  Random[Real, {-2, -0.0001}];
b  =  Random[Real, {a, -0.0001}];
(∫_a^b) DiracDelta[t] \[DifferentialD]t

Go with the orthogonal system of Sines and Cosines on [-π,π]; and see what you get when you try to approximate DiracDelta[t] with Sine and Cosine waves::

Clear[f, t];
f[t_]  =  DiracDelta[t];
a = -π;
b = π;
khigh  =  9; 
\(\(Clear[s, k, t];\)
\)
(Subscript[s, k])_[t_] := 
        If[EvenQ[k], Sin[((k t)/2)], 
          Cos[(1/2) (k - 1) t]];
Clear[fouriercoeff, SinCosapprox, k];
fouriercoeff[
          k_]  :=   fouriercoeff[
            k]  =  (((∫_a^b) f[t] (Subscript[s, k])[t] \[DifferentialD]t)/((∫_a^b) (Subscript[s, k])[t] (Subscript[s, k])[t] \[DifferentialD]t));
\(SinCosapprox[
        t_]\ \  = \ ∑\+\(k = 1\)\%khigh\(\( fouriercoeff[
            k]\)\(\ \)\(s\_k[t]\)\(
\)\)\)

See it:

Plot[SinCosapprox[t], {t, a, b}, AspectRatio -> 2, PlotRange -> All, 
        PlotStyle -> {{Thickness[0.01], Red}}];

Kick up khigh:

khigh  =  14; 
Clear[SinCosapprox];
\(SinCosapprox[
        t_]\ \  = \ ∑\+\(k = 1\)\%khigh\(\( fouriercoeff[
            k]\)\(\ \)\(s\_k[t]\)\(
\)\)\)

See it:

Plot[SinCosapprox[t], {t, a, b}, AspectRatio -> 2, PlotRange -> All, 
        PlotStyle -> {{Thickness[0.01], Red}}];

Kick up khigh even more:

khigh  =  30; 
Clear[SinCosapprox];
\(SinCosapprox[
        t_]\ \  = \ ∑\+\(k = 1\)\%khigh\(\( fouriercoeff[
            k]\)\(\ \)\(s\_k[t]\)\(
\)\)\)

See it:

Plot[SinCosapprox[t], {t, a, b}, AspectRatio -> 2, PlotRange -> All, 
        PlotStyle -> {{Thickness[0.01], Red}}];

Kick up khigh even more:

khigh  =  100; 
Clear[SinCosapprox];
\(\(SinCosapprox[
          t_]\ \  = \ ∑\+\(k = 1\)\%khigh fouriercoeff[k]\ s\_k[
              t];\)
\)

See it:

Plot[SinCosapprox[t], {t, a, b}, AspectRatio -> 2, PlotRange -> All, 
        PlotStyle -> {{Thickness[0.01], Red}}];

That's looking the way you might expect a decent approximation of DiracDelta[t] to look.
And these approximations are settling into a pattern.
\(\(\(\ \ \ \ \ \ \ 	\)\(DiracApprox[khigh, t]\  = 
        1\/\(2\ π\) + \(1\/π\) \(∑\+\(k = 1\)\%khigh Cos[
                k\ t]\)\)\)\):


When you integrate 
\(\(\(	\)\(DiracApprox[khigh, t]\  = 
        1\/\(2\ π\) + \(1\/π\) \(∑\+\(k = 1\)\%khigh Cos[
                k\ t]\)\)\)\) 
 from - h to h, you get

     (∫_(-h)^h) DiracApprox[khigh, t] dt =
                \(\(2  h\)\/\(2\ π\) + \(\(1\/π\) \(∑\+\(k = 1\)\%khigh\( Sin[k\ h] - Sin[k\ \((\(-h\))\)]\)\/k\)\(		\)\)\)
                 = \(\(1\/π\) \((h\  + \ ∑\+\(k = 1\)\%khigh\( 2  Sin[k\ h]\)\/k)\)\)
 See how
         (∫_(-h)^h) DiracApprox[khigh, t] dt 
turns out for h = 0.01 and various choices of khigh.

h  =  0.01;
khigh  =  50;
N[(1/π) (h  +  (∑_(k = 1)^khigh)(( 2  Sin[k h])/k)), 30]

khigh  =  100;
(1/π) (h  +  (∑_(k = 1)^khigh)(( 2  Sin[k h])/k))

khigh  =  500;
(1/π) (h  +  (∑_(k = 1)^khigh)(( 2  Sin[k h])/k))

khigh  =  5000;
N[(1/π) (h  +  (∑_(k = 1)^khigh)(( 2  Sin[k h])/k)), 30]

In fact when you run khigh all the way to infinity,you get:

(1/π) (h  +  (∑_(k = 1)^Infinity)(( 2  Sin[k h])/k))

When you integrate  (1/(2 π)) + (1/π) (∑_(k = 1)^khigh) Cos[
              k t] from a to b with 0 < a < b ≤ π, you get

     (∫_a^b) DiracApprox[khigh, t] dt =
                \(\((b\  - \ a)\)\/\(2\ π\) + \(1\/π\) \(∑\+\(k = 1\)\%khigh\( Sin[k\ b]\ \ \ \ \ \ \  - \ \ \ \ \ Sin[k\ a]\)\/\(\(\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \)\(\(k\)\(		\)\(.\)\)\)\)\)
                See how these turn out for random choices of a and b with  0 < a < b ≤ π,:

a  =  Random[Real, {0, Pi}];
b  =   Random[Real, {a, Pi}];
khigh  =  50;
N[(((b  -  a))/(2 π)) + (1/π) (∑_(k = 1)^khigh)(( Sin[k b] - Sin[k a])/k), 30]

\(
\(khigh\  = \ 500;\)\)
N[(((b  -  a))/(2 π)) + (1/π) (∑_(k = 1)^khigh)(( Sin[k b] - Sin[k a])/k), 30]

\(
\(khigh\  = \ 1000;\)\)
N[(((b  -  a))/(2 π)) + (1/π) (∑_(k = 1)^khigh)(( Sin[k b] - Sin[k a])/k), 30]

\(
\(khigh\  = \ Infinity;\)\)
N[(((b  -  a))/(2 π)) + (1/π) (∑_(k = 1)^khigh)(( Sin[k b] - Sin[k a])/k), 30]

When you integrate  (1/(2 π)) + (1/π) (∑_(k = 1)^khigh) Cos[
              k t] from a to b with  -π ≤ a < b < 0  you get

     (∫_a^b) DiracApprox[khigh, t] dt =
                \(\((b\  - \ a)\)\/\(2\ π\) + \(1\/π\) \(∑\+\(k = 1\)\%khigh\( Sin[k\ b] - Sin[k\ a]\)\/\(\(k\)\(		\)\(.\)\)\)\)
  See how these turn out for random choices of a and b with  -π ≤ a < b < 0,:

a  =  Random[Real, {-Pi, 0}];
b  =   Random[Real, {a, 0}];
khigh  =  50;
N[(((b  -  a))/(2 π)) + (1/π) (∑_(k = 1)^khigh)(( Sin[k b] - Sin[k a])/k), 30]

\(
\(khigh\  = \ 500;\)\)
N[(((b  -  a))/(2 π)) + (1/π) (∑_(k = 1)^khigh)(( Sin[k b] - Sin[k a])/k), 30]

\(
\(khigh\  = \ 1000;\)\)
N[(((b  -  a))/(2 π)) + (1/π) (∑_(k = 1)^khigh)(( Sin[k b] - Sin[k a])/k), 30]

\(
\(khigh\  = \ Infinity;\)\)
N[(((b  -  a))/(2 π)) + (1/π) (∑_(k = 1)^khigh)(( Sin[k b] - Sin[k a])/k), 30]

These experiments indicate that as khigh -> Infinity
            (∫_a^b) DiracApprox[khigh, t] dt -> (∫_a^b) DiracDelta[t] dt.
            
And if the functions
             DiracApprox[khigh,t] = (1/(2 π)) + (1/π) (∑_(k = 1)^khigh) Cos[
              k t]
converged to a normal function as khigh gets big, then DiracApprox[khigh,t] would converge to DiracDelta[t].
The trouble is that these DiracApprox[khigh, t]  do not converge to a normal function as khigh gets big. 

To get around this, the French mathematician Laurent Schwartz invented a new kind of convergence and a new kind of function so that it is now perfectly legitimate to say that for purposes of integration, DiracApprox[khigh,t]  does converge to DiracDelta[t].