Example Problem: Damped Spring Motion

To demonstrate how the rungeKutta4() method can be used to solve an initial value ODE problem, we will use as an example the motion of a damped spring. Consider the spring shown in Figure 20.1. The upper end of the spring is fastened to a solid object such that it can't move. The lower end of the spring is attached to a body of mass m .

Figure 20.1. Spring configuration

When the spring is stretched a distance x from its equilibrium position, the force exerted on the mass by the spring is given by Hooke's law.

Equation 20.9

The k parameter is the spring constant. The force on the body attached to the string can also be characterized by Newton's second law

Equation 20.10

Putting Eq. (20.9) and Eq. (20.10) together we obtain the general equation for the motion of an undamped spring.

Equation 20.11

Eq. (20.11) assumes that there are no other forces acting on the spring. In reality damping forces such as friction and air resistance will slow the spring's motion. Damping forces are a function of the velocity of the spring and a damping constant, µ. When damping forces are added to Eq. (20.11), we obtain the general equation of motion for a spring.

Equation 20.12

You can see that the spring equation is a second-order ODE with time as its independent variable. What makes the spring motion example a good test case for our ODE solver development is that there is an exact solution to Eq. (20.12). If µ ² > 4 mk , the spring system is overdamped. When the mass is moved from its equilibrium position and released it will simply return asymptotically to its equilibrium position. If µ ² < 4 mk , the system is underdamped and the result will be damped harmonic motion. The mass oscillates about its equilibrium position with asymptotically declining minimum and maximum values. The general solution for the position of a mass connected to an underdamped spring is shown in Eq. (20.13).

Equation 20.13

The A and B constants in Eq. (20.13) are defined by the expressions in Eq. (20.14).

Equation 20.14

The constants C ₁ and C ₂ are derived from the initial conditions. Assume that initially the spring mass is extended a distance x from its equilibrium position and the spring velocity is zero. Under these conditions the constants take the values in Eq. (20.15).

Equation 20.15

Before we can use the rungeKutta4() method to numerically solve for underdamped spring motion, Eq. (20.12) is recast in terms of two first-order ODEs.

Equation 20.16

Equation (20.16) is integrated with respect to time to solve for x and dx / dt .

SpringODE class

To solve the damped spring initial value problem, we will write a SpringODE class that is a subclass of the ODE class. The SpringODE class will represent the equations of motion for a damped spring. In addition to the members it inherits from the ODE class, the SpringODE class defines three new fields representing the spring constant, damping constant, and mass. While reading through this section focus on the process. Even if you don't need to solve for the motion of a damped spring, you can apply the solution process described in this section to your own ODE problems.

The SpringODE class declares one constructor. The equations of motion for a damped spring consist of two first-order differential equations and zero free variables . The first thing the SpringODE constructor does is to call the ODE class constructor passing it values of 2 and 0 respectively. The SpringODE constructor then initializes the k , mu , and mass fields.

The SpringODE class overrides the getFunction() and setInitialConditions() methods of the ODE class. Since SpringODE represents an initial value problem there is no need to override the getError() method. The getFunction() method is overridden to return the right-hand side of Eq. (20.16). This evaluation will occur at various points along the range of integration. The ytmp[] array holds the value of the dependent variables at the point currently being evaluated. The setInitialConditions() method is overridden to provide the initial conditions for the spring motion. At time t = 0 the spring is at rest so dx / dt = 0. The V[] array holds the initial displacement of the spring from its equilibrium position.

The SpringODE class source code is shown next .

 package TechJava.MathLib; public class SpringODE extends ODE {   double k, mu, mass;   //  The SpringODE constructor calls the ODE   //  class constructor passing it data for   //  a damped spring. There are two first   //  order ODEs and no free variables.   public SpringODE(double k, double mu, double mass) {     super(2,0);     this.k = k;     this.mu = mu;     this.mass = mass;   }   // The getFunction() method returns the right-hand   // sides of the two first-order damped spring ODEs   // y[0] = delta(dxdt) = delta(t)*(-k*x - mu*dxdt)/mass   // y[1] = delta(x) = delta(t)*(dxdt)   public void getFunction(double x, double dy[],                           double ytmp[]) {     dy[0] = -k*ytmp[1]/mass - mu*ytmp[0]/mass;     dy[1] = ytmp[0];   }   // This method initializes the dependent variables   // at the start of the integration range.   public void setInitialConditions(double V[]) {     setOneY(0, 0, 0.0);     setOneY(0, 1, V[0]);     setOneX(0, 0.0);   } } 

Solving the Spring Motion ODE

Now let us apply the fourth-order Runge -Kutta solver to compute the motion of a damped spring. We will write a driver program named RK4Spring.java that will create a SpringODE object and call the rungeKutta4() method on that object. The mass , mu , and k parameters are given values representing an underdamped spring. The initial conditions are set such that the spring is extended 0.2 meters from its equilibrium position. The ODE is integrated from t = 0 to t = 5.0 seconds using a step size of 0.1 seconds. The RK4Spring class source code is shown next.

 import TechJava.MathLib.*; public class RK4Spring {   public static void main(String args[]) {     //  Create a SpringODE object     double mass = 1.0;     double mu = 1.5;     double k = 20.0;     SpringODE ode = new SpringODE(k, mu, mass);     //  load initial conditions.  The spring is     //  initially stretched 0.2 meters from its     //  equilibrium  position.     double V[] = {-0.2};     ode.setInitialConditions(V);     //  Solve the ODE over the desired range using     //  a constant step size.     double dx = 0.1;     double range = 5.0;    int numSteps = ODESolver.rungeKutta4(ode, range, dx);     //  Print out the results     System.out.println("i   t   dxdt   x");     for(int i=0; i<numSteps; ++i) {       System.out.println(         ""+i+" " + ode.getOneX(i) +         " " + ode.getOneY(i,0) +" " + ode.getOneY(i,1));     }   } } 

Output ”

Rather than list a long table of output values, a plot was created that displays the spring position as a function of time as computed by the RK4Spring.java program. Also shown on the plot is the exact solution of the ODE. You can see from Figure 20.2 that the rungeKutta4() method did an excellent job of reproducing the exact solution of the spring equation. The spring oscillates around the position x = 0 with asymptotically diminishing maximum and minimum values.

Figure 20.2. Spring position as a function of time