(*^
::[	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, preserveAspect, groupLikeGraphics, M7, l34, w195, h198,  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; startGroup]
Example 1:  Un-Damped Linear Oscillator 
:[font = input; preserveAspect]
Clear["Global`*"];
:[font = subsubsection; inactive; preserveAspect]
Required packages
:[font = input; preserveAspect]
Needs["Algebra`Trigonometry`"];
:[font = subsection; inactive; preserveAspect; cellOutline; startGroup]
Analytic Solution:  
:[font = text; inactive; preserveAspect]
Let's start with the basic equation for SHO. 
Since this is the un-damped case, I'll set gam to zero. 
:[font = input; preserveAspect]
eq1= x''[t]+ 2 gam x'[t] + w0^2 x[t] ==0 /.{gam->0}
:[font = text; inactive; preserveAspect]
We use DSolve to obtain the analytical solution. 
:[font = input; preserveAspect]
dsol= DSolve[{eq1, x[0]==x0,x'[0]==v0},x[t],t][[1]] //Simplify //PowerExpand
:[font = text; inactive; preserveAspect]
Choosing some values we define x[t,w0] where w0 is the resonant frequency. 
:[font = input; preserveAspect]
values={x0->1, v0->0, gam->1};
x[t_,w0_]= x[t] /.dsol /.values
:[font = text; inactive; preserveAspect]
Now we plot a case for w0 < 1.
:[font = input; preserveAspect]
Plot[ x[t,0.7] ,{t,0,10}];
:[font = text; inactive; preserveAspect]
Now we plot a case for w0 = 1.
In the un-damped case, this limiting procedure is unnecessary, but I will include it here since we will need it for the general case. 
:[font = input; preserveAspect]
limit =Limit[ x[t,w0] ,w0->1]
:[font = input; preserveAspect]
(* Patch to take care of limit: *)
x[t_,w0_]:= limit /; 0.99 < w0 < 1.01
:[font = input; preserveAspect]
Plot[ x[t,1] ,{t,0,10}];
:[font = text; inactive; preserveAspect]
Now we plot a case for w0 > 1.
:[font = input; preserveAspect]
Plot[ x[t,10] ,{t,0,10},PlotRange->All];
:[font = text; inactive; preserveAspect]
Now we overlay the plots. 
:[font = input; preserveAspect]
Plot[ {x[t,10] 
      ,x[t,1]
      ,x[t,0.7]
      }
      ,{t,0,10}
      ,PlotRange->All
      ,PlotStyle->{{Thickness[0.004],RGBColor[1,0,0]}
                  ,{Thickness[0.008],RGBColor[0,1,0]}
                  ,{Thickness[0.012],RGBColor[0,0,1]}
                  }
      ];
:[font = text; inactive; preserveAspect]
We can also make a 3-D plot in {t,w0} space. 
:[font = input; preserveAspect]
Plot3D[ x[t,w0] 
      ,{t,0,2}
      ,{w0,0.1,20}
      ,PlotPoints->50
      ,PlotRange->All
      ,ViewPoint->{1.5,4,2}
      ];
:[font = text; inactive; preserveAspect]
We now define {q1.p1} so we can generate a phase-plot. 
:[font = input; preserveAspect]
q1[t_]= x[t,20] //Simplify
p1[t_]= q1'[t]  //Simplify
:[font = text; inactive; preserveAspect]
Here is a 2-D phase-plot. 
:[font = input; preserveAspect]
ParametricPlot[{q1[t],p1[t]} ,{t,0.01,2},PlotRange->All];
:[font = text; inactive; preserveAspect]
Here is a 3-D phase-plot. This provides us additional information about the motion. 
:[font = input; preserveAspect; endGroup]
ParametricPlot3D[{q1[t],p1[t],t},{t,0.01,2}
  ,BoxRatios->{1,1,1}
  ,PlotPoints->100
  ,PlotRange->All
  ];
:[font = subsection; inactive; preserveAspect; cellOutline; startGroup]
Numeric Solution:  
:[font = text; inactive; preserveAspect]
We substitute values into the eqs to be solved by NDSolve. 
:[font = input; preserveAspect]
eqs= {eq1, x[0]==x0,x'[0]==v0} /.values 
:[font = text; inactive; preserveAspect]
We choose w0=20, and generate a numeric solution. 
:[font = input; preserveAspect]
ndsol= NDSolve[eqs /.{w0->20}  ,x[t],{t,0,2}][[1]]
:[font = text; inactive; preserveAspect]
We assign values to {q2,p2}.  Note how we can easily take the derivative of a numeric function. 
:[font = input; preserveAspect]
q2[t_] = x[t] /.ndsol;
p2[t_] = q2'[t];
{q2[0.1],p2[0.1]}
:[font = text; inactive; preserveAspect]
Here is a 2-D phase-plot. This matches the above.
:[font = input; preserveAspect]
ParametricPlot[{q2[t],p2[t]},{t,0.01,2},PlotRange->All];
:[font = text; inactive; preserveAspect]
Here is a 3-D phase-plot. This matches the above. This provides us additional information about the motion. 
:[font = input; preserveAspect; endGroup; endGroup]
ParametricPlot3D[{q2[t],p2[t],t},{t,0.01,2}
  ,BoxRatios->{1,1,1}
  ,PlotPoints->100
  ,PlotRange->All
  ];
:[font = section; inactive; preserveAspect; cellOutline; startGroup]
Problem 2:  Damped Linear Oscillator 
:[font = text; inactive; preserveAspect; fontColorRed = 65535; endGroup]
Copy and modify from above. 
:[font = section; inactive; preserveAspect; cellOutline; fontColorBlue = 65535]
Optional Problems:
1) Consider a Sin[w t] Driving Force
2) Consider a Sin[w t] Driving Force with damping
3) Consider an additional b x^3 term and w0-> 2 I (Duffing's equation).
4) Consider the above with an additional Driving Force (Forced Duffing's equation).
;[s]
2:0,0;19,1;262,-1;
2:1,19,14,Times,1,18,0,0,65535;1,16,12,Times,1,14,0,0,65535;
:[font = section; inactive; Cclosed; preserveAspect; cellOutline; startGroup]
Example 3a:  Driving Force:
:[font = input; preserveAspect]
Clear["Global`*"];
:[font = subsubsection; inactive; preserveAspect]
Required packages
:[font = input; preserveAspect]
Needs["Algebra`Trigonometry`"];
:[font = subsection; inactive; preserveAspect; cellOutline; startGroup]
Numeric Solution:  
:[font = input; preserveAspect]
eq1= x''[t]+ 2 gam x'[t] + w0^2 x[t] == Q0 Cos[ w t]
:[font = text; inactive; preserveAspect]
We substitute values into the eqs to be solved by NDSolve. 
:[font = input; preserveAspect]
values={x0->1, v0->1, gam->0, Q0->1, w->1};

eqs= {eq1, x[0]==x0,x'[0]==v0} /.values 
:[font = text; inactive; preserveAspect]
We choose w0=20, and generate a numeric solution. 
:[font = input; preserveAspect]
ndsol= NDSolve[eqs /.{w0->Sqrt[2]}  ,x[t],{t,0,20}][[1]];
:[font = text; inactive; preserveAspect]
We assign values to {q2,p2}.  Note how we can easily take the derivative of a numeric function. 
:[font = input; preserveAspect]
q2[t_] = x[t] /.ndsol;
p2[t_] = q2'[t];
{q2[0.1],p2[0.1]}
:[font = text; inactive; preserveAspect]
Here is a 2-D phase-plot. This matches the above.
:[font = input; preserveAspect]
ParametricPlot[{q2[t],p2[t]},{t,0.01,20},PlotRange->All];
:[font = text; inactive; preserveAspect]
Here is a 3-D phase-plot. This matches the above. This provides us additional information about the motion. 
:[font = input; preserveAspect; endGroup; endGroup]
ParametricPlot3D[{q2[t],p2[t],t},{t,0.01,20}
  ,BoxRatios->{1,1,1}
  ,PlotPoints->100
  ,PlotRange->All
  ];
:[font = section; inactive; preserveAspect; cellOutline; startGroup]
Problem 3b:  Driving Force w/ damping
:[font = text; inactive; preserveAspect; fontColorRed = 65535; endGroup]
Copy and modify from above. 
:[font = section; inactive; Cclosed; preserveAspect; cellOutline; startGroup]
Example 3c: Non-linear term
:[font = input; preserveAspect]
Clear["Global`*"];
:[font = subsubsection; inactive; preserveAspect]
Required packages
:[font = input; preserveAspect]
Needs["Algebra`Trigonometry`"];
:[font = subsection; inactive; preserveAspect; cellOutline; startGroup]
Numeric Solution:  
:[font = input; preserveAspect]
eq1= x''[t]+ 2 gam x'[t] + w0^2 x[t] + b x[t]^3 == Q0 Cos[ w t]
:[font = text; inactive; preserveAspect]
We substitute values into the eqs to be solved by NDSolve. 
:[font = input; preserveAspect]
values={x0->0, v0->0.001, gam->0, Q0->0, w->1, b->+0.05};

eqs= {eq1, x[0]==x0,x'[0]==v0} /.values 
:[font = text; inactive; preserveAspect]
We choose w0=20, and generate a numeric solution. 
:[font = input; preserveAspect]
ndsol= NDSolve[eqs /.{w0->I 2}  ,x[t],{t,0,100}][[1]];
:[font = text; inactive; preserveAspect]
We assign values to {q2,p2}.  Note how we can easily take the derivative of a numeric function. 
:[font = input; preserveAspect]
q2[t_] = x[t] /.ndsol;
p2[t_] = q2'[t];
{q2[0.1],p2[0.1]}
:[font = text; inactive; preserveAspect]
Here is a 2-D phase-plot. This matches the above.
:[font = input; preserveAspect]
ParametricPlot[{q2[t],p2[t]},{t,0.01,100},PlotRange->All];
:[font = text; inactive; preserveAspect]
Here is a 3-D phase-plot. This matches the above. This provides us additional information about the motion. 
:[font = input; preserveAspect; endGroup; endGroup]
ParametricPlot3D[{q2[t],p2[t],t},{t,0.01,100}
  ,BoxRatios->{1,1,1}
  ,PlotPoints->500
  ,PlotRange->All
  ];
:[font = section; inactive; preserveAspect; cellOutline; startGroup]
Problem 3d: Non-linear term, driving force, w/ damping
:[font = text; inactive; preserveAspect; fontColorRed = 65535; endGroup]
Copy and modify from above. 
^*)