An Introduction to Second Order Linear Homogeneous Differential Equations and Harmonic Motion

By Debra Woods

Introduction

In this lesson, we will discuss one application of the use of second order linear differential equations with constant coefficients. 
a y'' + b y' + c y = f(x)
The application of interest in this case is harmonic motion.  We will limit our discussion to free vibrations, that is, when f(x) = 0.

We will begin by solving for the general case of 
a y'' + b y' + c y = 0  
To do this, we first form the characteristic equation, 	
a (r^2) + b r + c = 0  
 which can be solved by using the quadratic formula.  It is possible for Mathematica to solve this for us, but the form of the equation is not necessarily user friendly, that is, it may involve complex exponents.

Examples of General Solutions

Consider the following example, y'' + 4 y' + 3 y = 0.

a = 1; b = 4; c = 3;

Clear[y, x]sol = x /. Solve[a (x^2) + b x + c == 0, x]; disc = (b^2) - 4 a c; Which[disc > 0, 
      y[x_] = k[1] Exp[x sol[[1]]] + 
          k[2] Exp[x sol[[2]]], 
      disc == 0, 
      y[x_] = k[1] Exp[x sol[[1]]] + 
          k[2] x Exp[x sol[[1]]], 
      disc < 0, 
      y[x_] = Exp[Re[sol[[1]]]] 
          (k[1] Sin[
                  x Im[sol[[1]]]] + 
              k[2] Cos[x Im[sol[[1]]]])
            ]; Print["y[x] =", y[x]];

The above characteristic equation delivered two distinct real roots 
((Subscript[r,  1]) , r \(TraditionalForm\`\_2\) ) thus producing a solution to the corresponding differential equation of the form
			 y = k \(TraditionalForm\`\_1\) e \(TraditionalForm\`\^\(\(\(r\ \)\_1\) x\)\) + k \(TraditionalForm\`\_2\) e  \(TraditionalForm\`\^\(\(\(r\ \)\_2\) x\)\)

Next, we will consider the case involving only one distinct real root.  In this case, 
y'' + 2 y' + y = 0.

a = 1; b = 2; c = 1;

Clear[y, x]sol = x /. Solve[a (x^2) + b x + c == 0, x]; disc = (b^2) - 4 a c; Which[disc > 0, 
      y[x_] = k[1] Exp[x sol[[1]]] + 
          k[2] Exp[x sol[[2]]], 
      disc == 0, 
      y[x_] = k[1] Exp[x sol[[1]]] + 
          k[2] x Exp[x sol[[1]]], 
      disc < 0, 
      y[x_] = Exp[Re[sol[[1]]] x] 
          (k[1] Sin[
                  x Im[sol[[1]]]] + 
              k[2] Cos[x Im[sol[[1]]]])
            ]; Print["y[x] =", y[x]];

Notice that because only one solution is obtained, it had to be multiplied by x in order to obtain two linearly independent solutions in the fundamental solution set.  To convince yourself of this, review the section on reduction of order in your textbook.  Equations yielding only one solution will always produce solutions of the form 
 y = k \(TraditionalForm\`\_1\) e \(TraditionalForm\`\^\(\(\(r\ \)\_1\) x\)\) + k \(TraditionalForm\`\_2\)x e  \(TraditionalForm\`\^\(\(\(r\ \)\_1\) x\)\)'.

The last case is the case where the discriminant is negative.  This of course produces complex solutions, and is the case where Mathematica gives us solutions involving complex exponents, that is, E(iθ).

Clear[y, x, a]

By the use of Eulers formula E(iθ) = cos θ + i sin θ , we can change the form of the solution to
		 y = e ((Subscript[k, 1]) cos(qx) + (Subscript[k, 2]) sin(qx)) . 
where the complex roots to the characteristic equation are p ± q i.
The equation of the next example is y'' + y' + y = 0.

a = 1; b = 1; c = 1;

Clear[y, x]sol = x /. Solve[a (x^2) + b x + c == 0, x]; disc = (b^2) - 4 a c; Which[disc > 0, 
      y[x_] = k[1] Exp[x sol[[1]]] + 
          k[2] Exp[x sol[[2]]], 
      disc == 0, 
      y[x_] = k[1] Exp[x sol[[1]]] + 
          k[2] x Exp[x sol[[1]]], 
      disc < 0, 
      y[x_] = Exp[Re[sol[[1]]] x] 
          (k[1] Sin[
                  x Im[sol[[1]]]] + 
              k[2] Cos[x Im[sol[[1]]]])
            ]; Print["y[x] =", y[x]];

Harmonic Motion

Now we will divert our attention to the application in question, harmonic motion.  The differential equation of the free motion for a vibrating spring is of the form a y'' + b y' + c y = 0.  Physically, the coefficients, a, b, and c correspond to mass,damping constant, and the spring constant respectively. First, we will consider the case where b = 0.  This is the case of undamped motion.  When solving a y'' + c y = 0, if we divide by a, we get y'' + (c/a)y = 0, which we can rewrite as 
y'' + ωy = 0.  The equation will clearly yield a solution of the form y = k cos(ωx) + k' sin(ωx).  Let's test this with Mathematica. and see what the equation of motion looks like.

Example of Harmonic Motion

The Case of Undamped Motion

Example:  A 5 kg mass is attached to a spring of length .3 m, stretching the spring to a length of .7 m.  The mass is then pulled downward and additional .3 m and released.  Find the equation of motion for the mass.

To find the differential equation of motion, we must first find a, b, and c.  The mass a is given in the problem and is 5,  b = 0 for undamped motion, and to find c, the spring constant, we need to used the given information and Hooke's Law F = c x.

Clear[a, b, c, y, x]\(c = \(5\ 9.8\)\/ .4\)

The displacement of the spring due to attaching the mass is .4 m, so c = F/x = (5)(9.8)/.4 = 122.5 kg/sec¤.

a = 5; b = 0; Clear[y, x]sol = x /. Solve[a (x^2) + b x + c == 0, x]; disc = (b^2) - 4 a c; Which[disc > 0, 
      y[x_] = k[1] Exp[x sol[[1]]] + 
          k[2] Exp[x sol[[2]]], 
      disc == 0, 
      y[x_] = k[1] Exp[x sol[[1]]] + 
          k[2] x Exp[x sol[[1]]], 
      disc < 0, 
      y[x_] = Exp[Re[sol[[1]]] x] 
          (k[1] Sin[
                  x Im[sol[[1]]]] + 
              k[2] Cos[x Im[sol[[1]]]])
            ]; Print["y[x] =", y[x]];

Now, we will solve the IVP using the initial conditions,
 y(0) = .4,(initial position) , and y'(0) = 0 (initial velocity).

y0 =  .4; y1 = 0; ksols=Solve[{y[0] == y0,y'[0]==y1},{k[1], k[2]}];ysolved[x_] = 
      y[x] /. ksols[[1]]

The graph of the equation of motion is

Plot[ysolved[x], {x, 0, 2 π}, 
      PlotStyle -> {RGBColor[0, 0, 1]}];

This graph is a good example of undamped motion.  Now, let's consider damped motion.  For this, there are three cases.  We will first consider the case of overdamped motion, when the discriminant of the characteristic equation is positive.  In the initial work above, we saw that this lead to two distinct roots.  Let's change the above example to include a damping constant, b = 70.

Clear[a, b, c, y, x]\(c = \(5\ 9.8\)\/ .4\)a = 5; b = 70;

Clear[y, x]sol = x /. Solve[a (x^2) + b x + c == 0, x]; disc = (b^2) - 4 a c; Which[disc > 0, 
      y[x_] = k[1] Exp[x sol[[1]]] + 
          k[2] Exp[x sol[[2]]], 
      disc == 0, 
      y[x_] = k[1] Exp[x sol[[1]]] + 
          k[2] x Exp[x sol[[1]]], 
      disc < 0, 
      y[x_] = Exp[Re[sol[[1]]] x] 
          (k[1] Sin[
                  x Im[sol[[1]]]] + 
              k[2] Cos[x Im[sol[[1]]]])
            ]; Print["y[x] =", y[x]];

y0 =  .4; y1 = 0; ksols=Solve[{y[0] == y0,y'[0]==y1},{k[1], k[2]}];ysolved[x_] = 
      y[x] /. ksols[[1]]

The graph of the equation of motion is

Plot[ysolved[x], {x, 0, π}, 
      PlotStyle -> {RGBColor[0, 0, 1]}];

Notice that the motion is no longer oscillatory.

The next case is the case where the discriminant equals 0.  This is called critically damped motion.  Let's set the damping constant b equal to 4 * a * c, forcing the discriminant to 0.

Clear[a, b, c, y, x]\(\(c = \(5\ 9.8\)\/ .4; \)\)a = 5; b = (Sqrt[4 a c]);

Clear[y, x]sol = x /. Solve[a (x^2) + b x + c == 0, x]; disc = (b^2) - 4 a c; Which[disc > 0, 
      y[x_] = k[1] Exp[x sol[[1]]] + 
          k[2] Exp[x sol[[2]]], 
      disc == 0, 
      y[x_] = k[1] Exp[x sol[[1]]] + 
          k[2] x Exp[x sol[[1]]], 
      disc < 0, 
      y[x_] = Exp[Re[sol[[1]]] x] 
          (k[1] Sin[
                  x Im[sol[[1]]]] + 
              k[2] Cos[x Im[sol[[1]]]])
            ]; Print["y[x] =", y[x]];

y0 =  .4; y1 = 0; ksols=Solve[{y[0] == y0,y'[0]==y1},{k[1], k[2]}];ysolved[x_] = 
      y[x] /. ksols[[1]]

The graph of the equation of motion is

Plot[ysolved[x], {x, 0, π}, 
      PlotStyle -> {RGBColor[0, 0, 1]}];

The last case is that of under damped motion, that is, when the discriminant is less than zero.  Let's give our equation a small damping constant, say b = 30.

Clear[a, b, c, y, x, sol]\(\(c = \(5\ 9.8\)\/ .4; \)\)a = 5; b = 30;

Clear[y, x]sol = x /. Solve[a (x^2) + b x + c == 0, x]; disc = (b^2) - 4 a c; Which[disc > 0, 
      y[x_] = 
        k[1] Exp[x sol[[1]]] + 
          k[2] Exp[x sol[[2]]]; , 
      disc == 0, 
      y[x_] = 
        k[1] Exp[x sol[[1]]] + 
          k[2] x Exp[x sol[[1]]]; , 
      disc < 0, {resol = Re[sol[[1]]]; 
        , y[x_] = 
          Exp[resol x] 
            (k[1] Sin[
                    x Im[sol[[1]]]] + 
                k[2] Cos[
                    x Im[sol[[1]]]]); 
        }];

Print["y[x] = ", y[x]]

y0 =  .4; y1 = 0; ksols=Solve[{y[0] == y0,y'[0]==y1},{k[1], k[2]}];ysolved[x_] = 
      y[x] /. ksols[[1]]

The graph of the equation of motion is

Plot[ysolved[x], {x, 0, 6}, PlotRange -> All, 
      PlotStyle -> {RGBColor[0, 0, 1]}];

It is difficult on this graph to see the oscillatory motion.  Let's zoom in and view the graph.

Plot[ysolved[x], {x, 0, 6}, PlotStyle -> {RGBColor[0, 0, 1]}];

Here we can see that the under damping has caused a few oscillations.  The larger the damping constant, the fewer oscillations, and the closer the damping constant gets to zero, the closer the equation of motion gets to undamped motion.