Spring Mass Differential Equation

Routine Name: Spring Mass Equation

Author: Kyle Hovey

Language: C++ Description/Purpose: A spring-mass system with mass \( m \), damping coefficient \( \gamma \), and spring constant \( k \) is described by the differential equation

\[ m y''(t) + \gamma y'(t) + k y(t) = f(t) \]

where \( f(t) \) is defined to be the driving function for the system. This equation is derived from a sum of forces, so \( f(t) \) is equal to the force left over when all other forces cancel out.

Input:

Mass m
Damping Coefficient gamma
Spring Constant k
Initial Conditions \( y_o, y’_o \)

Output:

A function that takes time as input and returns the position of the mass at that time.

Usage/Example:

To use this function, you have to include it in the file you need to call it from.

#include "springmass/springmass.h"
#include <limits>
#include <iostream>

int main() {
  // Constants and initial values
  const double
      m = 10,
      gamma = 1,
      k = 1,
      yo = 0,
      dyo = 1;

  // Generate the analytic solution
  auto y = genSpringMass<double>(yo, dyo, m, gamma, k);

  // Output position at time t for 3 seconds (damped oscillation)
  for (double t = 0; t < 3; t += 0.1) {
    std::cout << y(t) << std::endl;
  }

  return EXIT_SUCCESS;
}

Output:

Implementation/Code:

The spring-mass differential equation is a second degree linear equation, which will have two roots to its corresponding characteristic equation. These roots should be allowed to have complex domain as a spring mass system is one that often tends to be oscillatory. The complimentary solution will have the form:

\[ y_c(t) = C_1 e^{\lambda_1 t} + C_2 e^{\lambda_2 t} \]

There is, however, a caveat. If \( \lambda_1 = \lambda_2 \) then there needs to be an additional \( u(t) \) that one of the solutions must be multiplied by to guarantee linear independence of the eigenfunctions of the differential equation.

\[ y_c(t) = C_1 e^{\lambda_1 t} + C_2 u(t) e^{\lambda_1 t} : u(t) = t \]

Caveat:

This covers all possible cases of the eigenvalues for this second order differential equation. Sometimes, people will split the former case into two parts, one where you have real non-distinct roots, one where you have complex non-distinct roots. The latter involves writing the analytic solution in terms of \( sin(x) \) and \( cos(x) \), which is equivalent to the real non-distinct case when using Euler’s Identity \( e^{ix} = cos(x) + i sin(x) \). Since the function I am writing only returns a value, and does not care how the solution is written, I am only using the two cases of equivalent roots and non-equivalent roots.

The particular solution would involve quite a bit of integration, and most likely some clever use of the Wronskian determinant via the method of Variation of Parameters. As we are not yet equipped to compute such integration in code in this class, I have left the particular solution as a stub to be filled out later.

Solving for the analytic particular solution to the homogenous part of the spring-mass differential equation, we can eliminate \( C_1 \) and \( C_2 \) by algebraicly determining their equivalents with respect to \( y_o \) and \( y’_o \). The coefficients will be different for distinct vs repeated roots, so a ternary statement is used in the following code to handle either case.

All this leaves for programming is determining the eigenvalues and leading coefficients, then to construct a function to determine the position of the mass at time t:

#include <functional>
#include <complex>
#include <iostream>

// We will be using this type a lot, make an alias for it
template <typename T>
using endomorphism = std::function<T(T)>;

// Make an alias for the zero function
template <typename T>
const endomorphism<T> zeroFunction = endomorphism<T>(
  /**
   * @param t Temporal input
   * @return Zero cast into the template type
   */
  [](T t) {
    (void) t;
    return (T) 0;
  }
);

template <typename T>
/**
 * Generate an analytic solution to a spring mass system that gives
 * the position of the mass at any time t.
 * @param yo Initial position (t = 0)
 * @param dyo Initial derivative (t = 0)
 * @param m The mass
 * @param gamma Damping constant (air resistance, etc...)
 * @param k Spring constant
 * @param endomorphism f Driving function (non-homogenous part)
 * @return Solution to spring-mass equation
 */
endomorphism<T> genSpringMass(
    T yo = 0,
    T dyo = 0,
    T m = 1,
    T gamma = 1,
    T k = 1,
    endomorphism<T> f = zeroFunction<T>
) {
  // Compute the eigenvalues and allow for complex solutions
  auto desc = std::sqrt(
      (std::complex<T>) (std::pow(gamma, 2) - 4 * m * k)
  );

  std::complex<T> eigens[2] = {
    (-gamma + desc) / std::complex<double>(2 * m),
    (-gamma - desc) / std::complex<double>(2 * m)
  };

  // Determine whether or not we have real or repeated roots
  auto distinct = eigens[0] != eigens[1];

  // Determine the leading coefficients of the solution
  auto denom = eigens[1] - eigens[0];

  std::complex<T> coeffs[2] = {
    distinct ? (eigens[1] * yo - dyo) / denom : yo,
    distinct ? (-eigens[0] * yo + dyo) / denom : dyo - eigens[0] * yo
  };

  // Determine additional u(t) required for linear independance
  const endomorphism<T> u(
    distinct ?
      (endomorphism<T>) [](T t) { (void) t; return 1; } :
      (endomorphism<T>) [](T t) { return t; }
  );

  // Stub for determining particular solution
  // TODO: Find particular solution
  (void) f;
  endomorphism<T> particular = zeroFunction<T>;

  return endomorphism<T>(
    /**
     * @param t Temporal input
     * @return Position of mass in system at t
     */
    [=] (T t) {
      /*
       * Since the characteristic equation has real coefficients, both
       * eigenvalues will be complex conjugates of each other. Since
       * conjugation is preserved under exponentiation, we are adding
       * two conjugates and thus cancelling the imaginary part. This
       * means that we only have to take the real part of the solution
       * (assuming that the particular solution has real range).
       */
      return (
        coeffs[0] * std::exp(eigens[0] * t) +
        coeffs[1] * std::exp(eigens[1] * t) * u(t) +
        particular(t)
      ).real();
    }
  );
}

Numerical Solutions (Math 5620)

Spring Mass Differential Equation

\[ m y''(t) + \gamma y'(t) + k y(t) = f(t) \]

\[ y_c(t) = C_1 e^{\lambda_1 t} + C_2 e^{\lambda_2 t} \]

\[ y_c(t) = C_1 e^{\lambda_1 t} + C_2 u(t) e^{\lambda_1 t} : u(t) = t \]