PAGE 1

28. Differential Equation Solution

\[ u'' + 2u = 0 \quad u(0) = 0, \quad u'(0) = 2 \]
\[ u = \sqrt{2} \sin \sqrt{2}t \quad u' = 2 \cos \sqrt{2}t \]

\[ \text{period} = \frac{2\pi}{\sqrt{2}} = \sqrt{2}\pi \]

9. Damped Oscillations

DE: \[ u'' + 20u' + 196u = 0 \quad u(0) = 2, \quad u'(0) = 0 \]

\[ u = e^{-10t} \left[ 2 \cos(4\sqrt{6}t) + \frac{5}{\sqrt{6}} \sin(4\sqrt{6}t) \right] \]

\[ \omega = 4\sqrt{6} \]

\[ u = \frac{7}{\sqrt{6}} e^{-10t} \cos\left( 4\sqrt{6}t - \tan^{-1}\left(\frac{5}{2\sqrt{6}}\right) \right) \]

\( |y| < 0.05 \)

\( -0.05 = u \)

A coordinate graph showing a damped sinusoidal wave decaying between dashed exponential envelope lines.

\[ T_d = \frac{2\pi}{4\sqrt{6}} \]

\[ \frac{T_d}{T} = 1.43 \]

undamped: \[ u'' + 196u = 0 \]

\[ \omega_0 = 14 \quad T = \frac{2\pi}{14} \]

PAGE 2

3.8 Forced Vibrations

\[ mu'' + \gamma u' + ku = f(t) \]

focus on \( \gamma = 0 \) cases

\( f(t) \) is a periodic function

"forcing function"

"disturbing function"

example

\[ u'' + 16u = 2 \cos 2t \quad u(0) = u'(0) = 0 \]

\[ u = C_1 \cos 4t + C_2 \sin 4t + Y \]

\[ Y = A \cos 2t + B \sin 2t \]

undetermined coeff is good.
\[ A = \frac{1}{6} \quad B = 0 \]

then w/ ICs, \[ C_1 = -\frac{1}{6}, \quad C_2 = 0 \]

\[ u(t) = \frac{1}{6} \cos 2t - \frac{1}{6} \cos 4t \]

Same amplitude, different frequencies

PAGE 3

Transform with Identities and a "Trick"

\[ \cos At - \cos Bt \]\[ = \cos \left( \frac{At - Bt}{2} + \frac{At + Bt}{2} \right) - \cos \left( \frac{At - Bt}{2} - \frac{At + Bt}{2} \right) \]
\[ u = \frac{1}{6} \cos 2t - \frac{1}{6} \cos 4t = \frac{1}{6} (\cos 2t - \cos 4t) \]\[ = \frac{1}{6} \left[ \cos(-t + 3t) - \cos(-t - 3t) \right] \]

Then use the identity: \( \cos(A \pm B) = \cos A \cos B \mp \sin A \sin B \)

\[ = \frac{1}{6} \left[ \cos(-t) \cos(3t) - \sin(-t) \sin(3t) - (\cos(-t) \cos(3t) + \sin(-t) \sin(3t)) \right] \]\[ = \frac{1}{6} \left[ -2 \sin(-t) \sin(3t) \right] = \frac{1}{3} \sin(t) \sin(3t) \]

Note: \( -\sin(-t) = \sin(t) \)

PAGE 4
\[ u = \frac{1}{3} \sin(t) \sin(3t) \]

In this expression, \( \sin(t) \) is the slow component and \( \sin(3t) \) is the fast component.

The "Beat" Phenomenon

A "beat" occurs when two waves with nearly the same frequencies interfere with each other.

Graph of a fast sine wave oscillation contained within a slower dashed envelope of  \pm \frac{1}{3} \sin(t) .
PAGE 5

Example: Forced Oscillations and Resonance

\[ u'' + 16u = 2 \cos 4t, \quad u(0) = \frac{1}{3}, \, u'(0) = 0 \]

The general solution is of the form:

\[ u = C_1 \cos 4t + C_2 \sin 4t + At \cos 4t + Bt \sin 4t \]

Solving for the constants, we obtain the specific solution:

\[ u = \frac{1}{3} \cos 4t + \frac{1}{4}t \sin 4t \]
The term \( \frac{1}{4}t \sin 4t \) goes to \( \infty \) as \( t \to \infty \).

Resonance

"Resonance" occurs when the applied force has the same or nearly the same frequency as the fundamental frequency of the system.

Graph of u versus t showing an oscillating wave with linearly increasing amplitude over time.
Figure 1: Oscillatory motion with increasing amplitude due to resonance.
PAGE 6

Computer Project 1. Nonlinear Springs

Goal: Investigate the behavior of nonlinear springs.

Tools needed: ode45, plot

Description

For certain (nonlinear) spring-mass systems, the spring force is not given by Hooke's Law but instead satisfies:

\[ F_{\text{spring}} = ku + \epsilon u^3, \]

where \( k > 0 \) is the spring constant and \( \epsilon \) is small but may be positive or negative and represents the "strength" of the spring (\( \epsilon = 0 \) gives Hooke's Law). The spring is called a hard spring if \( \epsilon > 0 \) and a soft spring if \( \epsilon < 0 \).

Diagram of a mass-spring system attached to a wall, showing a coiled spring and a rectangular mass.
Figure 2: Nonlinear spring-mass system diagram.

Questions

Suppose a nonlinear spring-mass system satisfies the initial value problem:

\[ \begin{cases} u'' + u + \epsilon u^3 = 0 \\ u(0) = 0, \quad u'(0) = 1 \end{cases} \]

Use ode45 and plot to answer the following:

  1. Let \( \epsilon = 0.0, 0.2, 0.4, 0.6, 0.8, 1.0 \) and plot the solutions of the above initial value problem for \( 0 \le t \le 20 \). Estimate the amplitude of the spring. Experiment with various \( \epsilon > 0 \). What appears to happen to the amplitude as \( \epsilon \to \infty \)? Let \( \mu^+ \) denote the first time the mass reaches equilibrium after \( t = 0 \). Estimate \( \mu^+ \) when \( \epsilon = 0.0, 0.2, 0.4, 0.6, 0.8, 1.0 \). What appears to happen to \( \mu^+ \) as \( \epsilon \to \infty \)?
  2. Let \( \epsilon = -0.1, -0.2, -0.3, -0.4 \) and plot the solutions of the above initial value problem for \( 0 \le t \le 20 \). Estimate the amplitude of the spring. Experiment with various \( \epsilon < 0 \). What appears to happen to the amplitude as \( \epsilon \to -\infty \)? Let \( \mu^- \) denote the first time the mass reaches equilibrium after \( t = 0 \). Estimate \( \mu^- \) when \( \epsilon = -0.1, -0.2, -0.3, -0.4 \). What appears to happen to \( \mu^- \) as \( \epsilon \to -\infty \)?
  3. Given that a certain nonlinear hard spring satisfies the initial value problem:
    \[ \begin{cases} u'' + \frac{1}{5}u' + (u + \frac{1}{5}u^3) = \cos \omega t \\ u(0) = 0, \quad u'(0) = 0 \end{cases} \]
    plot the solution \( u(t) \) over the interval \( 0 \le t \le 60 \) for \( \omega = 0.5, 0.7, 1.0, 1.3, 2.0 \). Continue to experiment with various values of \( \omega \), where \( 0.5 \le \omega \le 2.0 \), to find a value \( \omega^* \) for which \( |u(t)| \) is largest over the interval \( 40 \le t \le 60 \).
PAGE 7

ode45 Differential Equation Solver

This routine uses a variable step Runge-Kutta Method to solve differential equations numerically. The syntax for ode45 for first order differential equations and that for second order differential equations are basically the same. However, the .m files are quite different.

I. First Order Equations

\[ \begin{cases} y' = f(t, y) \\ y(t_0) = y_0 \end{cases} \]
  1. Create a .m file for \( f(t, y) \) (see the tutorial on numerical methods and .m files on how to do this). Save file as, for example, yp.m.
  2. Basic syntax for ode45. At a MATLAB prompt type:
    [t, y] = ode45('yp', [t0, tf], y0);

    (your version of ode45 may not require brackets around t0, tf)

    \( yp = \) the .m file of the function \( f(t, y) \) saved as yp.m

    \( t0, tf = \) initial and terminal values of \( t \)

    \( y0 = \) initial value of \( y \) at \( t0 \)

  3. For example, to numerically solve:
    \[ \begin{cases} t^2 y' = y + 3t \\ y(1) = -2 \end{cases} \quad \text{over } 1 \le t \le 4 \]
    • Create and save the file yp.m for the function \( \frac{1}{t^2}(y + 3t) \).
    • At a MATLAB prompt type:
      [t, y] = ode45('yp', [1, 4], -2);

      (your version of ode45 may not require brackets around 1, 4)

    • To print results type: [t, y]
    • To plot results type: plot(t, y)
    • To plot results type with a '+' symbol: plot(t, y, '+')

II. Second Order Equations

\[ \begin{cases} y'' + p(t)y' + q(t)y = g(t) \\ y(t_0) = y_0, y'(t_0) = y_1 \end{cases} \]
  1. First convert 2nd order equation to an equivalent system of 1st order equations. Let \( x_1 = y, x_2 = y' \):
    \[ \begin{cases} x_1' = x_2, \\ x_2' = -q(t)x_1 - p(t)x_2 + g(t) \\ x_1(t_0) = y_0, x_2(t_0) = y_1 \end{cases} \]
  2. Create and save a .m file which will return a vector-valued function. This is a little tricky so here is a specific example.
PAGE 8
  • Suppose the system is as below and \( 0 \le t \le 4 \)
    \[ \begin{cases} x_1' = x_2, \\ x_2' = -t x_1 - e^t x_2 + 3 \sin 2t \\ x_1(0) = 2, x_2(0) = 8 \end{cases} \]
  • Create the following function file and save it as F.m:
    function xp=F(t,x)
    xp=zeros(2,1); % since output must be a column vector
    xp(1)=x(2);
    xp(2)=-t*x(1)+exp(t)*x(2)+3*sin(2*t);
    

    (Don't forget the ";" after each line.)

  • Basic syntax for ode45. At a MATLAB prompt, type:
    [t, x] = ode45('F', [t0, tf], [x10, x20]);

    \( F = \) the .m file of the vector-function saved as above

    \( t0, tf = \) initial and terminal values of \( t \)

    \( x10 = \) initial value of \( x_1 \) at \( t0 \): \( x10 = x_1(t_0) \)

    \( x20 = \) initial value of \( x_2 \) at \( t0 \): \( x20 = x_2(t_0) \)

(The example above becomes: [t, x] = ode45('F', [0, 4], [2, 8]);)

  • Since \( x_1(t) = y \), to print out the values of the solution \( y \) for \( t_0 \le t \le t_f \), at a MATLAB prompt type: [t, x(:, 1)]
  • To plot the solution on a graph \( t \) vs \( y \), type: plot(t, x(:, 1)) (since the vector \( x \) has 1st component \( x_1 = y \) and 2nd component \( x_2 = y' \).)
  • To plot \( x_1 \) vs \( x_2 \) (phase plane) type: plot(x(:, 1), x(:, 2))