Unprotect[en,fit,En,n,Ne,o,poly,scalar];
Needs["Utilities`CleanSlate`"]
CleanSlate[] ;
Remove["Global`*"];
Protect[en,fit,En, n,Ne,o,poly,scalar];

Section I. Representation theory

A test of the completeness relation: expand the energy eigenfunction of a harmonic oscillator in terms of the energy eigenfunctions of a particle in a box.

The energy of the particle in a box is an observable and therefore the pure states of the energy |(Subscript[e, n])⟩  satisfy the completeness relation

	(∑_(n=1)^∞)|(Subscript[e, n])⟩⟨(Subscript[e, n])| =Î						(1)

The particle in a box does not have a continuous spectrum. 
	Denote the pure states of the energy of a harmonic oscillator by |(Subscript[E, m])⟩.  Using (1) I can write

	|(Subscript[E, m])⟩ = Î|(Subscript[E, m])⟩ =(∑_(n=1)^∞)|(Subscript[e, n])⟩⟨(Subscript[e, n])|(Subscript[E, m])⟩  				(2)
	
In this exercise I want to test this equation. 

If you  paid attention you might have foreseen  some trouble brewing. The particle-in-the-box is confined inside a box of length L, whose infinitely repulsive walls are located at x=0 and x=L.  The harmonic oscillator vibrates around a position (Subscript[x, 0]).  This information is not contained in |(Subscript[e, n])⟩ and |(Subscript[E, m])⟩.  These two kets, |(Subscript[e, n])⟩ and |(Subscript[E, m])⟩, tell us what the energy of the state is but give no positional information. We don't know where the box is located and we do not know around which point the oscillator oscillates. It should be clear that this missing information is essential.  If the box is in London and the oscillator is in New York we can not express the wave function of one system as a linear combination of the wave function of the other system. To construct such a linear combination we must pay attention to these positions. How do we do that?

We know the pure states of the harmonic oscillator and of the particle in a box in the coordinate representation (Schrodinger representation). So let us act on Eq. 2 with the bra ⟨x| to obtain:

	⟨x|(Subscript[E, m])⟩ =(∑_(n=1)^∞)⟨x|(Subscript[e, n])⟩⟨(Subscript[e, n])|(Subscript[E, m])⟩  					(3)
	
Here ⟨x|(Subscript[e, n])⟩ is the energy wave function of the particle in a box, which is:

	⟨x|(Subscript[e, n])⟩ =(Sqrt[(2/L)]) Sin[((n π x)/L)]      n=1, 2, 3, .... x∈[0,L] 		(4)
	
	⟨x|(Subscript[e, n])⟩ = 0        x∉[0,L]						(5)
	
L is the length of the box and x is the position of the particle. This expression was obtained for a box whose walls are located at x=0 and x=L (see Metiu, Quantum Mechanics). If the box is positioned differently the wave functions are different. 

If we use Eqs. 4 and 5 in Eq. 3 we obtain

	⟨x|(Subscript[E, m])⟩ =(∑_(n=1)^∞)(Sqrt[(2/L)]) Sin[((n π x)/L)] ⟨(Subscript[e, n])|(Subscript[E, m])⟩        x∈[0,L]  		(6)
	
and 

	⟨x|(Subscript[E, m])⟩ = 0                        x ∉ [0, L]					(7)

Now the limitations of Eq. 2 are brought to the fore.  Eqs. 6 and 7, which result directly from Eq. 2 (by acting on it with the bra ⟨x| ) can only represent an oscillator whose motion is such that the position of the oscillating particle satisfies x∈[0,L]. There is a fake generality in Eq. 2. If we were to take it literally and assume it to be true under all circumstances we can make some tragic (well, as quantum mechanics goes) mistakes. We can try to represent an oscillator located in New York with functions that are zero outside London and therefore are equal to zero in New York. A bit of common sense and analysis needs to be injected when we use Eq. 2.  
	The solution to this “dilemma” is to make sure that the oscillator vibrates around the  point (Subscript[x, 0]) =L/2, located in the middle of the box, and that the box is bigger than the amplitude of the oscillations. We also must make sure that the box is not too big because constructing a function ⟨x|(Subscript[E, n])⟩ that is located in a narrow range of x, around (Subscript[x, 0]), as  a linear combination of function ⟨x|(Subscript[e, n])⟩    that are non-zero in a very large interval around (Subscript[x, 0]) is possible but is not easy. 

We can test Eq. 6 because we know that the harmonic oscillator wave function is

⟨x|(Subscript[E, n])⟩ =(Sqrt[(1/λ)]) (1/((π^(1/4))(Sqrt[(2^n)n!])))exp[-((((x-(Subscript[x, 0])))^2)/(2(λ^2)))] (Subscript[H, n])(((x-(Subscript[x, 0]))/λ))        n=0,1,2,...  	(8)

with 

	λ = (Sqrt[(ℏ/(μ ω))])							     	 (9)
and 

	(Subscript[x, 0]) = (L/2)									(10)

I took this from Metiu, Quantum Mechanics and I shifted the position around which the oscillator  vibrates to (Subscript[x, 0])=(L/2); this is why (x-(Subscript[x, 0])) appears in Eq. (8)  (Subscript[H, n])(x) is the Hermite polynomial of order n and λ is a length. 

The test of Eq. 6 proceeds as follows:

1. Calculate  ⟨(Subscript[e, n])|(Subscript[E, m])⟩ in
	 ⟨x|(Subscript[E, m])⟩ =(∑_(n=1)^∞)(Sqrt[(2/L)]) Sin[((n π x)/L)] ⟨(Subscript[e, n])|(Subscript[E, m])⟩        x∈[0,L]   	(11)
 by using 
	
		⟨(Subscript[e, n])|(Subscript[E, m])⟩ = (∫_(-∞)^∞)⟨(Subscript[e, n])|x⟩⟨x|(Subscript[E, m])⟩dx			(12)
		
2. Test whether the right hand side of Eq. 11 (which is the same as Eq. 6)  is equal to the left hand side,  by plotting them on the same graph. Note that if you want a very accurate representation you cannot test it by comparing the plots of the two functions because plots have low accuracy. For a convincing comparison you need to compare the numerical values. 
	
To perform these calculations we need to choose the parameters L, and λ carefully. In what follows I will let them be symbolic until I perform numerical calculations.

Here is the plan for the Mathematica program:

1. Create the function 

	en[x][n,L] 
	
    which gives  ⟨x|(Subscript[e, n])⟩
    
2. Create a function 

	En[n,λ,L][x]
	
    which gives ⟨x|(Subscript[E, n])⟩	
    
3. Calculate  ⟨(Subscript[e, n])|(Subscript[E, m])⟩ from Eq. 11:

		⟨(Subscript[e, n])|(Subscript[E, m])⟩ = (∫_(-∞)^∞)⟨(Subscript[e, n])|x⟩⟨x|(Subscript[E, m])⟩dx		 

4. Use these results to evaluate the right hand side of  

⟨x|(Subscript[E, m])⟩ =(∑_(n=1)^∞)(Sqrt[(2/L)]) Sin[((n π x)/L)] ⟨(Subscript[e, n])|(Subscript[E, m])⟩        x∈[0,L]  		 

5. Compare the right-hand-side of the above equation to the left-hand-side (namely ,⟨x|(Subscript[E, m])⟩), to see if they are equal. Study how this representation depends on the number of terms in the sum and on the values of L and λ.

Cell 1. The wave function for the particle in the box
⟨x|(Subscript[e, n])⟩ =(Sqrt[(2/L)]) Sin[((n π x)/L)]      n=1, 2, 3, .... x∈[0,L] 	(13)
	
	⟨x|(Subscript[e, n])⟩ = 0        x∉[0,L]				(14)

I will not use Eq. 14 in the definition but I will be careful to use the function only in the range x ∈ [0, L]

Unprotect[en];
Clear[en];
en[x_][n_,L_]:=Sqrt[2/L] Sin[ n Pi x/L]
Protect[en];

Check that the function works properly

en[x][n,L]

Cell 2. Define the function En[n, λ, L][x] which gives the harmonic oscillator wave function: 

⟨x|(Subscript[E, n])⟩ =(Sqrt[(1/λ)]) (1/((π^(1/4))(Sqrt[(2^n)n!])))exp[-((((x-(Subscript[x, 0])))^2)/(2(λ^2)))] (Subscript[H, n])(((x-(Subscript[x, 0]))/λ))     	(15)
In this equation 
	n = 0, 1, 2, ...
and 
	(Subscript[x, 0])=(L/2)						(16)

Unprotect[En];
Clear[En];
En[n_,λ_,L_][r_]:=(1/Sqrt[λ]) (1/Pi)^(1/4) (1/Sqrt[2^n n!]) Exp[-(1/2)(r-L/2)^2/λ^2] HermiteH[n,(r-L/2)/ λ]
Protect[En];

Check that it works properly

En[n,λ,L][x]

Note that the representation has nothing to do with the fact that the function we are trying to represent is the wave function of the energy ground state of an oscillator. We could have tried just as well to fit some arbitrary function.

Cell 3: Compare ⟨x|(Subscript[e, n])⟩  to ⟨x|(Subscript[E, m])⟩

I will call the set containing those values of the variable x for which the function f(x) differs from zero, the support of f. The support of ⟨x|(Subscript[e, n])⟩  is  x∈[0,L].  The support of ⟨x|(Subscript[E, 0])⟩ is controlled by the magnitude of the parameter λ (which appears at the exponent of ⟨x|(Subscript[E, n])⟩); this support is roughly x∈[(L/2)-3λ, ((L )/2) + 3λ].
We can test that this is true by making a plot. 

I take n=0, λ=0.42, L=0.5. These lengths are in Å and therefore x must be in Å also.

Plot the ground state of the harmonic oscillator for L=0.5 and λ=0.042

Plot[En[0,0.042,0.5][r] ,{r,0,0.5},PlotRange->All,AxesLabel->{"x in Å", "⟨x|(Subscript[E, 0])⟩"},
PlotStyle->Thickness[0.008]]

The harmonic oscillator wave function is centered around (Subscript[x, 0])=(L/2) =0.25  and the support is approximately  x∈[0.1, 0.4]  which is consistent with the rough estimate I gave above.

It is very important that the support of the function being represented ( i.e. ⟨x|(Subscript[E, n])⟩) is contained in the supports of the functions doing the representing (i.e.  ⟨x|(Subscript[e, n])⟩)

The values we have chosen for L and λ satisfy this condition.

Here is the plot for a few particle in a box wave functions for L=0.5 :

Plot[{en[x][1,0.5],en[x][2,0.5],en[x][3,0.5]},{x,0,0.5},PlotStyle->Thickness[0.008],AxesLabel->{"x A","⟨x|(Subscript[e, n])⟩ box"}, PlotLabel->"wave functions particle in box"]

Let us plot these three functions together with the ground state of the oscillator (which is blue)

Plot[{En[0,0.042,0.5][x],en[x][1,0.5],en[x][2,0.5],en[x][3,0.5]},{x,0,0.5},PlotStyle->Thickness[0.008]]

You see that the blue function (harmonic oscillator) is not at all similar to the three particle in a box wave functions. Nevertheless, the completeness theorem says that we can represent the oscillator function as a linear combination of particle in a box functions. We want to test whether this is true.  

Notice that the harmonic oscillator wave function is zero when x∈[0,0.1] and when x∈[0.4,0.5], while the particle in the box wave functions differ from zero in these regions. In the sum representing the oscillator wave function the particle in the box functions will have to cancel each other in x∈[0,0.1] and x∈[0.4,0.5] and to build up a peak around x=0.25. This seems like a tough task. Let us see if they can do it.

Cell 4: Calculate the coefficients ⟨(Subscript[e, n])|(Subscript[E, 0])⟩ by using	
	
⟨(Subscript[e, n])|(Subscript[E, 0])⟩ = (∫_(-∞)^∞)⟨(Subscript[e, n])|x⟩⟨x|(Subscript[E, 0])⟩dx				 (17)

Calculate the integral for arbitrary n, L and λ

int=Integrate[en[x][n,L]*En[0,λ,L][x],{x,0, L}]//FullSimplify

This looks very complicated, but it is manageable. The Erf function is calculated by Mathematica automatically when the argument is a number. The result of the integration seems to be a complex number but it is not. Here are a few values of the integral for n = 1, 2, 3... and L=0.5, λ=0.042. Note that they are all real numbers (because n is integer)

Table[int/.{L->0.5,λ->0.042},{n,1,10}]//Chop

Create a function 
coef[n][L,λ]
that evaluates the expression int we have obtained for the integral (which is ⟨(Subscript[e, n])|(Subscript[E, m])⟩ ). If we don't do that we would have to make a function coef which will do the integral for every value of n, which is possible but very inefficient.

coef[m_][L1_,λ1_]:=(((-(1/2))*I*(-1 + E^(I*n*Pi))*
   Sqrt[1/L]*Pi^(1/4)*
   Sqrt[λ]*
   (Erf[(L^2 - 2*I*n*Pi*
        λ^2)/(2*Sqrt[2]*L*
       λ)] + 
    Erf[(L^2 + 2*I*n*Pi*
        λ^2)/(2*Sqrt[2]*L*
       λ)]))/
  E^((n*Pi*(I*L^2 + 
      n*Pi*λ^2))/(2*L^2)))/.{λ->λ1,L->L1,n->m}

Check that it gives the same values for the coefficients as int

Table[{coef[n][0.5,0.042]//Chop, int/.{L->0.5,λ->0.042}}//Chop,{n,1,10}]

The function coef works correctly. Note that every other coefficient is zero. Could you have proven easily that this has to be true?

Cell 5. Create a function 

	fit[m][L,λ]

 that calculates

	(∑_(n=1)^m)(Sqrt[(2/L)]) Sin[((n π x)/L)] ⟨(Subscript[e, n])|(Subscript[E, 0])⟩  		(18)
	
which is the representation of the harmonic oscillator ground state.
Test the expansion.

Unprotect[fit];
Clear[fit]
fit[m_][L_,λ_]:=Module[{ co, funcList, sum},    
co=Table[coef[i][L,λ],{i,1,m}];
funcList=Table[en[x][n,L],{n,1,m}];
sum=co.funcList;
Plot[{sum,En[0,λ,L][x]},{x,0,L},PlotRange->All, PlotStyle-> {{Thickness[0.008],RGBColor[1,0,0]},{Thickness[0.008],Dashing[{0.04,0.04}],RGBColor[0,0,1]}}]
] 
Protect[fit];

I have designed the function fit[n][L,λ] to plot the sum representing the harmonic oscillator ground state wave function by a dashed, blue line and to plot the sum that is supposed to represent the harmonic wave function by a thick red line. If the fit is good that two plots should coincide.  The function fit is used in the function Manipulate, which provides two sliders that can be moved to change the number of terms in the sum and the size of the box. Move the slider and you can see the sum fitting the harmonic oscillator wave function change. For a box size of 0.5 Angstrom the sum with 13 terms gives a decent fit. 

Play with the sliders to get a feel for the role of the parameters L (as compared to λ, which is fixed). Essentially the larger the box the larger of number of terms in the sum. If the box is much larger than the support of the oscillator wave function, even 40 terms in the sum do not give an accurate representation.  There are no mechanical rules for getting the fit right other than using common sense and testing the results

Manipulate[ fit[n][L,0.042],{{n,2,"nbTerms"},2,40,1,Appearance->"Labeled"},{{L,0.5,"box size"},0.3,10,Appearance->"Labeled"}]

Section II:  The Gram-Schmidt orthogonalization method: how to construct good basis sets.

Cell 1: The use of Gram-Schmidt to orthogonalize a set of three vectors in (R^3)

Define three vectors, denoted o[1], o[2] and o[3] by using a random number generator to generate the coordinates

Unprotect[o];
Clear[o];
SeedRandom[2];
Table[o[i]=RandomReal[{-1,1},3],{i,1,3}]
Protect[o];

Test to see if they are orthonormal (it would be a miracle if they were)

Table[o[i].o[j],{i,1,3},{j,1,3}]

They are not. If they were orthogonal o[i].o[i] would be equal to 1 and o[i].o[j] would be zero when i ≠j

Start the Gram-Schmidt iteration scheme: 
	 |(Subscript[n, 1])⟩ =  |(Subscript[o, 1])⟩							(1)
	 |(Subscript[n, 2])⟩ =  |(Subscript[o, 2])⟩ - (( |(Subscript[n, 1])⟩)/(⟨(Subscript[n, 1])| (Subscript[n, 1])⟩))⟨(Subscript[n, 1])|(Subscript[o, 2])⟩				(2)
	 |(Subscript[n, 3])⟩ =  |(Subscript[o, 3])⟩ -(( |(Subscript[n, 1])⟩)/(⟨(Subscript[n, 1])| (Subscript[n, 1])⟩))⟨(Subscript[n, 1])|(Subscript[o, 3])⟩-(( |(Subscript[n, 2])⟩)/(⟨(Subscript[n, 2])| (Subscript[n, 2])⟩))⟨(Subscript[n, 2])|(Subscript[o, 3])⟩		(3)

Implement Eq. (1)

Unprotect[n];
Clear[n];
n[1]=o[1]

Implement Eq. (2)

n[2]=o[2]-n[1] (n[1].o[2]/(n[1].n[1]))

Test that n[1] is orthogonal to n[2]

n[1].n[2]

Implement Eq. (3)

n[3]=o[3]-n[1](n[1].o[3])/(n[1].n[1])-n[2](n[2].o[3])/(n[2].n[2])

Check that n[3] is orthogonal with n[2]

n[2].n[3]

Check the orthogonality of n[3] with n[1]

n[1].n[3]

Build a matrix with the elements n[i].n[j]

Table[n[i].n[j],{i,1,3},{j,1,3}]//Chop//MatrixForm

The vectors n[i], i = 1, 2, 3 are not normalized. Define the normalized vectors 

	ν[i] = n[i]/Sqrt[n[i].n[i]]

Unprotect[ν];
Clear[ν];
{ν[1],ν[2],ν[3]}=Table[n[i]/Sqrt[n[i].n[i]],{i,1,3}]

Test that they are orthogonal

Table[ν[i].ν[j],{i,1,3},{j,1,3}]//Chop//MatrixForm

Here is an interesting question: what happens within the Gram-Schmidt scheme if  o[3] is linearly dependent on o[1] and o[2]?

Let us consider perform the Gram-Schmidt orthonormalization in the case when 
	o[3] = a o[1] + b  o[2]

Generate o[1] and o[2]

Unprotect[o];
Clear[o]
SeedRandom[2];
Table[o[i]=RandomReal[{-1,1},3],{i,1,2}]

Generate o[3] = a o[1] + b o[2]

Clear[a,b]
o[3]=a*o[1] +b*o[2]

Now start the Gram-Schmidt iteration scheme:

n[1]=o[1]

n[2]=o[2]-n[1] (n[1].o[2]/(n[1].n[1]))

n[3]=o[3]-n[1](n[1].o[3])/(n[1].n[1])-n[2](n[2].o[3])/(n[2].n[2])//FullSimplify//Chop

n[3] generated by the Gram-Schmidt procedure is equal to zero if o[3] is a linear combination of the preceding vectors in the initial basis set (o[1] and o[2] in our case).  The Gram-Schmidt procedure automatically eliminates the term that is a linear combination of the others, from the orthonormal basis set being constructed. This property is general.

Cell 2: The application of the Gram-Schmidt procedure to create an orthonormal basis set polynomials.

In the general theory of the Gram-Schmidt procedure I have a set of old functions:

	⟨x|(Subscript[o, 1])⟩, ⟨x|(Subscript[o, 2])⟩, ⟨x|(Subscript[o, 3])⟩, ⟨x|(Subscript[o, 4])⟩

The orthogonal set constructed from this set is: 

	⟨x|(Subscript[n, 1])⟩ = ⟨x|(Subscript[o, 1])⟩							(1)
	⟨x|(Subscript[n, 2])⟩ = ⟨x|(Subscript[o, 2])⟩ - ((⟨x | (Subscript[n, 1])⟩)/(⟨(Subscript[n, 1]) |(Subscript[n, 1])⟩))⟨(Subscript[n, 1])|(Subscript[o, 2])⟩					(2)
	 ⟨x|(Subscript[n, 3])⟩ = ⟨x|(Subscript[o, 3])⟩ -((⟨x | (Subscript[n, 1])⟩)/(⟨(Subscript[n, 1])| (Subscript[n, 1])⟩))⟨(Subscript[n, 1])|(Subscript[o, 3])⟩-((⟨x | (Subscript[n, 2])⟩)/(⟨(Subscript[n, 2]) | (Subscript[n, 2])⟩))⟨(Subscript[n, 2])|(Subscript[o, 3])⟩		(3)
	.⟨x|(Subscript[n, 4])⟩ = ⟨x|(Subscript[o, 4])⟩ -((⟨x | (Subscript[n, 1])⟩)/(⟨(Subscript[n, 1])| (Subscript[n, 1])⟩))⟨(Subscript[n, 1])|(Subscript[o, 4])⟩-((⟨x | (Subscript[n, 2])⟩)/(⟨(Subscript[n, 2]) | (Subscript[n, 2])⟩))⟨(Subscript[n, 2])|(Subscript[o, 4])⟩ -((⟨x | (Subscript[n, 3])⟩)/(⟨(Subscript[n, 2]) | (Subscript[n, 3])⟩))⟨(Subscript[n, 3])|(Subscript[o, 4])⟩      (4)
	.
	⟨x|(Subscript[n, j])⟩  =  ⟨x|(Subscript[o, j])⟩ - (∑_(i=1 )^(j-1)){⟨x|(Subscript[n, i])⟩((⟨x| (Subscript[n, i])⟩⟨ (Subscript[n, i])| (Subscript[o, j])⟩)/(⟨(Subscript[n, i])| (Subscript[n, i])⟩))}  				     (5)
	
The functions ⟨ x|(Subscript[n, j])⟩ produced by this algorithm are not normalized; the functions

	⟨x|(Subscript[N, i])⟩ = ((⟨x|(Subscript[n, i])⟩)/(Sqrt[⟨(Subscript[n, i])| (Subscript[n, i])⟩]))							(5)
are. 

Apply this theory to the set 

	⟨x|(Subscript[o, n])⟩ = (x^n)  , n = 0, 1, 2, 3						(6)
	
The scalar product is defined by 

	⟨a|b⟩ ≡(∫_0^π)a(x)b(x)dx							(7)

The orthonormalized basis set ⟨x|(Subscript[N, i])⟩ will be used to represent functions of argument x, with x∈ [0, π]

To show how the Gram-Schmidt procedure works I will construct four orthogonal polynomials, starting with the "old basis set" {1, x, (x^2), (x^3)}.

Define the starting ("old") linearly independent set

Unprotect[o];
Clear[o]
o[1][x_]:=1
o[2][x_]:=x
o[3][x_]:=x^2
o[4][x_]:=x^3
Protect[o];

Define the scalar product 

	⟨a|b⟩ = (∫_0^π)a(x)b(x)dx

Unprotect[scalar];
Clear[scalar]
scalar[ a_,b_]:=Integrate[a*b,{x,0,Pi}]
Protect[scalar];

Check that the function scalar works

scalar[o[1][x],o[3][x]]

The matrix having the elements 

	(Subscript[S, i,j]) = ,⟨(Subscript[o, i])|(Subscript[o, j])⟩   
	
is called the overlap matrix. Obviously, if the functions {(Subscript[o, i])(x)} are orthogonal the overlap matrix is diagonal. If the functions were orthonormal the overlap matrix would be equal to the unit matrix. Calculate the overlap matrix:

S=Table[scalar[o[j][x],o[i][x]],{i,1,4},{j,1,4}]//N

Check that the functions ⟨x|(Subscript[o, i])⟩  are linearly independent (the determinant of the overlap matrix is not close to zero)

Det[S]

Use the Gram-Schmidt algorithm to define new functions, which are orthogonal. The first new function is 
⟨x|(Subscript[n, 1])⟩ = ⟨x|(Subscript[o, 1])⟩
and it is given by n[1][x]:

Unprotect[n];
Clear[n]
n[1][x_]=o[1][x]
Protect[n];

Define the second  new function is
	⟨x|(Subscript[n, 2])⟩ = ⟨x|(Subscript[o, 2])⟩ - ((⟨x | (Subscript[n, 1])⟩)/(⟨(Subscript[n, 1]) |(Subscript[n, 1])⟩))⟨(Subscript[n, 1])|(Subscript[o, 2])⟩
and it is given by n[2][x]:

Unprotect[n];
n[2][x_]=o[2][x]-n[1][x]*scalar[n[1][x],o[2][x] ]/scalar[n[1][x],n[1][x]] 
Protect[n];

The third new function is 
	⟨x|(Subscript[n, 3])⟩ = ⟨x|(Subscript[o, 3])⟩ -((⟨x | (Subscript[n, 1])⟩)/(⟨(Subscript[n, 1])| (Subscript[n, 1])⟩))⟨(Subscript[n, 1])|(Subscript[o, 3])⟩-((⟨x | (Subscript[n, 2])⟩)/(⟨(Subscript[n, 2]) | (Subscript[n, 2])⟩))⟨(Subscript[n, 2])|(Subscript[o, 3])⟩
and it is given by n[3][x]:

Unprotect[n];
n[3][x_]=o[3][x]-n[1][x]*scalar[n[1][x],o[3][x] ]/scalar[n[1][x],n[1][x]] -n[2][x]*scalar[n[2][x],o[3 ][x]]/scalar[n[2][x],n[2][x]] 
Protect[n];

Define the fourth new function 
⟨x|(Subscript[n, 4])⟩ = ⟨x|(Subscript[o, 4])⟩ -((⟨x | (Subscript[n, 1])⟩)/(⟨(Subscript[n, 1])| (Subscript[n, 1])⟩))⟨(Subscript[n, 1])|(Subscript[o, 4])⟩-((⟨x | (Subscript[n, 2])⟩)/(⟨(Subscript[n, 2]) | (Subscript[n, 2])⟩))⟨(Subscript[n, 2])|(Subscript[o, 4])⟩ -((⟨x | (Subscript[n, 3])⟩)/(⟨(Subscript[n, 2]) | (Subscript[n, 3])⟩))⟨(Subscript[n, 3])|(Subscript[o, 4])⟩  
is given by n[4][x]:

Unprotect[n];
n[4][x]=o[4][x]-n[1][x]*scalar[n[1][x],o[4 ][x]]/scalar[n[1][x],n[1] [x]]-n[2][x]*scalar[n[2][x],o[4 ][x]]/scalar[n[2][x],n[2][x]] -n[3][x]*scalar[n[3][x],o[4][x] ]/scalar[n[3][x],n[3][x]] 
Protect[n];

Simplify them and  display them:

Table[n[i][x]//Expand,{i,1,4}]//TableForm

Calculate the overlap matrix for the new functions:

Table[scalar[n[i][x],n[j][x]],{i,1,4},{j,1,4}]//MatrixForm

The functions n[i][x] are orthogonal. Normalize them:

	(Subscript[N, i])(x) = (Subscript[n, i])(x)/(Sqrt[⟨(Subscript[n, i])|(Subscript[n, i ⟩])])

Unprotect[Ne];
Clear[Ne]
{Ne[1],Ne[2],Ne[3],Ne[4]}=Table[n[i][x]/Sqrt[scalar[n[i][x],n[i][x]]],{i,1,4}]
Protect[Ne];

Display them:

Table[Ne[i],{i,1,4}]//Expand//TableForm

Calculate the overlap matrix

Table[Integrate[Ne[i] Ne[j],{x,0,Pi}],{i,1,4},{j,1,4}]//MatrixForm

This example shows how the Gram-Schmidt procedure is implemented. It clearly does what is supposed to.

Cell 3. Test how well a basis set of orthogonal polynomials fits the function (x^2)sin(x) in the interval x∈[0, π].

Because writing a program to implement Gram-Schmidt for a large basis set requires a lot of typing, I will use the Mathematica function Orthogonalize to create a basis set consisting of orthogonal, real (as opposed to complex) polynomials (Subscript[ν, 0])(x), (Subscript[ν, 1])(x), ...  (Subscript[ν, n])(x),  defined in the interval  x∈[0, π].  The scalar product is defined as

	⟨a|b⟩  ≡ (∫_0^π)a(x) b(x)dx					(1)
		
I want to test how well these polynomials represent a function f(x),  x∈[0, π].  The general equation is 

	f(x) =(∑_(i=0)^n)(Subscript[a, i]) (Subscript[ν, i])(x)						(2)
with 
	 (Subscript[a, j]) = ⟨(Subscript[ν, j])|f⟩ = (∫_0^π)(Subscript[ν, j])(x) f(x)dx				(3)

8.1  The Mathematica function Orthogonalize

The function Orthogonalize takes as arguments a list of expressions and the definition of the scalar product. It returns a list of expressions which are orthonormal with respect to the scalar product given in the argument of Orthogonalize. To see how it works I create a list of three three-dimensional vectors having components chosen at random. The scalar product is the function Dot which gives the dot product of the vectors:

Dot[{a,b,c},{A, B, C}]

Create three vectors by giving them components chosen at random.

SeedRandom[5]
ta=Table[{RandomReal[],RandomReal[],RandomReal[]},{i,1,3}]

Check that they are linearly independent by calculating the overlap matrix and its determinant

S=Table[Dot[ta[[i]],ta[[j]]],{i,1,3},{j,1,3}]

Det[S]

The three vectors are linearly independent. Use Orthogonalize to generate three orthonormal vectors

ve=Orthogonalize[ta,Dot]

Test whether these vectors are orthonormal:

Table[Dot[ve[[i]],ve[[j]] ],{i,1,3},{j,1,3}]//Chop//MatrixForm

They are orthonormal. The function Orthogonalize gave the correct result.

8.2  Use the function Orthogonalize to create a set of orthonormal polynomials from the set 

	1, x, (x^2), (x^3), (x^4)

with respect to the scalar product  

	⟨a|b⟩ = (∫_0^π)a(x)b(x)dx				(1)

νν=Orthogonalize[{1,x,x^2,x^3, x^4},Integrate[#1 *#2,{x,0,Pi}]&]//Expand

The function 
	Integrate[#1*#2,{x,0,π}]&  
used as an argument in Orthogonalize defines the scalar product with respect to which the orthonormalization is performed.

Test that the polynomials are orthonormal.

Table[Integrate[νν[[i]]*νν[[j]],{x,0,Pi}],{i,1,5},{j,1,5}]//MatrixForm

They are orthonormal.

Define the function poly[n,x] which generates a set of orthogonal polynomials from {1, x, (x^2), ..., (x^n)}. The scalar product is 

	⟨a|b⟩ = (∫_0^π)a(x)b(x)dx

Unprotect[poly];
Clear[poly]
poly[n_Integer,x_]:=Module[{pElem,a,tb},
pElem=Table[x^i,{i,0,n}];
a=Orthogonalize[pElem,Integrate[#1*#2,{x,0,Pi}]&];
tb =Table[polynomial[i-1]=a[[i]],{i,1,n }];
tb//N//Expand
    ]
Protect[poly];

Test that the function executes

poly[5,x]//MatrixForm

The result is the same as the one obtained when I did the calculation using Gram Schmidt one step at a time.

8.3  Use the orthogonal polynomials (Subscript[ν, n])(x) as a basis set to expand (x^2) Sin(x).  We write 

⟨x|f⟩ ≡ (x^2)sin(x) = (∑_(n=0)^4)(Subscript[a, n])⟨x|(Subscript[ν, n])⟩ =  (∑_(n=0)^4)(Subscript[a, n])(Subscript[ν, n])(x)		(2)
	
with

(Subscript[a, n]) =⟨(Subscript[ν, n])|f⟩ =(∫_0^π)⟨(Subscript[ν, n])|x⟩⟨x|f⟩dx=

= (∫_0^π)(Subscript[ν, n])(x)(x^2)sin(x) dx          				(3)

To begin with, use as a basis set polynomials up to the fourth order. Use poly to generate the list po containing orthonormal polynomials up to fourth order and evaluate the right hand side of 

	 (x^2)sin(x) = (∑_(n=0)^4)(Subscript[a, n])(Subscript[ν, n])(x)

The first four coefficients (Subscript[a, n])  are (use Eq. (3))

Integrate[#*x^2 Sin[x],{x,0,Pi}]&/@poly[4,x]

m=4; Clear[a];
ν=poly[m,x];
Table[a[i]=Integrate[x^2Sin[x]*ν[[i]],{x,0,Pi}]//N//Chop[#,10^-6]&,{i,0,m}];
rhs=Sum[a[i]*ν[[i]],{i,1,m}]//Expand//Chop[#,10^-6]&

The coefficients a[i] are

Table[a[i],{i,1,4}]

Test how good the fit is by plotting (x^2)sin(x)  (in red) and (∑_(n=0)^4)(Subscript[c, n])(Subscript[p, n])(x) (dashed black) in the same graph.

pl3=Plot[{x^2 Sin[x],rhs},{x,0,Pi},PlotStyle-> {{Thickness[0.006],RGBColor[1,0,0]},{Dashing[{0.04,0.04}],Thickness[0.01],RGBColor[0,0,0]}} ]

The fit is not perfect because we used only a few polynomials. Let's try 10.

m=10;
po=poly[m,x];
Table[c[i]=Integrate[x^2Sin[x]*po[[i]],{x,0,Pi}]//N//Chop[#,10^-6]&,{i,0,m}];
rhs=Sum[c[i]*po[[i]],{i,1,m}]//Expand//Chop[#,10^-6]&

Plot the function and the representation. As you can see below the fit is rather good.

pl10=Plot[{x^2 Sin[x],rhs},{x,0,Pi},PlotStyle-> {{Thickness[0.006],RGBColor[1,0,0]},{Dashing[{0.04,0.04}],Thickness[0.01],RGBColor[0,0,0]}} ]

8.4  Isn't this just a power series expansion?  Almost, but not quite. Let us compare the representation we got by using orthonormal polynomials, to a power series representation of x^2 sin(x).

powerS=Series[x^2Sin[x],{x,0,4}]//Normal

Compare this to the expansion of the same order in orthonormal polynomials

m=4;
po=poly[m,x];
Table[c[i]=Integrate[x^2Sin[x]*po[[i]],{x,0,Pi}]//N//Chop[#,10^-6]&,{i,0,m}];
rhs=Sum[c[i]*po[[i]],{i,1,m}]//Expand//Chop[#,10^-6]&

Plot the power series expansion (red), the representation as a sum of orthonormal polynomials (dashed black line) and (x^2)sin(x)  (thin, blue line)

pl0=Plot[{powerS,rhs,x^2Sin[x]},{x,0,Pi},PlotStyle-> {{Thickness[0.008],RGBColor[1,0,0]},{Dashing[{0.04,0.04}],Thickness[0.008],RGBColor[0,0,0]},{RGBColor[0,0,1],Thickness[0.005]} } ]

The power series (red line) is good near the origin because we expanded around the point x=0. The polynomial expansion is good on the whole interval. The thin line is (x^2)sin(x) and the dashed black line is the polynomial fit.

Expanding in power series around x=0 was not too smart. What would happen if we expand around x=π/2 ? How good is the power series then?

The power series expansion is middleS:

middleS= Series[x^2Sin[x],{x,Pi/2,3}]  //Normal//Expand//N

plp=Plot[{middleS ,rhs,x^2Sin[x]},{x,0,Pi},PlotStyle-> {{Thickness[0.008],RGBColor[0.5,0,0.5]},{Dashing[{0.04,0.04}],Thickness[0.008],RGBColor[1,0,0]},{RGBColor[0,0,1],Thickness[0.005]} } ]

(x^2)sin(x)is in blue, the polynomial expansion is the dashed, red line and the power series expansion is in red. A power series expansion around x=π/2  is better than expanding around x=0, but not as good as the polynomial expansion. This is very interesting because both expansions are polynomials of the same order. This behavior is observed because the power series is designed to give an accurate representation around the expansion point, while the expansion in orthogonal polynomials is designed to give a good fit everywhere.

Cell 9: Test linear independence and show how one can fix a bad basis set.

I tested whether the set

	{sin(x), 1, x, (x^2), (x^3),(x^4),(x^5)}				(1)
	
is linearly independent and I found that it is nearly dependent (the determinant of the overlap matrix is very small).  Usually this gives trouble in numerical work.   However, if I use Gram-Schmidt to convert this basis set into an orthonormal set I get a very good basis set.

The function basis[n] defined below creates the set shown in Eq. 1.

Clear[basis]
basis[n_]:= Module[{} ,
tb=Table[ x^i,{i,0,n}];
{Sin[x],tb}//Flatten]

Look at the result:

b=basis[5]

The scalar product is 

	⟨a|b⟩  ≡ (∫_0^π)a(x)b(x)dx

Calculate the overlap matrix
	
	(Subscript[S, ij])=⟨(Subscript[b, i])|(Subscript[b, j])⟩

S=Table[Integrate[b[[i]]*b[[j]],{x,0,Pi}],{i,1,Length[b]},{j,1, Length[b]}]//N

Calculate the determinant of S

Det[S]

This is an anxiety-causing, low value. Not linear dependence, but too close for comfort.  I suspected that the use of sin(x) in the set is responsible for the small value of the determinant. I removed it to test the remaining set for linear independence

rb=Table[x^i,{i,0,5}]

S=Table[Integrate[b[[i]]*b[[j]],{x,0,Pi}],{i,1,Length[b]},{j,1, Length[b]}]//N

Det[S]

The determinant is very small. It appears that the set is nearly linearly dependent.  Let us construct an orthonormal set out of the initial set that contained sin(x)

{ba[1],ba[2],ba[3],ba[4],ba[5],ba[6],ba[7]}=Orthogonalize[b,Integrate[#1 #2, {x,0,Pi}]&]//N//FullSimplify//Chop//Expand

Display the new basis set more clearly by using MatrixForm

{ba[1],ba[2],ba[3],ba[4],ba[5],ba[6],ba[7]}//MatrixForm

Look at the overlap matrix for this orthonormal basis set

S=Table[Integrate[ba[i]*ba[j],{x,0,Pi} ],{i,1,7},{j,1,7}]//N//Chop[#,10^-6]&//MatrixForm

This is extremely interesting and shows the tremendous utility of the Gram-Schmidt procedure. We started with a miserable basis set that was almost linearly dependent and ended up with an orthonormal basis set of good quality.