Preface

This tutorial was made solely for the purpose of education and it was designed for students taking Applied Math 0330. It is primarily for students who have very little experience or have never used Mathematica before and would like to learn more of the basics for this computer algebra system. As a friendly reminder, don't forget to clear variables in use and/or the kernel.

Finally, the commands in this tutorial are all written in bold black font, while Mathematica output is in regular fonts. This means that you can copy and paste all commands into Mathematica, change the parameters and run them. You, as the user, are free to use the scripts to your needs for learning how to use the Mathematica program, and have the right to distribute this tutorial and refer to this tutorial as long as this tutorial is accredited appropriately.

Finally, the commands in this tutorial are all written in bold black font, while Mathematica output is in normal font. This means that you can copy and paste all commands into Mathematica, change the parameters and run them. You, as the user, are free to use the scripts for your needs to learn the Mathematica program, and have the right to distribute this tutorial and refer to this tutorial as long as this tutorial is accredited appropriately.

Return to computing page for the first course APMA0330
Return to computing page for the second course APMA0340
Return to Mathematica tutorial for the second course APMA0330
Return to Mathematica tutorial for the first course APMA0340
Return to the main page for the course APMA0340
Return to the main page for the course APMA0330
Return to Part IV of the course APMA0330

Spring Problems

Sometimes it is necessary to consider the second derivative when constructing a mathematical model. Suppose that we have a mass lying on a flat, frictionless surface and that this mass is attached to one end of a spring with the other end of the spring attached to a wall. We denote the spring displacement by x. If x>0, then the spring is stretched. If x<0, the spring is compressed. If x = 0, then the spring is in a state of equilibrium. If we pull on the mass, then the mass will oscillate back and forth across the table.

spring1 = ParametricPlot[{t + Cos[3*t], 2*Sin[3*t]}, {t, 0, 4*Pi}]
wall = Graphics[{Pink, Polygon[{{-5, -5}, {-5, 5}, {-3, 5}, {-3, -5}}]}]
line1 = Graphics[Line[{{-3, 0}, {1, 0}}]]
line2 = Graphics[Line[{{4*Pi + 1, 0}, {20, 0}}]]
box = Graphics[{LightBlue, Polygon[{{20, -2.5}, {20, 2.5}, {25, 2.5}, {25, -2.5}}]}]
ground = Graphics[ Polygon[{{17, -2.5}, {26, -2.5}, {26, -3.5}, {17, -3.5}}]]
Show[box, line1, line2, wall, spring1, ground]
We can construct a differential equation that models our oscillating mass. First, we must consider the restorative force on the spring. We will make the assumption that this force depends on the displacement of the spring, F(x). Using Taylor's Theorem from calculus, we can expand F to obtain
\begin{align*} F(x) &= F(0) + F' (0) \,x + \frac{1}{2}\,x^2 F''(0) + \cdots \\ &= -k\,x + \frac{1}{2}\,x^2 F''(0) + \cdots , \end{align*}
where \( F' (0) =-k \) and F(0) = 0. If the displacement is not too large, then xn will be small for \( n \ge 2, \) and we can ignore higher ordered terms. Thus, we can consider the restorative force on the spring to be proportional to displacement of the spring from its equilibrium length,
\[ F(x) \approx - k\,x . \]
This equation is known as is Hooke's Law. We can test this law experimentally, and it is reasonably accurate if the displacement of the spring, x, is not too large.

By Newton's second law of motion, the force on the mass m must be

\[ F(x) = m\, \ddot{x} . \]
Setting the two forces equal, we have a second order differential equation,
\[ m\, \ddot{x} = -k\, x \qquad\mbox{or} \qquad \ddot{x} + \omega^2 x =0 , \]
where \( \omega^2 = k/m , \) which describes our oscillating mass. The spring-mass system is an example of a harmonic oscillator.

Now let us add a damping force to our system. For example, we might add a dashpot, a mechanical device that resists motion, to our system. Think of a dashpot as the small cylinder that keeps your screen door from slamming shut. The simplest assumption would be to take the damping force of the dashpot to be proportional to the velocity of the mass, \( \dot{x}(t) . \) Thus, we have will have an additional force, \( F(x) = -b\,\dot{x}(t) , \) acting on our mass, where b > 0. Our new equation for the spring-mass system is

\[ m\,\ddot{x} + b\,\dot{x} + k\, x =0 . \]

Consider a linear damped oscillator:
\[ m\,\ddot{x} + b\,\dot{x} + k\, x =0 ,\qquad \mbox{with initial conditions: } x(0) = x_0 , \ \dot{x}(0) = v_0 , \]
for some positive parameters m, b, and k. It is convenient to transfer the given differential equation to a system by introducing the velocity variable \( v= \dot{x} . \) The result is
\[ v= \dot{x} , \quad \dot{v} = - \frac{b}{m}\,v - \frac{k}{m}\, x , \qquad \mbox{subject to } x(0) = x_0 , \ v(0) = v_0 . \]
We define the equation for Mathematica:
osc[m_, b_, k_, x0_, v0_] :=
{x'[t] == v[t], v'[t] == -(b / m) * v[t] - (k / m) * x[t], x[0] == x0, v[0] == v0}
We define a routine which gives the result of integrating the differential equation with NDSolve, starting at t = 0, and integrating out to t = tfinal.

1.4.9. Applications

Example 4.9.2: Pure oscillations give the solutions of the following differential equation subject to the given initial conditions (displacement is 1, but velocity is zero):

\[ \ddot{x} (t) + 25\,x(t) =0 , \qquad x(0) =1, \quad \dot{x}(0) =0. \]

Clear[x,y,t,z];
soln = DSolve[{x''[t] + 25 x[t] == 0, x[0] == 1, x'[0] == 0}, x[t], t]
Plot[x[t] /. soln, {t, -1, 2.5}]
s[t_] = x[t]/.soln[[1]]
Plot[s[t],{t,0,3},AxesLabel->{"t","Displacement"}]

Example 4.9.3: Consider the initial value problem

\[ \ddot{x} + 17\,\dot{x} + 16\,x(t) =0, \qquad x(0) =1, \quad \dot{x} (0) =5. \]
soln = DSolve[{x''[t] + 17 x'[t] + 16 x[t] == 0, x[0] == 1, x'[0] == 5}, x[t], t];
s1[t_] = x[t] /. soln[[1]]
Out[32]= 1/5 E^(-16 t) (-2 + 7 E^(15 t))
Plot[s1[t], {t, 0, 3}, AxesLabel -> {"t", "Displacement"}]

When does the maximum excursion occur?

FindRoot[s1'[t] == 0, {t, 0}]
Out[9]= {t -> 0.101322}
What is the maximum excursion?
s1[t /. %]
Out[10]= 1.18603

Example 4.9.4: Now we turn to a differential equation whose characteristic polynomial has a double root:

\[ \ddot{x} + 6\,\dot{x} + 9\,x =0, \qquad x(0) =-1, \quad \dot{x}(0) =5. \]

soln = DSolve[{x''[t] + 6 x'[t] + 9 x[t] == 0, x[0] == -1, x'[0] == 5}, x[t], t]
s2[t_] = x[t] /. soln[[1]]
Plot[s2[t], {t, 0.2, 3}, PlotStyle->{Thick,Magenta}]

Out[11]= {{x[t] -> E^(-3 t) (-1 + 2 t)}}
Out[12]= E^(-3 t) (-1 + 2 t)
Out[13]=

When does the extreme occur?
Solve[s2'[t] == 0, t]
out[48]= {{t -> 5/6}}
Where is the extreme?
s2[t /. %]
out[49]= {2/(3 E^(5/2))}
N[%]
out[50]= {0.0547233}

Example 4.9.5: We consider another initial value problem

\[ \ddot{x} + 2\,\dot{x} + 37\,x =0, \qquad x(0) =-1, \quad \dot{x} (0) =5. \]
Since the characteristic polynomial has complex roots \( \lambda = -1 \pm 6{\bf j} ,\) the solutions decay with oscillations.
soln = DSolve[{x''[t] + 2 x'[t] + 37 x[t] == 0, x[0] == -1,
x'[0] == 5}, x[t], t]
Out[2]= -(1/3) E^-t (3 Cos[6 t] - 2 Sin[6 t])}}
s3[t_] = x[t] /. soln[[1]]

Clear[x,t];
soln = DSolve[{x''[t] + 2 x'[t] + 37 x[t] == 0, x[0] == -1, x'[0] == 5}, x[t], t]
s3[t_] = x[t] /. soln[[1]]
Plot[s3[t], {t, 0, 3}, PlotRange -> {-1, 1}]

Out[19]= {{x[t] -> -(1/3) E^-t (3 Cos[6 t] - 2 Sin[6 t])}}
Out[20]= -(1/3) E^-t (3 Cos[6 t] - 2 Sin[6 t])
Out[21]=

Now we define the amplitude of these oscillations:
amplitude = s3[t] /. c3_ Exp[c4_] (c1_ Cos[a_] + c2_ Sin[a_]) -> c3*Sqrt[c1^2 + c2^2]
Out[19]= -(Sqrt[13]/3)
Plot[{amplitude*Exp[-t], -amplitude*Exp[-t], s3[t]}, {t, 0, 3}, PlotStyle -> {Blue, Blue, Black}]

Example 4.9.6: Forced oscillations are modeled by the following nonhomogeneous equation:

\[ \ddot{x} + 10\,\dot{x} + 41\,x =25\,\sin 4t, \qquad x(0) =-1, \quad \dot{x}(0) =5. \]

soln = DSolve[{x''[t] + 10 x'[t] + 41 x[t] == 25 Sin[4 t], x[0] == -1, x'[0] == 5}, x[t], t]
Out[1]= {{x[t] ->
1/712 E^(-5 t) (-392 Cos[4 t] - 445 E^(5 t) Cos[4 t] +
125 E^(5 t) Cos[4 t] Cos[8 t] + 200 Sin[4 t] -
200 E^(5 t) Cos[8 t] Sin[4 t] + 200 E^(5 t) Cos[4 t] Sin[8 t] +
125 E^(5 t) Sin[4 t] Sin[8 t])}}

We extract the transient solution and the steady-state solution. First, we define the solution function:

ss[t_] = x[t] /. soln[[1]]
Out[2]= 5/712 E^(-5 t) (64 Cos[4 t] - 89 E^(5 t) Cos[4 t] + 25 E^(5 t) Cos[4 t] Cos[8 t] + 40 Sin[4 t] -
40 E^(5 t) Cos[8 t] Sin[4 t] + 40 E^(5 t) Cos[4 t] Sin[8 t] + 25 E^(5 t) Sin[4 t] Sin[8 t])
Then we find the steady-state solution by finding the general solution of the given differential equation
gensol = DSolve[{x''[t] + 10 x'[t] + 41 x[t] == 25 Sin[4 t]}, x[t], t]
Out[3]= {{x[t] -> E^(-5 t) C[2] Cos[4 t] + E^(-5 t) C[1] Sin[4 t] - 5/712 (89 Cos[4 t] - 25 Cos[4 t] Cos[8 t] +
40 Cos[8 t] Sin[4 t] - 40 Cos[4 t] Sin[8 t] - 25 Sin[4 t] Sin[8 t])}}
gs[t_] = x[t] /. gensol[[1]]
Out[4]= E^(-5 t) C[2] Cos[4 t] + E^(-5 t) C[1] Sin[4 t] - 5/712 (89 Cos[4 t] - 25 Cos[4 t] Cos[8 t] +
40 Cos[8 t] Sin[4 t] - 40 Cos[4 t] Sin[8 t] - 25 Sin[4 t] Sin[8 t])
and then vanishing arbitrary constants:
steadyState[t_] = gs[t] /. {C[1] -> 0, C[2] -> 0}
Out[5]= -(5/712) (89 Cos[4 t] - 25 Cos[4 t] Cos[8 t] + 40 Cos[8 t] Sin[4 t] - 40 Cos[4 t] Sin[8 t] - 25 Sin[4 t] Sin[8 t])
Then the transient solution will be their difference:
transient[t_] = Simplify[ss[t] - steadyState[t]]
Out[6]= 1/89 E^(-5 t) (-49 Cos[4 t] + 25 Sin[4 t])
Then we check the initial conditions:
steadyState[0]
Out[7]= -(40/89)
D[steadyState[t], t] /. t -> 0
Out[8]= 100/89
transient[0]
Out[9]= -(49/89)
D[transient[t], t] /. t -> 0
Out[10]= 345/89
ss[0]
Out[11]= -1
D[ss[t], t] /. t -> 0
Out[12]= 5

Now we consider a similar problem but with another input function:

Clear[soln,L,x]
L[t_,x_] = x''[t] + 10 x'[t] + 41 x[t];
solnRule = DSolve[{L[t, x] == 8 *Exp[-5*t] Cos[2 t], x[0] == -1, x'[0] == 5}, x[t], t]
Out[3]= {{x[t] ->
1/6 E^(-5 t) (-13 Cos[4 t] + 6 Cos[t]^2 Cos[4 t] +
Cos[4 t] Cos[6 t] + 3 Sin[2 t] Sin[4 t] + Sin[4 t] Sin[6 t])}}
s[t_] = x[t] /. solnRule[[1]]
Out[4]= 1/6 E^(-5 t) (-13 Cos[4 t] + 6 Cos[t]^2 Cos[4 t] +
Cos[4 t] Cos[6 t] + 3 Sin[2 t] Sin[4 t] + Sin[4 t] Sin[6 t])
Then we find the general solution
gensol = DSolve[{x''[t] + 10 x'[t] + 41 x[t] == 8*Exp[-5*t] Cos[2 t]}, x[t], t]
Out[5]= {{x[t] ->
E^(-5 t) C[2] Cos[4 t] + E^(-5 t) C[1] Sin[4 t] +
1/6 E^(-5 t) (6 Cos[t]^2 Cos[4 t] + Cos[4 t] Cos[6 t] +
3 Sin[2 t] Sin[4 t] + Sin[4 t] Sin[6 t])}}
gs[t_] = x[t] /. gensol[[1]]
Out[6]= E^(-5 t) C[2] Cos[4 t] + E^(-5 t) C[1] Sin[4 t] +
1/6 E^(-5 t) (6 Cos[t]^2 Cos[4 t] + Cos[4 t] Cos[6 t] +
3 Sin[2 t] Sin[4 t] + Sin[4 t] Sin[6 t])
Now the steady-state solution is
steadyState[t_] = gs[t] /. {C[1] -> 0, C[2] -> 0}
Out[6]= 1/6 E^(-5 t) (6 Cos[t]^2 Cos[4 t] + Cos[4 t] Cos[6 t] + 3 Sin[2 t] Sin[4 t] + Sin[4 t] Sin[6 t])

and the transient solution becomes

transient[t_] = Simplify[s[t] - steadyState[t]]
Out[7]= -(13/6) E^(-5 t) Cos[4 t]
Now we plot two transient solutions:
Plot[{1/89 E^(-5 t) (-49 Cos[4 t] + 25 Sin[4 t]) , transient[t]},{t,0,3}, PlotStyle->{Thick,{Blue,Pink}}]

We also plot two solutions

Plot[{ss[t], s[t]}, {t, 0, 2}, PlotStyle -> {Thick, {Cyan, Purple}}]

Transient solutions could be plotted in a similar way.

Plot[x[t] - steadyState[t], {t, 0, 1.5}, PlotRange -> {-1/2, 0.2}]

Example 4.9.7: Undamped Spring with External Forcing.

Enter the coefficients and the differential equation:

F0 = 1; omega = 2; k = 5;
ode = y''[t] + k^2*y[t] == F0*Cos[omega*t];
Solve the initial value problem:
DSolve[{ode, y[0] == 0, y'[0] == 0}, y[t], t]
Out[3]= {{y[t] ->
1/210 (-10 Cos[5 t] + 7 Cos[3 t] Cos[5 t] + 3 Cos[5 t] Cos[7 t] +
7 Sin[3 t] Sin[5 t] + 3 Sin[5 t] Sin[7 t])}}
Copy the output from the preceding step to define the solution function soln[t]:
soln[t_] = y[t] /. %[[1]]
Out[4]= 1/210 (-10 Cos[5 t] + 7 Cos[3 t] Cos[5 t] + 3 Cos[5 t] Cos[7 t] +
7 Sin[3 t] Sin[5 t] + 3 Sin[5 t] Sin[7 t])
Plot the solution function.
Plot[soln[t], {t, 0, 10}, AxesLabel -> {"t", "y[t]"}, PlotStyle->{Thick,Pink}]
Define the derivative function as v[t], and plot the trajectory in the phase plane. (You may need to enlarge the plot to see details.)

Example 4.9.8: Beats

We reconsider Example 1.4.7 for periodic input, when resonance is observed.
Enter the coefficients and the differential equation.

F0 = 1; omega = 4; k = 5;
ode = y''[t] + k^2*y[t] == F0*Cos[omega*t];
Solve the initial value problem:
DSolve[{ode, y[0] == 0, y'[0] == 0}, y[t], t]
Out[3]= {{y[t] ->
1/90 (-10 Cos[5 t] + 9 Cos[t] Cos[5 t] + Cos[5 t] Cos[9 t] +
9 Sin[t] Sin[5 t] + Sin[5 t] Sin[9 t])}}
Copy the output from the preceding step to define the solution function y[t]:

y[t_] = y[t] /. %[[1]]
Out[4]= 1/90 (-10 Cos[5 t] + 9 Cos[t] Cos[5 t] + Cos[5 t] Cos[9 t] +
9 Sin[t] Sin[5 t] + Sin[5 t] Sin[9 t])
Plot the solution function:
Plot[y[t], {t, 0, 15}, AxesLabel -> {"t", "y[t]"}, PlotStyle->{Thick,Green}]

Define the derivative function as v[t], and plot the trajectory in the phase plane:

v[t_] = D[y[t], t]
Out[6]= 1/90 (36 Cos[5 t] Sin[t] + 50 Sin[5 t] - 36 Cos[t] Sin[5 t] +
4 Cos[9 t] Sin[5 t] - 4 Cos[5 t] Sin[9 t])
ParametricPlot[{y[t], v[t]}, {t, 0, 15}, AspectRatio -> 1, AxesLabel -> {"y[t]", "v[t]"}, PlotStyle->{Thick,Gray}]

Example 4.9.9: Resonance

We repeat the commands (after quitting the kernel) here for defining the coefficients and the differential equation---now with \[omega] very close to k:

F0 = 1; omega = 4.9; k = 5;
ode = y''[t] + k^2*y[t] == F0*Cos[omega*t];
DSolve[{ode, y[0] == 0, y'[0] == 0}, y[t], t]
Out[3]= {{y[t] ->-1.0101 Cos[5. t] + 1. Cos[0.1 t] Cos[5. t] +
0.010101 Cos[5. t] Cos[9.9 t] + 1. Sin[0.1 t] Sin[5. t] +
0.010101 Sin[5. t] Sin[9.9 t] }}
y[t_] = y[t] /. %[[1]]
Out[4]= -1.0101 Cos[5. t] + 1.0101 Cos[4.9 t] Cos[5. t]^2 +
1.0101 Cos[4.9 t] Sin[5. t]^2
Plot[y[t], {t, 0, 15}, AxesLabel -> {"t", "y[t]"}, PlotStyle->{Thick,Cyan}]
Now we find the velocity
v[t_] = D[y[t], t]
Out[6]= 0. - 4.94949 Cos[5. t]^2 Sin[4.9 t] + 5.05051 Sin[5. t] -
4.94949 Sin[4.9 t] Sin[5. t]^2
and plot v versus y:
ParametricPlot[{y[t], v[t]}, {t, 0, 15}, AspectRatio -> 1, AxesLabel -> {"y[t]", "v[t]"}, PlotStyle->{Thick,Orange}]

F0 = 1; omega = 2; k = 5;
ode = y''[t] + k^2*y[t] == F0*Cos[omega*t];
DSolve[{ode, y[0] == 0, y'[0] == 0}, y[t], t]
y[t_] = y[t] /. %[[1]]
Out[3]= {{y[t] ->
1/210 (-10 Cos[5 t] + 7 Cos[3 t] Cos[5 t] + 3 Cos[5 t] Cos[7 t] +
7 Sin[3 t] Sin[5 t] + 3 Sin[5 t] Sin[7 t])}}
Out[4]= 1/210 (-10 Cos[5 t] + 7 Cos[3 t] Cos[5 t] + 3 Cos[5 t] Cos[7 t] +
7 Sin[3 t] Sin[5 t] + 3 Sin[5 t] Sin[7 t])
Plot[y[t], {t, 0, 10}, AxesLabel -> {"t", "y[t]"}, PlotStyle->{Thick,Dashed}]

v[t_] = D[y[t], t]
Out[6]= 1/210 (14 Cos[5 t] Sin[3 t] + 50 Sin[5 t] - 14 Cos[3 t] Sin[5 t] + 6 Cos[7 t] Sin[5 t] - 6 Cos[5 t] Sin[7 t])
ParametricPlot[{y[t], v[t]}, {t, 0, 10}, AspectRatio -> 1, PlotStyle->{Thick,Brown}, AxesLabel -> {"y[t]", "v[t]"}, PlotLabel -> "v[t] vesus y[t]"]

 

 

II. Aging Spring Equation

Consider the following spring equation with aging:

\[ \ddot{x} + 4\,e^{-t/4}\,x(t) =0 \qquad x(0) =1, \quad \dot{x} (0)=2 . \]

We solve it and plot in the following sequence of steps.

sol = DSolve[{x''[t] + 4*Exp[-t/4]*x[t] == 0, x[0] == 1, x'[0] == 2},
x, t] // InputForm
Out[1]=
{{x -> Function[{t}, (BesselJ[0, 16*Sqrt[E^(-t/4)]]*
BesselY[0, 16] - BesselJ[0, 16]*BesselY[0,
16*Sqrt[E^(-t/4)]] + BesselJ[1, 16]*
BesselY[0, 16*Sqrt[E^(-t/4)]] -
BesselJ[0, 16*Sqrt[E^(-t/4)]]*BesselY[1, 16])/
(BesselJ[1, 16]*BesselY[0, 16] - BesselJ[0, 16]*
BesselY[1, 16])]}}
Plot[x[t] /. sol[[1]], {t, 0, 30}]
Wrap sol =... in a bracket for clarity to obtain the same result:
(sol = DSolve[{x''[t] + 4*Exp[-t/4]*x[t] == 0, x[0] == 1, x'[0] == 2},
x, t]) // InputForm

f[t_] = x[t] /. sol[[1]]
f[0.5]

Out[5]= {1.40733}
f[0.5][[1]]
Out[6]= 1.40733
This gives 1.40733 without the curly brackets.

Another way to remove curly brackets:

f[t_] = x[t] /. sol[[1]][[1]]
Out[7]=
(BesselJ[0, 16 Sqrt[E^(-t/4)]] BesselY[0, 16] -
BesselJ[0, 16] BesselY[0, 16 Sqrt[E^(-t/4)]] +
BesselJ[1, 16] BesselY[0, 16 Sqrt[E^(-t/4)]] -
BesselJ[0, 16 Sqrt[E^(-t/4)]] BesselY[1, 16])/(BesselJ[1,
16] BesselY[0, 16] - BesselJ[0, 16] BesselY[1, 16])
f[0.5]
Out[8]= 1.40733

 

 

 

Differential Equations of higher order

Fundamental Sets of Solutions

General Solutions

Complex Roots

Reduction of order

Variation of Parameters

Method of Undetermined Coefficients

Operator Methods

Numerical Solutions

Spring Problems

Pendulum

Electric Circuits

Dyer Model

Applications

Boundary Value Problems

 

Return to Mathematica page

Return to the main page (APMA0330)
Return to the Part 1 (Plotting)
Return to the Part 2 (First Order ODEs)
Return to the Part 3 (Numerical Methods)
Return to the Part 4 (Second and Higher Order ODEs)
Return to the Part 5 (Series and Recurrences)
Return to the Part 6 (Laplace Transform)