Excercise 7

1. A baseball in free-fall

Continue the example of the lecture notes on the example for a baseball thrown straight up. We solved the equation
	(d^2) (y/d)(t^2) = a = 
          g,
with the initial conditions
	y(t = 0) = 0.

1.1. We were left with one constant undefined. Fix that constant by using the additional initial condition
	\(TraditionalForm\`d\(y(0)\)\/dt = \(\(y\& . \)(0) = 10  m/s\)\).

1.2. Find the highest point of the ball's path and the moment when the ball returns to the height it was thrown from.

Solutions

1.1 Simply plug into the equation in the lecture notes the equation \(TraditionalForm\`\(y\&, \)(t = 0) = 10\).

Sol11 = 
      DSolve[{y''[t] == a, y[0] == 0, y'[0] == 10}, 
        y[t], t]

1.2 Take the solution of (1.1), and define a new function to be that solution. Solve the time at which the ball reaches its maximum height, dy((Subscript[t, max]))/dt = 
        0. Plug (Subscript[t, max]) into the solution of previous part, and you get the maximum height.

y12[t_] = y[t] /. Sol11
tmax = 
      t /. Solve[D[y12[t], t] == 0, t]
y12[tmax] /. a -> -9.81
Clear[tmax]

2. A brick falling in water

Continue the example on the brick dropped from a boat. Our differential equation for the position y of the brick was
	\(TraditionalForm\`F = \(m\ g - k\ y\& .  = m\ \(y\& .. \)\)\).

2.1. Set some sensible initial conditions for the brick's position, (i) solve the equation and give the constant some numeric values. (ii) Plot the solution. What do you find? Study the effect of changing the mass of the brick without chaning the coefficent k. (iii) Find the terminal velocity of the falling brick (hint: use the Limit-command, to tell that m and k are positive, give Assumptions->{k>0, m>0} to the command).

2.2. Solve the brick example numerically. Compare the numerical solution to the analytic one by plotting both solutions, or by comparing the values of ygiven by the two solutions.

Solutions

2.1 Firstly, just solve the equations with initial conditions \(TraditionalForm\`y(0) = 0, \ \(y\& . \)(0) = 0\).

Sol21 = 
      DSolve[{m y''[t] == m g - k y'[t], y'[0] == 
              0, y[0] == 0}, y[t], t] // Simplify

For convenience, let's first make a regular function out of the rule given by DSolve. Then, plot that with the constants g,  m and k substituted with some numerical values.

y21[t_] = y[t] /. Sol21[[1]]
Plot[
      y21[t] /. {k -> 1, g -> -9.81, m -> 1}, {t, 0, 
        4}]

Finally, take the limit t -> ∞ of the derivative using the Limit-command. We need to give additional assumptions to Limit.

Limit[D[y21[t], t], t -> Infinity, 
      Assumptions -> {m > 0, k > 0}]

2.2. To get the numerical solution, one needs to first substitute in the numerical values for the parameters, and then use NDSolve in place of DSolve. This is done by inserting the /. and a rule right after the equation in the first argument to DSolve. Like in (2.1), we make a function out of the solution and plot that.

Sol22 = 
      NDSolve[{m y''[t] == m g - k y'[t], y'[0] == 
                0, y[0] == 0} /. {g -> -9.81, k -> 1, 
              m -> 1}, y[t], {t, 0, 10}] // 
        Simplify
y22[t_] = y[t] /. Sol22
Plot[y22[t], {t, 0, 10}]

3. Flying rocket

3.1. Solve the rocket example, and plot the solution, but this time take into account the force due to air-resistence, \(TraditionalForm\`F\_air = \(\(-k\)\ v\^2 = \(-k\)\ \(\(x\& . \)(t)\)\^2\)\). Use k = 0.3 in your calculations. Is the equation of motion still linear? Compare the path of this rocket to the one studied in the examples.

Solutions

3.1 No biggie, just add the given term (Subscript[F, air]) into the equations and you're all set to solve the equations using NDSolve. Both solutions, with and without friction are solved and then plotted.

Clear[m];
g = 9.81;
m0 = 100;
a = 5;
f = 3000;
k = 0.3
m[t_] = m0 - a t;
Sol31nofriction = 
      NDSolve[{f  -  m[t] g == m[t] x''[t], 
          x[0] == 0, x'[0] == 0}, 
        x[t], {t, 0, 10}]
Sol31 = 
      NDSolve[{f  -  m[t] g - k x'[t]^2 == m[t] x''[t], 
          x[0] == 0, x'[0] == 0}, x[t], {t, 0, 10}]
Plot[{x[t] /. Sol31nofriction, x[t] /. Sol31}, {t, 0, 10}, 
      PlotStyle -> {Hue[0], Hue[1/3]}]

4. Frictionless and damped oscillators

4.1 Solve the frictionless
	\(TraditionalForm\`m\ \(x\& .. \) + k\ x = 0\)
and damped 
	\(TraditionalForm\`m\ \(x\& .. \) + γ\ x\& .  + k\ x = 0\)
oscillators with the same initial conditons. Plot both of them in the same coordinates. What do you find? Does the damping affect the frequency of the oscillations?

Solutions

4.1 Trivial. Use NDSolve as usual, and plot the solutions.

m = 1;
k = 0.3;
g = 0.1;
Sol41Frictionless = 
      NDSolve[{m x''[t] + k x[t] == 0, 
          x[0] == 0, x'[0] == 1}, 
        x[t], {t, 0, 100}]
Sol41Friction = 
      NDSolve[{m x''[t] + g x'[t] + k x[t] == 0, 
          x[0] == 0, x'[0] == 1}, 
        x[t], {t, 0, 100}]
Plot[{x[t] /. Sol41Frictionless, x[t] /. Sol41Friction}, {t, 0, 
        50}, PlotStyle -> {Hue[0], Hue[1/3]}]

5. LC- and LCR-circuits

5.1. Continue the example on the LCR-circuits. The general analytic solution was given in the lecture notes. Now fix the constants and study the behaviour of the system by plotting the solutions. Vary the values for the components and see what happens.

Solutions

5.1. Just solve the LCR-circuit equations using DSolve for example, and see how the solutions behave.

Clear[l, r, c, q, t]
Sol51 = 
      DSolve[{l q''[t] + r q'[t] + q[t]/c == 0}, q[t], t] // 
        Simplify
Simplify[%, Assumptions -> {-4 l + c (r^2) > 0}]

6. Competion between populations

6.1 Study the solutions of equation dN(
              t)/dt = ϵ N(
              t) - β ((N(t))^2) by varying the values of the parameters ϵ (birth-rate) and β (shadowing), and the initial populations N(t = 0). Plot the population as a function of time.

6.2 Study the behaviour of the equations describing the life of the hares and foxes and of spruces and pines by using different values for the parameters. Plot the populations of foxes as a function of the hare population, and also the spruce population as a function of the number of spruces. Using the figures, find the stable solutions where the populations stay constant in time.

Additional excercises

You can try these if you've time and interest. These contain rather advanced stuff, so you can safely ignore these.

1. Nonlinear oscillator. Consider a body attached to a spring again. If the amplitude of the oscillations is not small, the force excerted by the spring will not follow Hooke's law, (Subscript[F, spring]) = -k x. Instead, in general higher powers of x appear. Suppose the force is given by the function (Subscript[F, spring]) = -k x - ϵ (x^3). Solve the equations of motion without friction for m = 1kg, k = 0.5N/m and
	 (i) ϵ = 0
	 (ii) ϵ = 0.1N/m\(TraditionalForm\`\(\( \)\(\^3\)\)\)
Plot the solutions of (i) and (ii) in the same figure. Study the effect of changing the initial displacement, x(t = 0).

2. Coupled oscillators. Consider two mass points of masses (Subscript[m, 1]) and (Subscript[m, 2]). The first is connected to a ceiling with a spring (Subscript[k, 1])and the second is connected to the first by a spring (Subscript[k, 2]). Let (Subscript[x, 1])and (Subscript[x, 2]) be the displacements from the equilibrium positions of the two masspoints. We can write the equations of motion
	\(TraditionalForm\`m\_1\ \(x\& .. \)\_1 = \(-k\_1\) x\_1 - \(k\_2\)(
            x\_1 - x\_2)\)
\(TraditionalForm\`\(m\_2\) \(x\& .. \)\_2 = \(-\(\(k\_2\)(
            x\_2 - x\_1)\)\)\)
Solve numerically the above equations and plot the solutions. Set initial conditions to your liking and try different values of (Subscript[m, 1]), (Subscript[m, 2]), (Subscript[k, 1]) and (Subscript[k, 2]).

Solutions

Additional.1. The equations of motion are no longer linear for ϵ > 0. For linear equations, you can always multiply the solution by a constant, and get a new solutions. Therefore the amplitude of the oscillations does not affect their frequency. If we introduce the nonlinear response from the spring, we find that large oscillations have a higher frequency than small oscillations. This behaviour is realistic for real systems, obviously one cannot make the oscillations arbitrarily large, and in some point one must consider the ϵ (x^3)term in the force of the spring. Experiment with the initial displacement!

eq := {m x''[t] + k x[t] + e x[t]^3 == 0, 
        x[0] == 1, x'[0] == 0}
m = 1;
k = 0.5;
e = 0;
SolA1a = NDSolve[eq, x[t], {t, 0, 100}]
e = 0.1;
SolA1b = NDSolve[eq, x[t], {t, 0, 100}]
Plot[{x[t] /. SolA1a, x[t] /. SolA1b}, {t, 0, 50}, 
      PlotStyle -> {Hue[0], Hue[1/3]}]

Additional.2. This turned out to be a fun excercise. First, construct the solutions.

Clear[m1, m2, k1, k2, m, k]

    
eq = {m1 x1''[t] == -k1 x1[t] - 
              k2 (x1[t] - x2[t]), 
          m2 x2''[t] == -k2 (x2[t] - x1[t]), 
          x1[0] == 1, 
          x2[0] == -1, x1'[0] == 0, x2'[0] == 
            0} /. {m1 -> 1, m2 -> 2, k1 -> 1, 
          k2 -> 1}

    
{x1n[t_], 
        x2n[t_]} = {x1[t], 
          x2[t]} /. NDSolve[
            eq, {x1[t], x2[t]}, {t, 0, 100}][[1]]

Regular plot of the displacements,red is masspoint 1 and green is 2.

Plot[{x1n[t], x2n[t]}, {t, 0, 40}, 
      PlotStyle -> {Hue[0], Hue[1/3]}]

Finally, a movie! Run the following cell, and then double-click any of the figures that appear. That should start the animation. Experiment with different values of the parameters!

Table[
      ListPlot[{{0, x1n[t]}, {0, x2n[t] - 5}}, 
        PlotStyle -> PointSize[0.1], PlotRange -> {3, -8}], {t, 
        0, 100, 0.5}]