Physics 229& 100 Homework #7
Name: Nasser Abbasi

1. Consider the problem of a mass m on a spring of spring constant k subject to a force which is zero until t=0 and is equal to Sin[ωt] for t>0. We want to consider cases where ω=ω\(TraditionalForm\`\_0\)=(Sqrt[k/m]), and cases where ω is slightly different from (Subscript[ω, 0]). This is a simple mechanical model for an atom illuminated by near-resonant laser light.

a) To exercise your intuition, first sketch what you think the two solutions will look like; use PowerPoint or some other drawing package ( Paint if you have to) , and include the sketch of your expectations here.

For the case when the forcing frequency marches the natural frequency, resonance will occure. The mass will obtain maximum amplitude. When the forcing frequnecy is different from the natural frequency, then the mass will vibrate in sinuosidel manner (harmonic motion). This is sketch using paint (sorry about the bad quality, hard to write/draw using a mouse).
the y-axis is the displacement y(t). We see that when w=x_n, then the displacement will continue to grow.

b) Then use the greens function derived in class to make a plot of the position as a function of time for the case  ω=ω\(TraditionalForm\`\_0\)=2, m=1. Contrast this with the case in which the force is Sin[0.9*(Subscript[ω, 0])t] for t>0, i.e. off resonance illumination.

Remove["Global`*"];
Needs["Graphics`Colors`"];

greenFunc = 
      Sin[(((Sqrt[k]) (t  -  τ))/(Sqrt[m]))]/((Sqrt[k]) (Sqrt[m]));    
forcingFunction  =  
      Sin[ω τ];
y[
        t_] := (∫_0^t) 
          forcingFunction  greenFunc \[DifferentialD]τ;  (*this is the response to the forcing function*)

Due to the divide by zero problem when w=2 exactly (resonance), will plot value very near resonance first. 
Since (Subscript[w, n]) = ((Sqrt[k])/(Sqrt[m])) then k=m ((Subscript[w, n])^2)m , hence for w=2, m=1, we get k=4

p1 = 
        Plot[Evaluate[
            y[t] /. {m -> 1, ω -> 1.99999, 
                k -> 4}], {t, 0, 60}, 
          AxesLabel -> {"t", "y(t)"}, PlotStyle -> Red, 
          DisplayFunction -> Identity];
Show[p1, 
      DisplayFunction -> $DisplayFunction]

Notice also there is what is called 'beats' within the exponentional increasing motion. We see this by plotting for large t

Plot[
      Evaluate[y[t] /. {m -> 1, ω -> 1.99999, 
            k -> 4}], {t, 0, 400}, 
      AxesLabel -> {"t", "y(t)"}]

Now plot for 0.9 w

p2 = 
        Plot[Evaluate[
            y[t] /. {m -> 1, ω -> 0.9*2, k -> 4}], {t, 
            0, 65}, AxesLabel -> {"t", "y(t)"}, 
          DisplayFunction -> Identity];
Show[p2, DisplayFunction -> $DisplayFunction]

Put both plots together

Show[{p1, p2}, 
      PlotLabel -> "Comparing y(t) for resonance (RED) to non-resonance", 
      DisplayFunction -> $DisplayFunction]

2. Consider a linear  molecule as shown above  ( like acetylene). Write the Lagrangian for motions along the line in terms of the masses and spring constants. There are two different kinds of bonds, so there are two different spring constants. To make this very clear, copy the graphic into PowerPoint and add text to denote spring constants kHC and kCC , and masses mC and mH, and then re-import the graphic by copying and pasting. 



a) Use EulerEquations to generate the equations of motion.

Remove["Global`*"]

    
v1 = x[1]'[t];
v2 = x[2]'[t];
v3 = x[3]'[t];
v4 = x[4]'[t];

    
 KE = (1/2) mH (v1^2)  +  (1/2) mC (v2^2)  + (1/2) 
          mC (v3^2)  + (1/2) mH (v4^2)

PE =  (1/2) kHC (((x[1][t] - x[2][t]))^2) + (1/2) 
          kCC (((x[2][t] - x[3][t]))^2) + (1/2) 
          kHC (((x[3][t] - x[4][t]))^2)

L = KE - PE // Simplify

Needs["Calculus`VariationalMethods`"]
deps = {x[1][t], x[2][t], x[3][t], x[4][
            t]};
(eqs = EulerEquations[L, deps, t]) // MatrixForm

I can also display the above EQM is a more standard format as follows (even though this is not needed to procceed)

Needs["LinearAlgebra`MatrixManipulation`"]
vars = {x[1]''[t], x[2]''[t], x[3]''[t], x[4]''[
            t]};

    
{massMatrix, b} = 
        LinearEquationsToMatrices[eqs, vars];

    
{stiffnessMatrix, doNotCare} = 
        LinearEquationsToMatrices[b, deps];

    
Print[MatrixForm[massMatrix], 
      MatrixForm[Transpose[List[vars]]], "+", 
      MatrixForm[stiffnessMatrix], 
      MatrixForm[Transpose[List[deps]]], "== 0"]

b) Try to use DSolve to solve the resulting system of equations.

(sol = Flatten[DSolve[eqs, deps, t]]) // MatrixForm

c) Depending on which version of Mathematica you are running, you will either paralyze the machine or generate quite a lot of gibberish. Now convert the systems of equations to an eigenvalue problem and determine the eigenfrequencies and normal modes of the molecule. This problem is slightly different from the previous ones we discussed because the masses are not all identical, so be careful about the relationship between the eigenvalues and the eigenfrequencies.

Substitute trial solution x(t)= a exp(i w t) into the system of equations

(trialsol = 
          Thread[Table[x[i], {i, 1, 4}] -> 
              Table[Function[t, Evaluate[a[i] E^(I ω t)]], {i, 
                  1, 4}]]) // TableForm

(eqOfMotionAfterSub = eqs /. trialsol) // TableForm

Eliminate common factor (e^iwt)

(eqOfMotionAfterSub = 
          eqOfMotionAfterSub /. 
            qx_ == 0 :> 
              Simplify[qx*E^(-I ω t)] == 0) // 
      TableForm

eigVectorComponents = {a[1], a[2], a[3], 
          a[4]};
{massStiffnessMatrix, doNotCare} = 
        LinearEquationsToMatrices[eqOfMotionAfterSub, 
          eigVectorComponents];
Print[MatrixForm[massStiffnessMatrix], 
      MatrixForm[
        Transpose[List[eigVectorComponents]]], "== 0"]

Now Solve Det[K-w^2 M]==0 to find eigen frequencies. Then for each eigenFrequency, solve for eigenMode V from equation (K-w^2M)*V=0
The matrix I called massStiffnessMatrix above is the same as   [ k-w^2 M ] 
But first replace w^2 by λ

charEquation = 
        massStiffnessMatrix /. (ω^2) -> λ;
charEquation = Det[charEquation] == 0

eigFreq = 
        Flatten[Solve[charEquation, λ]];
eigFreq = eigFreq[[All, 2]]

eigFreq = (Sqrt[eigFreq])

For each eigenfrequency, sub back into m2 and find eigen vectors. Use NullSpace for this

eigVector = Table[0, {4}];

    
Do[Print["Normal Mode for freq ", eigFreq[[i]], " is ", 
        eigVector[[i]] = 
          Simplify[
            NullSpace[
              massStiffnessMatrix /. ω -> 
                  eigFreq[[i]]]]], {i, 1, 4}]

To see the eigenmodes better, I use some numerical values for mass and stiffness and then do an animation

(*pick some values *) 

kHC = 5;

    kCC = 2;

    mC = 50;

    mH = 20;


    N[eigVector] // MatrixForm

This below is extra. It is animation of each eigen mode.