## Mass-Spring-Damper System Modeling

### Introduction

Consider an ideal mass-spring-damper system with mass m (in kg

, spring constant k (in N/m) and damping coefficient R (in N-s/m).

This system can be described with the following formula:

. We use the ode function to define the equation of motion and the initial conditions.

f1 := ode({m*x''(t) + R*x'(t) + k*x(t), x(0) = 0 , x'(0) = 1}, x(t))

Assuming m = 10, R = 1 and k = 10, the equation can be re-written as

. We substitute these values in our differential equation expression and solve the equation.

f2 := subs(f1, m = 10 , R = 1 , k = 10)

f3 := solve(f2);

As we can see above, the solution has an exponential decay term and a sinusoidal term. This implies that the system will show an underdamped response to a perturbation, as can be seen in the Figure 1 below.

p := plot::Point2d([t, f3[1]], t=0..50):

c := plot::Curve2d([t,f3[1]], t = 0..a, a= 0..50):

plot(p,c);

Note that Figure 1 was customized by modifying line properties and adding a title within the notebook interface’s object browser. These properties could have also been specified programmatically in the plot commands above.

Figure 2 shows an animation of the oscillating mass. This animation was created in a separate notebook and copied into this notebook. The ability to copy graphics and animations between notebooks is useful when users wish to document supporting analysis without having the code displayed in the notebook.

The code used to create this animation is available here.

### Frequency response analysis

We will now analyze the frequency response of our mass-spring-damper system. We begin by deriving the Laplace transform of our mass-spring-damper system, and then generate Bode plots (magnitude vs. frequency and phase vs. frequency). We perform these tasks within a custom function called FreqResp, which outputs the bode plots for a given ODE input. The Bode plot indicates that the cutoff frequency of our mass-spring-damper system is approximately 1 rad/sec.

fn := 10*x''(t) + x'(t) + 10*x(t) - h(t):

FreqResp(fn)