- Details
- Parent Category: Mathematics Assignments' Solutions

# We Helped With This Differential Equations Assignment: Have A Similar One?

Category | Math |
---|---|

Subject | Differential Equations |

Difficulty | College |

Status | Solved |

More Info | Differential Equations Help |

## Short Assignment Requirements

## Assignment Description

**Lab
Project Week 7: Second order equations, systems of Equations and Runge-**

**Kutta
methods**

In the last project, you coded up Euler’s method for scalar first order equations:

(1) *x*^{0 }= *f*(*t,x*)*, x*(*t*_{0})
= *x*_{0}*.*

In class you are currently considering second order scalar equations which have the form

(2) *x*^{00 }= *f*(*t,x,x*^{0})*, x*(*t*_{0})
= *x*_{0 }and *x*^{0}(*t*_{0}) = *v*_{0}*.*

For instance, the equation for a damped harmonic oscillator looks like

So how do we fix up EM to handle the “second derivative” which appears in the above? The answer is to convert the second order scalar equation into a pair of first order equations.

**Conversion from second order
scalar equation to a first order system: **Let *v*(*t*) = *x*^{0}(*t*).
Then observe that the differential equation in (2) is the same as *v*^{0 }= *f*(*t,x,v*). So all together we have

*x*^{0 }= *v*

(3) *v*^{0 }= *f*(*t,x,v*)

*x*(*t*_{0})
= *x*_{0 }and *v*(*t*_{0}) = *v*_{0}*.*

In the above only first derivatives appear. If we put **x **= (*x,v*), **f**(*t,***x**) = (*v,f*(*t,x,v*))
and **x**_{0 }= (*x*_{0}*,v*_{0}) then we
can compress the above system to

(4) **x**^{0 }= **f**(*t,***x**)*, ***x**(*t*_{0})
= **x**_{0}*.*

This looks a lot like (1) except now **x **and **f **are vectors in **R**^{2 }instead of scalars. So now you apply
EM, but adjusted to handle vectors instead of scalars.

**EM for the system **(4)**: **Pick *h *∈ **R **kind of small. Note that in (4) you are given *t*_{0 }and **x**_{0}. For integers *j *≥ 0 put

(5) *t _{j}*

_{+1 }=

*t*+

_{j }*h*and

**x**

_{j}_{+1 }=

**x**

*+*

_{j }*h*

**f**(

*t*

_{j},**x**

*)*

_{j}*.*

Then, if *h *is small enough, one expects

**x**(*t _{j}*) ∼

**x**

_{j}.That is to say, the true solution of (4) is
approximated at *t *= *t _{j }*by

**x**

*.*

_{j}**Problem 1: ***Adjust
your EM code along the lines described above to simulate solutions of the
equation*

(6) *x*^{00 }= −*x, x*(0) = 1 *and x*^{0}(0)
= 0

*for t *∈ [0*,*20]*.
Note that the true solution of this equation is x*(*t*) = cos(*t*)*.
Run your code for h *= 1*,*1*/*4*,*1*/*16*,*1*/*64 *and compare the output with the true solution.*

What you will find is that, even
for small values of *h *the numerical method does not work particularly
well. The approximation will always grow in amplitude while the true solution
does not. The reason for this is that Euler’s Method is not really very
accurate. A better method, and one that is widely used in industry, is called
“Fourth-Order Runge-Kutta” or “RK4”.

1

2

**RK4 for the system **(4)**:**^{[1]}^{ }Pick *h *∈ **R **sort of small. Note that in (4) you are given *t*_{0 }and **x**_{0}. For integers *j *≥ 0 put

(7) *t _{j}*

_{+1 }=

*t*+

_{j }*h*and

**x**

where

(8)

Then, as in EM, we have **x**(*t _{j}*)
∼

**x**

_{j}.The error |**x**(*t _{j}*)
−

**x**

*| will be smaller, for the same*

_{j}*h*, if you use RK4 vs EM, and that is the point of the next problem.

**Problem 2: ***Implement
RK4 for *(6) *for t *∈ [0*,*20]*. Run your code for h *= 1*,*1*/*4*,*1*/*16*,*1*/*64 *and compare the output with the true solution and the output from EM (with
the same step size). You will see that RK4 is remarkably better even at very
large step sizes.*

At this point you now have a RK4 solver and the ability to apply it to second order scalar equations. So let’s do it.

**Problem 3: **Implement RK4 for the “ideal
pendulum equation”

*x*^{00 }= −sin(*x*)*, x*(0) = *x*_{0 }and *x*^{0}(0)
= *v*_{0}*.*

This is the equation of motion for a rigid pendulum
of length 1, with mass equal to 1, in a place where the acceleration due to
gravity is exactly 1 and there is no friction. Here *x *is the angle the
pendulum makes with respect to vertical. Your TA will draw a picture for you.

Use your solver to simulate the above with the following initial conditions. For each, describe in words (to each other and to your TA) what is going on with pendulum.

• *x*_{0 }= 0, *v*_{0 }= 0. • *x*_{0 }= 0, *v*_{0 }= 1.

• *x*_{0 }= 0, *v*_{0 }= 2.

• *x*_{0 }= 0, *v*_{0 }= 2*.*1.

• *x*_{0 }= *π*, *v*_{0 }= 0.

• *x*_{0 }= *π*, *v*_{0 }= 0*.*1.

Use a step size of *h *= 1*/*128 and
simulate each on the interval *t *= [0*,*50].

[1] It might seem that RK4 is pretty strange; why it works is beyond the scope of this class but probably not beyond your understanding. It has something to do with Simpson’s rule for approximating integrals: you should read about it if you are curious. The Wikipedia page on Runge-Kutta is a pretty good reference.

## Assignment Description

**Lab Project 2 (Week 8): Resonance and the
Tacoma Narrows Bridge**^{[1]}

You have been extensively studying driven simple harmonic oscillators, which are modeled by the equation

(1) *mx*^{00 }+ *cx*^{0 }+ *kx *= *F*(*t*)*.*

Here, *x *is
the dispacement of the oscillator from equilbrium, *m *is the mass, *c *is
the damping constant, *k *is the spring constant and *F*(*t*) is
the driving force. This project is about **resonance**: at **resonant
frequencies**, small periodic driving forces have the ability to produce
large amplitude oscillations, due to the storage of vibrational energy.

No engineering student gets through their life without being told that resonance is why the Tacoma Narrows Bridge (TNB) fell over. If you don’t know what the TNB is, then get on your phone right now, go to YouTube and type “Tacoma Narrows Bridge” into the search bar and watch some videos. In this project, you are going to do some numerical experiments for a (slightly!) more sophisticated model of the TNB.

Here’s the standard story about the TNB: wind blew across
the bridge and the frequency of the oscillations of the wind was close to the
natural frequency of the bridge. Hence there was resonance, which caused the
oscillations of the bridge to get huge which tore it apart. It’s a good story,
and not totally wrong. But, obviously, a bridge is not a mass on a spring. So
here we’ll study a *nonlinear *modification of (1), developed by Lazer
& McKenna in 1990,^{[2]}which some people feel gives an at least
slightly more realistic model for TNB.^{[3]}

Here is the idea: the cables on a suspension bridge only exert a force when they are stretched, not when they are slack. The spring/mass system modeled by (1) has a spring which pushes back when it is compressed and pulls when stretched. To “fix up” (1) to be more like a suspension bridge, Lazer & McKenna proposed using the following equation:

(2) *y*^{00 }+ *by*^{0 }+ *k*_{1}*y *+ *k*_{2}*H*(*y*)
= 10 + *F*(*t*)*.*

In the above *y *is the
horizontal displacement of the bridge, measured “downward.” That is to say, if *y
> *0 then the bridge has moved closer to the ground. Note that we are
treating the whole bridge as one piece, a pretty big assumption! The constant *b *is the damping constant. The number “10” is meant to be the acceleration
due to gravity.^{[4]}^{ }The constants *k*_{1 }and *k*_{2 }are spring constants. Again, *F*(*t*)
is a forcing function. Lastly,

*.*

Note that *H *is
only “on” if *y *≥ 0 (which is “down”). If *y < *0, which is to
say if the bridge is elevated, then it exerts no force. Which is to say, it
models slightly more realistically how the cables would exert force.

2

It turns out this model has an interesting feature which is that there are large amplitude steady state solutions even when the driving force is not at the resonant frequency. And you will investigate this now.

**Problem 1: **Modify the RK4 code you wrote in last
project to simulate second order systems so that it works for (2). Keep it
flexible enough so that you can easily change the constants *k*_{1}, *k*_{2 }and *b*, the initial conditions and the forcing
function *F*(*t*). You might find the “heaviside” function to be useful
in getting MATLAB to implement the function *H*(*y*). (Use a step
size *h *= 1*/*16 throughout this project.)

**Problem 2: **Put *k*_{2 }= 0, *k*_{1 }= 17 and *b *= 0*.*01. Note that here
the “nonlinear part” is zeroed out, so this is really just a driven simple
harmonic oscillator. Suppose that *F*(*t*) = 0*.*1cos(*ωt*).
Use the approach of Section 4.2 in the textbook to find the general solution,
and describe its steady- state behavior (see the definition in Section 4.4).

The presence of the “10” on the
right hand side of (2) will slightly complicate things, but a simple trick
reduces the problem to exactly the form in the book. The trick is this: let *y *= *C *+ *z *where *z *is a new function and *C *is
cleverly chosen constant. Your TA should be able to help you with figure out
this special constant! It has something to do with

“equilbrium.”

**Problem 3: **Put *k*_{2 }=
0, *k*_{1 }= 17 and *b *= 0*.*01. Put *y*(0) = 0
and *y*^{0}(0) = 0. Let

*F*(*t*)
= 0*.*1cos(*ωt*). For *ω *∈ [3*,*5] (with increments of 0*.*1),
run your solver from Problem 1 from *t *= 0 to *t *= 1500. This
should be enough time for the solution to settle into its “steady-periodic
state.” Compare your numerically computed solution with the formula for the
steady state you got in Problem 2. They should be in pretty close agreement
when *t *is large (but necessarily near *t *= 0). At what value of *ω *is the amplitude of the steady state biggest? This is, of course, the
resonant frequency!

**Problem 4: **Put *k*_{2 }= 0, *k*_{1 }= 17 and *b *= 0*.*01. Let *F*(*t*)
= 0*.*1cos(4*t*). You can determine from the formulas in 2.6 choices
for *y*(0) and *y*^{0}(0) so that the “transient solution” is
exactly zero. Which is to say, you won’t need to “wait” to see that steady
state part. Use these special initial conditions in the RK4 solver. Run your
solver and see if your RK4 method with these special initial conditions and see
if it agrees with the analytic solution (it should agree almost exactly
instantly).

**Problem 5: **Take your
special initial conditions from the last problem and use them as the initial
conditions for the system with *k*_{2 }= 5, *k*_{1 }=
12, *b *= 0*.*01 and *F*(*t*) = 0*.*1cos(4*t*).
If you have done things correctly, NOTHING should change from the previous
problem. Why? This solution is called a “small amplitude periodic solution.”

**Problem 6: **With the
constants and *F*(*t*) as in Problem 5, vary *y*(0) between −6
and 0 (with increments of 0*.*1) and *y*^{0}(0) between 0 and
1 (again with increments of 0*.*1). Let your solver run from *t *= 0
to *t *= 1500 and record the amplitude of the solution at that time. Some
of the initial conditions will have reach a steady state which is exactly the
same small amplitude solution you saw in Problem 4. Other ones will be much
larger in amplitude: these are “large amplitude steady states” and a nonlinear
phenomenon, not a linear one. (That is to say, if you put *k*_{2 }=
0 there is only the small amplitude solution.) Make a plot of the resulting
steady state amplitude vs *y*(0) and *y*^{0}(0). It is
conjectured that that these sort of “large amplitude nonlinear periodic”
solutions are the real culprit behind the TNB disaster! The point is: if the
initial condition is set is just so, the solution is attracted to a

3

very large frequency oscillation even when the driving frequency is not right on the resonant frequency. NONLINEARITY = DOOM.

[1] This project is largely based on—and in some parts wholly borrowed from—Chapter
4.5 of *Differential Equations *by Blanchard, Devaney and Hall.

[2] Lazer
& McKenna, “Large-amplitude periodic oscillations in suspension bridges:
some new connections with nonlinear analysis” *SIAM Review*, 1990 pp
537-587

[3] Clearly civil engineers do way more sophisticated stuff than this when they study bridges! So if you really want a deep understand of TNB, take a class in that department!

[4] Mathematicians are lazy and can’t be bothered getting accurate physical constants; of course engineers know that acceleration due to gravity is really 9.807. But 10 is close to 9.807 so we’ll stick with it.