(*^
::[	Information =

	"This is a Mathematica Notebook file.  It contains ASCII text, and can be
	transferred by email, ftp, or other text-file transfer utility.  It should
	be read or edited using a copy of Mathematica or MathReader.  If you 
	received this as email, use your mail application or copy/paste to save 
	everything from the line containing (*^ down to the line containing ^*)
	into a plain text file.  On some systems you may have to give the file a 
	name ending with ".ma" to allow Mathematica to recognize it as a Notebook.
	The line below identifies what version of Mathematica created this file,
	but it can be opened using any other version as well.";

	FrontEndVersion = "Macintosh Mathematica Notebook Front End Version 2.2";

	MacintoshStandardFontEncoding; 
	
	fontset = title, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeTitle, center, M7, bold, e8,  24, "Times"; 
	fontset = subtitle, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeTitle, center, M7, bold, e6,  18, "Times"; 
	fontset = subsubtitle, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeTitle, center, M7, italic, e6,  14, "Times"; 
	fontset = section, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeSection, grayBox, M22, bold, a20,  18, "Times"; 
	fontset = subsection, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeSection, blackBox, M19, bold, a15,  14, "Times"; 
	fontset = subsubsection, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeSection, whiteBox, M18, bold, a12,  12, "Times"; 
	fontset = text, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7,  14, "Times"; 
	fontset = smalltext, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7,  10, "Times"; 
	fontset = input, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeInput, M42, N23, bold, L-5,  12, "Courier"; 
	fontset = output, output, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, L-5,  12, "Courier"; 
	fontset = message, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, R65535, L-5,  12, "Courier"; 
	fontset = print, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, L-5,  12, "Courier"; 
	fontset = info, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, B65535, L-5,  12, "Courier"; 
	fontset = postscript, PostScript, formatAsPostScript, output, inactive, noPageBreakInGroup, nowordwrap, groupLikeGraphics, M7, l34, w405, h197,  12, "Courier"; 
	fontset = name, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, italic,  10, "Geneva"; 
	fontset = header, inactive, noKeepOnOnePage, preserveAspect, M7,  12, "Times"; 
	fontset = leftheader, inactive, L2,  12, "Times"; 
	fontset = footer, inactive, noKeepOnOnePage, preserveAspect, center, M7,  12, "Times"; 
	fontset = leftfooter, inactive, L2,  12, "Times"; 
	fontset = help, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7,  10, "Times"; 
	fontset = clipboard, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7,  12, "Times"; 
	fontset = completions, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7,  12, "Times"; 
	fontset = special1, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7,  12, "Times"; 
	fontset = special2, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7,  12, "Times"; 
	fontset = special3, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7,  12, "Times"; 
	fontset = special4, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7,  12, "Times"; 
	fontset = special5, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7,  12, "Times"; 
	paletteColors = 128; currentKernel; 
]
:[font = section; inactive; Cclosed; preserveAspect; cellOutline; backColorRed = 58981; backColorGreen = 58981; backColorBlue = 58981; startGroup]
Example  1:   Two Coupled Harmonic Oscillators, 1 spring
:[font = input; preserveAspect]
Clear["Global`*"];
:[font = subsubsection; inactive; preserveAspect; startGroup]
Solution
:[font = text; inactive; preserveAspect; startGroup]
The kinetic energy and potential energy for the coupled oscillators are:
:[font = input; preserveAspect]
T= (m/2) ( (x[1])'[t]^2 + (x[2])'[t]^2 ); 
 
V=( (k1/2)( x[2][t] - x[1][t])^2 
   + k2/2   x[1][t]^2 
   + k2/2   x[2][t]^2
  ) /.{k2->0} 
:[font = text; inactive; preserveAspect]
where the coordinates and momenta
:[font = input; preserveAspect]
q={x[1][t], x[2][t]};
p=D[q,t]
:[font = text; inactive; preserveAspect]
are the displacements from the equilibrium configuration.  

Note here, the use of the x[i][t] notation. This allows us to use the i index to label the coordinates and momenta, while still taking advantage of the derivative shorthand notation (x[i])'[t].

The Tij and Vij matrices for small oscillations follow from the coefficients in front of the quadratic terms in x'[i][t] and x[i][t],
:[font = input; preserveAspect]
Tij=Outer[D[T,#1,#2]&,p,p];
Vij=Outer[D[V,#1,#2]&,q,q]; 
:[font = text; inactive; preserveAspect]
The eigenvector matrix  is
:[font = input; preserveAspect]
matrix = Vij - w2 Tij;
matrix //MatrixForm
:[font = text; inactive; preserveAspect]
where w2 is the square of the eigenfrequency. The eigenfrequencies follow from the setting the determinate of matrix equal to zero and solving for w2, 
:[font = input; preserveAspect]
w2Sol= Solve[ Det[matrix]==0, w2] //Simplify
:[font = text; inactive; preserveAspect]
The square root of w2 gives the two eigenfrequencies.The eigenvector {a1,a2} for  the first eigenfrequency  follows from solving  the  algebraic equations (M_{ij} v_{j}=0):
:[font = input; preserveAspect]
vector= Array[a,Length[q]]
:[font = input; preserveAspect]
eq= Thread[0==matrix.vector]
:[font = text; inactive; preserveAspect]
To obtain the appropriate eigenvector, we must substitute in the appropriate eigenvalue for w2. Here, we solve for the first eigenvector by substituting with w2Sol[[1]]. 
:[font = input; preserveAspect]
eq1 = eq /.w2Sol[[1]] //ExpandAll //Simplify
:[font = text; inactive; preserveAspect]
The normalization equation is:
:[font = input; preserveAspect]
eqNorm={vector.vector==1};
:[font = text; inactive; preserveAspect]
Solving eq1 for the eigenvector (i.e., the components of vector), we have:
:[font = input; preserveAspect]
ev1rule=Solve[Join[eq1,eqNorm],vector][[2]] 
:[font = text; inactive; preserveAspect]
The eigenvector for the second normal frequency follows by substituting with w2Sol[[2]]. 
:[font = input; preserveAspect]
eq2 = eq /.w2Sol[[2]] //ExpandAll //Simplify
:[font = input; preserveAspect; endGroup]
ev2rule=Solve[Join[eq2,eqNorm],vector][[2]] 
:[font = text; inactive; preserveAspect; startGroup]
The general solution is a linear combination of terms of the form  
Cos[w1 t + phase1]  and   Cos[w2 t + phase2]   where w1 and w2 are the two eigenfrequencies.  
Here, we will identify w1 with freq1, and w2 with freq2. 
If we define the eigenmatrix: 
:[font = input; preserveAspect]
 eigenmatrix= vector //.{ev1rule,ev2rule}
:[font = text; inactive; preserveAspect]
 the  general solution is  
:[font = input; preserveAspect]
xRule={x1[t],x2[t]} ->
       ( Transpose[eigenmatrix]. {c1 Cos[ freq1 t + phase1] 
                                 ,c2 Cos[ freq2 t + phase2]} 
        )//Thread //Simplify
:[font = text; inactive; preserveAspect]
Note, we have a special case for a zero frequency as this corresponds to pure translation. 
:[font = input; preserveAspect; startGroup]
xRule={x1[t],x2[t]} ->
       ( Transpose[eigenmatrix]. {c1  t + phase1
                                 ,c2 Cos[ freq2 t + phase2]} 
        )//Thread //Simplify
:[font = text; inactive; preserveAspect; endGroup]
where c1 and c2 are the relative amplitudes of the two normal modes.  The eigenfrequencies,  freq1 and freq2,  are  
:[font = input; preserveAspect]
freqrule= Thread[{freq1,freq2} -> (Sqrt[w2] /.w2Sol) ]
:[font = text; inactive; preserveAspect]
 The two normal modes correspond to the special cases c2=0 and  c1=0.
 We consider the value of the parameters given by 
:[font = input; preserveAspect]
 values= {k1->1, k2->1, m->1, phase1->0, phase2->0};
:[font = text; inactive; preserveAspect]
For simplicity, we collect all the appropriate rules into allRules
:[font = input; preserveAspect]
allRules=Join[xRule,freqrule,values];
:[font = text; inactive; preserveAspect]
Now, we examine the first normal mode {c1=1, c2=0}. We use Plot to  plot the motion of the two coordinates {x1[t],x2[t]} on the same graph separated by a distance of +/-1. We use the command ParametricPlot  to map the trajectory of this system in the  {x1[t],x2[t]} plane. 
 (We write out each plot command separately here. In the next example, we will define a function to do this step.)
:[font = input; preserveAspect]
plot1a= 
Plot[{x1[t]+1, x2[t]-1} //.allRules //.{c1->1, c2->0} //Evaluate
       ,{t,0,20}
       ,PlotStyle->{{Thickness[0.010]},{Thickness[0.005]}}
    ];
:[font = input; preserveAspect]
plot1b= ParametricPlot[ 
  {x1[t], x2[t]} //.allRules //.{c1->1, c2->0} //Evaluate
 ,{t,0,20}];
:[font = input; preserveAspect]
Show[GraphicsArray[{plot1a,plot1b}]];
:[font = text; inactive; preserveAspect]
The thick line corresponds to x[1], and the thiner line to x[2]. From both plots, we see that the motion is simple translation. 
:[font = text; inactive; preserveAspect]
Next, we examine the second normal mode {c1=0, c2=1}. 
:[font = input; preserveAspect]
plot2a= 
Plot[{x1[t]+1, x2[t]-1} //.allRules //.{c1->0, c2->1} //Evaluate
       ,{t,0,20}
       ,PlotStyle->{{Thickness[0.010]},{Thickness[0.005]}}
    ];
:[font = input; preserveAspect]
plot2b= ParametricPlot[ 
  {x1[t], x2[t]} //.allRules //.{c1->0, c2->1} //Evaluate
 ,{t,0,20}];
:[font = input; preserveAspect]
Show[GraphicsArray[{plot2a,plot2b}]];
:[font = text; inactive; preserveAspect; plain; italic]
Comparing the second normal mode to the first, we note that this has a lower eigenfrequency. Also, we note that the object are in phase.
:[font = text; inactive; preserveAspect]
Finally, we examine a mixed mode {c1=1, c2=1}. 
:[font = input; preserveAspect]
plot3a= 
Plot[{x1[t]+1, x2[t]-1} //.allRules //.{c1->1, c2->1} //Evaluate
       ,{t,0,20}
       ,PlotStyle->{{Thickness[0.010]},{Thickness[0.005]}}
    ];
:[font = input; preserveAspect]
plot3b= ParametricPlot[ 
  {x1[t], x2[t]} //.allRules //.{c1->1, c2->1} //Evaluate
 ,{t,0,20}];
:[font = input; preserveAspect; endGroup; endGroup; endGroup]
Show[GraphicsArray[{plot3a,plot3b}]];
:[font = section; inactive; preserveAspect; cellOutline; backColorRed = 58981; backColorGreen = 58981; backColorBlue = 58981; startGroup]
Problem  2:   Two Coupled Harmonic Oscillators, 3 springs
:[font = text; inactive; preserveAspect; fontColorRed = 65535; endGroup]
Copy and modify from above. 
:[font = section; inactive; preserveAspect; cellOutline; backColorRed = 58981; backColorGreen = 58981; backColorBlue = 58981; startGroup]
Problem  3:   Three Coupled Harmonic Oscillators, 2 springs
:[font = text; inactive; preserveAspect; fontColorRed = 65535; endGroup]
Copy and modify from above. 
^*)