Chapter 28, Problem 49P

The dynamics of a forced spring-mass-damper system can be represented by the following second-order ODE:

$m\frac{d^{2}x}{d{t}^2}+c\frac{dx}{dt}+k_{1}x+k_3{x}^3=P\cos \left(\omega t\right)$

where $m=1$ kg, $c=0.4\text{ N}$? s/m, $P=0.5\text{ N}$, and $\omega =0.5\text{/s}$. Use a numerical method to solve for displacement $\left(x\right)$ and velocity $\left(v=dx\text{/}dt\right)$ as a function of time with the initial conditions $x=v=0$. Express your results graphically as time-series plots (x and v versus t) and a phase plane plot (v versus x). Perform simulations for both (a) linear $\left(k_1=1;k_3=0\right)$ and (b) nonlinear $\left(k_1=1;k_3=0.5\right)$ springs.

To determine

To calculate: The displacement and velocity as a function of time for linear system where $k_1=1\text{ and }{k}_3=0$ if the dynamic of a forced spring-mass-damper system is given as $m\frac{d^{2}x}{d{t}^2}+c\frac{dx}{dt}+k_{1}x+k_3{x}^3=P\cos \left(\omega t\right)$.

Answer to Problem 49P

Solution:

The first few solutions for displacement and velocity as a function of time for liner system is,

t	x	v
0	0	0
0.125	0.003836	0.060764
0.25	0.015019	0.117402
0.375	0.032976	0.169015
0.5	0.05703	0.214828
0.625	0.086414	0.254195
0.75	0.120289	0.286602
0.875	0.157759	0.311673
1	0.197891	0.329164
1.125	0.239729	0.338967
1.25	0.282313	0.341104
1.375	0.324691	0.335717
1.5	0.365939	0.323069
1.625	0.405171	0.303527
1.75	0.441553	0.277555
1.875	0.474314	0.245705
2	0.50276	0.208601
2.125	0.526274	0.16693
2.25	0.544332	0.121428
2.375	0.556503	0.072863
2.5	0.562454	0.02203
2.625	0.56195	-0.03027
2.75	0.554859	-0.08324
2.875	0.541146	-0.13608
3	0.520874	-0.18805
3.125	0.494199	-0.23842
3.25	0.461364	-0.28651
3.375	0.422694	-0.33168
3.5	0.378588	-0.37339
3.625	0.329513	-0.41112
3.75	0.275993	-0.44444
3.875	0.218601	-0.47301
4	0.157952	-0.49653
4.125	0.094688	-0.5148
4.25	0.029475	-0.52771
4.375	-0.03701	-0.53518
4.5	-0.1041	-0.53725
4.625	-0.1711	-0.53399
4.75	-0.23738	-0.52558
4.875	-0.30229	-0.51221
5	-0.36524	-0.49416

Explanation of Solution

Given Information:

The dynamic of a forced spring-mass-damper system is given as,

$m\frac{d^{2}x}{d{t}^2}+c\frac{dx}{dt}+k_{1}x+k_3{x}^3=P\cos \left(\omega t\right)$

The values,

$\begin{array}{l}m=1\text{ kg}\\ c=0.4\text{ N}\cdot \text{s/m}\\ P=0.5\text{ N}\\ \omega =0.5\text{ /s}\end{array}$

The initial condition, $x=v=0$.

Formula used:

The fourth-order RK method for $\frac{dy}{dt}=f\left(t,y\right)$ is,

$\begin{array}{c}{y}_{n+1}=y_n+\frac{1}{6}\left(k_1+2k_2+2k_3+k_4\right)\\ t_{n+1}=t_n+h\end{array}$

Where,

$\begin{array}{c}{k}_1=hf\left(t_n,y_n\right)\\ k_2=hf\left(t_n+\frac{h}{2},y_n+\frac{k_1}{2}\right)\\ k_3=hf\left(t_n+\frac{h}{2},y_n+\frac{k_2}{2}\right)\\ k_4=hf\left(t_n+h,y_n+k_3\right)\end{array}$

Calculation:

Consider the dynamic of a forced spring-mass-damper system,

$m\frac{d^{2}x}{d{t}^2}+c\frac{dx}{dt}+k_{1}x+k_3{x}^3=P\cos \left(\omega t\right)$

As $\frac{dx}{dt}=v$, replace $\frac{dx}{dt}$ by $v$ in the above equation,

$m\frac{dv}{dt}+cv+k_{1}x+k_3{x}^3=P\cos \left(\omega t\right)$

Divide both the sides of above equation by m,

$\frac{dv}{dt}=\frac{P\cos \left(\omega t\right)}{m}-\frac{cv}{m}-\frac{k_{1}x}{m}-\frac{k_3{x}^3}{m}$

Now, substitute the values $m=1\text{ kg, }c=0.4\text{ N}\cdot \text{s/m, }P=0.5\text{ N and }\omega =0.5\text{ /s}$ in the above equation,

$\begin{array}{c}\frac{dv}{dt}=\frac{0.5\cos \left(0.5t\right)}{1}-\frac{0.4v}{1}-\frac{k_{1}x}{1}-\frac{k_3{x}^3}{1}\\ \frac{dv}{dt}=0.5\cos \left(0.5t\right)-0.4v-k_{1}x-k_3{x}^3\end{array}$

For linear, substitute $k_1=1$ and $k_3=0$ in the equation,

$\frac{dv}{dt}=0.5\cos \left(0.5t\right)-0.4v-x$

Use VBA code for RK4 method as below to solve for x and v,

Code:

OptionExplicit

' Generate subfunction RK4SysTest()

Sub RK4SysTest()

'Declare the variables as integer

DimiAsInteger,m AsInteger,n AsInteger,j AsInteger

'Declare the variables as double

Dimx_iAsDouble,y_i(10)AsDouble,x_fAsDouble

DimdiffxAsDouble,x_outAsDouble

Dimx_p(200)AsDouble,y_p(200, 10)AsDouble

'Set the initial values

n =2

x_i=0

x_f=25

y_i(1)=0

y_i(2)=0

diffx=0.125

x_out=0.125

'move the values at a specific cell

Range("a3").Select

ActiveCell. Value="RK4 method"

'name each columns

ActiveCell. Offset(1, 0).Select

ActiveCell. Value="t"

ActiveCell. Offset(0, 1).Select

ActiveCell. Value="x"

ActiveCell. Offset(0, 1).Select

ActiveCell. Value="v"

' Call the function ODESolver

CallODESolver(x_i,y_i,x_f,diffx,x_out,x_p,y_p,m,n)

'Display the results in Sheet1

Sheets("Sheet1").Select

Range("a5:n205").ClearContents

Range("a5").Select

'Use for loop to store the different values of x_p

Fori=0To m

ActiveCell. Value=x_p(i)

'Use for loop to store the different values of y_p

For j =1To n

'Define the offset

ActiveCell. Offset(0, 1).Select

ActiveCell. Value=y_p(i,j)

Next j

'Define the offset

ActiveCell. Offset(1,-n).Select

Nexti

Range("a5").Select

EndSub

' Generate the subfunction ODESolver

SubODESolver(x_i,y_i,x_f,diffx,x_out,x_p,y_p,m,n)

'Declare the variables as integer

DimiAsInteger

'Declare the variables as double

Dim x AsDouble,y(10)AsDouble,x_endAsDouble

Dim h AsDouble

'Set the variables

m =0

'Store the value of x_i to x

x =x_i

'Use for loop to store the different values of y_i in y

Fori=1To n

y(i)=y_i(i)

Nexti

'Store the value of x in x_p

x_p(m)= x

'Use for loop to store the different values of y in y_p

Fori=1To n

y_p(m,i)= y(i)

Nexti

Do

'Display the result of x_end

x_end= x +x_out

If(x_end>x_f)Thenx_end=x_f

h =diffx

'Call the Integrator function

CallIntegrator(x,y,h,n,x_end)

m = m +1

'Store the value of x in x_p

x_p(m)= x

Fori=1To n

'Store the value of y in y_p

y_p(m,i)= y(i)

Nexti

'Use condition to exit the loop

If(x >=x_f)ThenExitDo

Loop

EndSub

' Generate a subfunction Integrator

SubIntegrator(x,y,h,n,x_end)

'Declare the variables

Dim j AsInteger

Dimynew(10)AsDouble

Do

'use loop for calculations

If(x_end- x < h)Then h =x_end- x

'Call the RK4Sys function

Call RK4Sys(x,y,h,n,ynew)

For j =1To n

'Store the different values of ynew to y

y(j)=ynew(j)

Next j

'Use condition to exit the loop

If(x >=x_end)ThenExitDo

Loop

EndSub

' Generate the subfunction RK4Sys

Sub RK4Sys(x,y,h,n,ynew)

'Declare the variables as integer

Dim j AsInteger

'Declare the variables as double

Dimy_m(10)AsDouble,y_e(10)AsDouble

Dim k_1(10)AsDouble,k_2(10)AsDouble,k_3(10)AsDouble,k_4(10)AsDouble

Dimslope(10)

'Call the Derivs function

CallDerivs(x,y,k_1)

For j =1To n

'Store the value of y(j) + k_1(j) * h / 2 in y_m

y_m(j)= y(j)+ k_1(j)* h /2

Next j

'Call the Derivs function

CallDerivs(x + h /2,y_m,k_2)

For j =1To n

'Store the value of y(j) + k_2(j) * h / 2 in y_m

y_m(j)= y(j)+ k_2(j)* h /2

Next j

'Call the Derivs function

CallDerivs(x + h /2,y_m,k_3)

For j =1To n

'Store the value of y(j) + k_3(j) * h in y_e

y_e(j)= y(j)+ k_3(j)* h

Next j

'Call the Derivs function

CallDerivs(x + h,y_e,k_4)

For j =1To n

'Store the value of k_1(j) + 2 * (k_2(j) + k_3(j)) + k_4(j)) / 6 in slope

slope(j)=(k_1(j)+2*(k_2(j)+ k_3(j))+ k_4(j))/6

Next j

For j =1To n

'Use formula to find ynew

ynew(j)= y(j)+ slope(j)* h

Next j

'Wrtite the next value of x

x = x + h

EndSub

' Generate the subfunction Derivs

SubDerivs(x,y,dydx)

'Write the ordinary differential equations

dydx(1)= y(2)

dydx(2)=-0.4* y(2)- y(1)+0.5* Cos(0.5* x)

EndSub

Output:

To draw the graph of the above results, follow the steps in excel sheet as given below,

Step 1: Select the cell from A4 to A205 and cell B4 to B205. Then, go to the Insert and select the scatter with smooth lines from the chart.

Step 2: Select the cell from A4 to A205 and cell C4 to C205. Then, go to the Insert and select the scatter with smooth lines from the chart.

Step 3: Select one of the graphs and paste it on another graph to merge the graphs.

The graph obtained is,

And, to draw the phase plane plot follow the steps as below,

Step 4: Select the column B and column C. Then, go to the Insert and select the scatter with smooth lines from the chart.

The graph obtained is,

To determine

To calculate: The displacement and velocity as a function of time for non-linear system where $k_1=1\text{ and }{k}_3=0.5$ if the dynamic of a forced spring-mass-damper system is given as $m\frac{d^{2}x}{d{t}^2}+c\frac{dx}{dt}+k_{1}x+k_3{x}^3=P\cos \left(\omega t\right)$.

Answer to Problem 49P

Solution:

t	x	v
0	0	0
0.125	0.003836	0.060764
0.25	0.015019	0.117402
0.375	0.032976	0.169014
0.5	0.05703	0.214821
0.625	0.086412	0.254166
0.75	0.120279	0.286506
0.875	0.157729	0.311417
1	0.19781	0.328581
1.125	0.239542	0.337784
1.25	0.28192	0.338923
1.375	0.323936	0.332004
1.5	0.36459	0.31716
1.625	0.402907	0.294658
1.75	0.437953	0.264919
1.875	0.468859	0.228524
2	0.494837	0.186221
2.125	0.515205	0.138911
2.25	0.5294	0.087637
2.375	0.536996	0.033543
2.5	0.537718	-0.02217
2.625	0.531437	-0.07829
2.75	0.518177	-0.13366
2.875	0.498099	-0.18721
3	0.47149	-0.23801
3.125	0.438743	-0.28531
3.25	0.400333	-0.32852
3.375	0.3568	-0.36723
3.5	0.308725	-0.40117
3.625	0.256712	-0.43023
3.75	0.201374	-0.45436
3.875	0.143325	-0.47361
4	0.083172	-0.48804
4.125	0.021513	-0.49771
4.25	-0.04106	-0.50268
4.375	-0.10396	-0.50296
4.5	-0.1666	-0.49854
4.625	-0.2284	-0.48938
4.75	-0.28875	-0.47544
4.875	-0.34705	-0.45667
5	-0.40271	-0.43306

Explanation of Solution

Given Information:

The dynamic of a forced spring-mass-damper system is given as,

$m\frac{d^{2}x}{d{t}^2}+c\frac{dx}{dt}+k_{1}x+k_3{x}^3=P\cos \left(\omega t\right)$

The values,

$\begin{array}{l}m=1\text{ kg}\\ c=0.4\text{ N}\cdot \text{s/m}\\ P=0.5\text{ N}\\ \omega =0.5\text{ /s}\end{array}$

And,

$k_1=1\text{ and }{k}_3=0.5$

The initial condition, $x=v=0$.

Formula used:

The fourth-order RK method for $\frac{dy}{dt}=f\left(t,y\right)$ is,

$\begin{array}{c}{y}_{n+1}=y_n+\frac{1}{6}\left(k_1+2k_2+2k_3+k_4\right)\\ t_{n+1}=t_n+h\end{array}$

Where,

$\begin{array}{c}{k}_1=hf\left(t_n,y_n\right)\\ k_2=hf\left(t_n+\frac{h}{2},y_n+\frac{k_1}{2}\right)\\ k_3=hf\left(t_n+\frac{h}{2},y_n+\frac{k_2}{2}\right)\\ k_4=hf\left(t_n+h,y_n+k_3\right)\end{array}$

Calculation:

Consider the dynamic of a forced spring-mass-damper system,

$m\frac{d^{2}x}{d{t}^2}+c\frac{dx}{dt}+k_{1}x+k_3{x}^3=P\cos \left(\omega t\right)$

As $\frac{dx}{dt}=v$, replace $\frac{dx}{dt}$ by $v$ in the above equation,

$m\frac{dv}{dt}+cv+k_{1}x+k_3{x}^3=P\cos \left(\omega t\right)$

Divide both the sides of above equation by m,

$\frac{dv}{dt}=\frac{P\cos \left(\omega t\right)}{m}-\frac{cv}{m}-\frac{k_{1}x}{m}-\frac{k_3{x}^3}{m}$

Now, substitute the values $m=1\text{ kg, }c=0.4\text{ N}\cdot \text{s/m, }P=0.5\text{ N and }\omega =0.5\text{ /s}$ in the above equation,

$\begin{array}{c}\frac{dv}{dt}=\frac{0.5\cos \left(0.5t\right)}{1}-\frac{0.4v}{1}-\frac{k_{1}x}{1}-\frac{k_3{x}^3}{1}\\ \frac{dv}{dt}=0.5\cos \left(0.5t\right)-0.4v-k_{1}x-k_3{x}^3\end{array}$

For linear, substitute $k_1=1$ and $k_3=0.5$ in the equation,

$\frac{dv}{dt}=0.5\cos \left(0.5t\right)-0.4v-x-0.5x^3$

Use VBA code for RK4 method as below to solve for x and v,

Code:

OptionExplicit

' Generate subfunction RK4SysTest()

Sub RK4SysTest()

'Declare the variables as integer

DimiAsInteger,m AsInteger,n AsInteger,j AsInteger

'Declare the variables as double

Dimx_iAsDouble,y_i(10)AsDouble,x_fAsDouble

DimdiffxAsDouble,x_outAsDouble

Dimx_p(200)AsDouble,y_p(200, 10)AsDouble

'Set the initial values

n =2

x_i=0

x_f=25

y_i(1)=0

y_i(2)=0

diffx=0.125

x_out=0.125

'move the values at a specific cell

Range("a3").Select

ActiveCell. Value="RK4 method"

'name each columns

ActiveCell. Offset(1, 0).Select

ActiveCell. Value="t"

ActiveCell. Offset(0, 1).Select

ActiveCell. Value="x"

ActiveCell. Offset(0, 1).Select

ActiveCell. Value="v"

' Call the function ODESolver

CallODESolver(x_i,y_i,x_f,diffx,x_out,x_p,y_p,m,n)

'Display the results in Sheet1

Sheets("Sheet1").Select

Range("a5:n205").ClearContents

Range("a5").Select

'Use for loop to store the different values of x_p

Fori=0To m

ActiveCell. Value=x_p(i)

'Use for loop to store the different values of y_p

For j =1To n

'Define the offset

ActiveCell. Offset(0, 1).Select

ActiveCell. Value=y_p(i,j)

Next j

'Define the offset

ActiveCell. Offset(1,-n).Select

Nexti

Range("a5").Select

EndSub

' Generate the subfunction ODESolver

SubODESolver(x_i,y_i,x_f,diffx,x_out,x_p,y_p,m,n)

'Declare the variables as integer

DimiAsInteger

'Declare the variables as double

Dim x AsDouble,y(10)AsDouble,x_endAsDouble

Dim h AsDouble

'Set the variables

m =0

'Store the value of x_i to x

x =x_i

'Use for loop to store the different values of y_i in y

Fori=1To n

y(i)=y_i(i)

Nexti

'Store the value of x in x_p

x_p(m)= x

'Use for loop to store the different values of y in y_p

Fori=1To n

y_p(m,i)= y(i)

Nexti

Do

'Display the result of x_end

x_end= x +x_out

If(x_end>x_f)Thenx_end=x_f

h =diffx

'Call the Integrator function

CallIntegrator(x,y,h,n,x_end)

m = m +1

'Store the value of x in x_p

x_p(m)= x

Fori=1To n

'Store the value of y in y_p

y_p(m,i)= y(i)

Nexti

'Use condition to exit the loop

If(x >=x_f)ThenExitDo

Loop

EndSub

' Generate a subfunction Integrator

SubIntegrator(x,y,h,n,x_end)

'Declare the variables

Dim j AsInteger

Dimynew(10)AsDouble

Do

'use loop for calculations

If(x_end- x < h)Then h =x_end- x

'Call the RK4Sys function

Call RK4Sys(x,y,h,n,ynew)

For j =1To n

'Store the different values of ynew to y

y(j)=ynew(j)

Next j

'Use condition to exit the loop

If(x >=x_end)ThenExitDo

Loop

EndSub

' Generate the subfunction RK4Sys

Sub RK4Sys(x,y,h,n,ynew)

'Declare the variables as integer

Dim j AsInteger

'Declare the variables as double

Dimy_m(10)AsDouble,y_e(10)AsDouble

Dim k_1(10)AsDouble,k_2(10)AsDouble,k_3(10)AsDouble,k_4(10)AsDouble

Dimslope(10)

'Call the Derivs function

CallDerivs(x,y,k_1)

For j =1To n

'Store the value of y(j) + k_1(j) * h / 2 in y_m

y_m(j)= y(j)+ k_1(j)* h /2

Next j

'Call the Derivs function

CallDerivs(x + h /2,y_m,k_2)

For j =1To n

'Store the value of y(j) + k_2(j) * h / 2 in y_m

y_m(j)= y(j)+ k_2(j)* h /2

Next j

'Call the Derivs function

CallDerivs(x + h /2,y_m,k_3)

For j =1To n

'Store the value of y(j) + k_3(j) * h in y_e

y_e(j)= y(j)+ k_3(j)* h

Next j

'Call the Derivs function

CallDerivs(x + h,y_e,k_4)

For j =1To n

'Store the value of k_1(j) + 2 * (k_2(j) + k_3(j)) + k_4(j)) / 6 in slope

slope(j)=(k_1(j)+2*(k_2(j)+ k_3(j))+ k_4(j))/6

Next j

For j =1To n

'Use formula to find ynew

ynew(j)= y(j)+ slope(j)* h

Next j

'Wrtite the next value of x

x = x + h

EndSub

' Generate the subfunction Derivs

SubDerivs(x,y,dydx)

'Write the ordinary differential equations

dydx(1)= y(2)

dydx(2)=-0.4* y(2)- y(1)-0.5*(y(1))^3+0.5* Cos(0.5* x)

EndSub

Output:

Few data are shown below,

To draw the graph of the above results, follow the steps in excel sheet as given below,

Step 1: Select the cell from A4 to A205 and cell B4 to B205. Then, go to the Insert and select the scatter with smooth lines from the chart.

Step 2: Select the cell from A4 to A205 and cell C4 to C205. Then, go to the Insert and select the scatter with smooth lines from the chart.

Step 3: Select one of the graphs and paste it on another graph to merge the graphs.

The graph obtained is,

And, to draw the phase plane plot follow the steps as below,

Step 4: Select the column B and column C. Then, go to the Insert and select the scatter with smooth lines from the chart.

The graph obtained is,