Excercise sheet 7

Competion between populations

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.

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 pines. Using the figures, find the stable solutions where the populations stay constant in time.

Solutions

Here is demonstrated the effect of changing the environmental limitation factor β.

Clear[n, 
        t, ϵ, β];

    ϵ = Table[i, {i, 0.1, 1, 0.1}];

    β = Table[0.05, {i, 1, 10, 1}];

    Neb = Length[ϵ];

    Table[
        Part[NDSolve[{n'[
                  t] == ϵ[[i]] n[
                      t] - β[[i]] n[t]^2, n[0] == 10}, 
            n[t], {t, 0, 10}], 1, 1], {i, 1, Neb}];

    Plot[Evaluate[Table[n[t] /. %[[i]], {i, 1, Neb}]], {t, 0, 10}, 
      PlotStyle -> Table[Hue[(i - 1)/Neb], {i, 1, Neb}], 
      PlotRange -> All]

You can find a steady state solution by setting the right hand side of the equations to zero at time t = 0:
	((Subscript[ϵ, h]) - (Subscript[β, h]) (Subscript[N, f])(
                  0)) (Subscript[N, h])(0) = 0
(-(Subscript[ϵ, f]) + (Subscript[β, f]) (Subscript[N, h])(
                  0)) (Subscript[N, f])(0) = 0
This can be solved for (Subscript[N, f])(0) and (Subscript[N, h])(0) to obtain the stability conditions
	(Subscript[N, h])(0) = ((Subscript[ϵ, f])/(Subscript[β, f])) and (Subscript[N, f])(0) = ((Subscript[ϵ, h])/(Subscript[β, h]))
Below is an example of a steady solution. A variation of the parameters leads to oscillatory behaviour. The function Unset can be used to clear variables which have super- or subscripts (regular Clear doesn't work).

deq1 = (Subscript[N, h])'[
            t] == ((Subscript[ϵ, h]) - (Subscript[β, h]) (Subscript[N, f])[t]) 
            (Subscript[N, h])[t];

    deq2 = (Subscript[N, f])'[
            t] == (-(Subscript[ϵ, f]) + (Subscript[β, f]) (Subscript[N, h])[t]) 
            (Subscript[N, f])[t];

    (Subscript[ϵ, h]) = 2; (Subscript[ϵ, f]) = 0.5; (Subscript[β, h]) = 0.2; (Subscript[β, f]) = 
      0.05;

    NDSolve[{deq1, deq2, (Subscript[N, h])[0] == 10, 
          (Subscript[N, f])[0] == 10}, {(Subscript[N, h])[t], (Subscript[N, f])[t]}, {t, 0, 
          30}];

    Plot[{(Subscript[N, h])[t] /. %, (Subscript[N, f])[t] /. %}, {t, 0, 30}, 
      PlotStyle -> {Hue[0], Hue[0.7]}, 
      PlotRange -> All]

    Unset[(Subscript[ϵ, h])]; Unset[(Subscript[ϵ, f])];

    Unset[(Subscript[β, h])]; Unset[(Subscript[β, f])];

Again the stable situation can be analytically discovered to be
	\(TraditionalForm\`\(N\_p\)(
          0) = \(\(β\_p\) ϵ\_s - \(ϵ\_p\) γ\_s\)\/\(\(-β\_p\) β\_s + \(γ\_p\) γ\_s\)\), and \(TraditionalForm\`\(N\_s\)(
          0) = \(\(β\_s\) ϵ\_p - \(ϵ\_s\) γ\_p\)\/\(\(-β\_p\) β\_s + \(γ\_p\) γ\_s\)\).
A stable solution could be obtained if one of the above would be zero, but that's not so interesting. Note also that for some values of the parameters stable situations are impossible (the initial sizes of the populations would have to assume negative values).

deq1 = (Subscript[N, s])'[
            t] == ((Subscript[ϵ, s]) - (Subscript[γ, s]) 
                  (Subscript[N, s])[t] - (Subscript[β, s]) (Subscript[N, p])[t]) 
            (Subscript[N, s])[t];

    deq2 = (Subscript[N, p])'[
            t] == ((Subscript[ϵ, p]) - (Subscript[γ, p]) 
                  (Subscript[N, p])[t] - (Subscript[β, p]) (Subscript[N, s])[t]) 
            (Subscript[N, p])[t];

    (Subscript[ϵ, s]) = 0.2; (Subscript[ϵ, p]) = 0.12; (Subscript[β, s]) = 0.1; (Subscript[β, p]) = 
      0.1; (Subscript[γ, s]) = 0.05; (Subscript[γ, p]) = 0.02;

    NDSolve[{deq1, deq2, (Subscript[N, s])[0] == 1.33, 
          (Subscript[N, p])[0] == 1.33}, {(Subscript[N, s])[t], (Subscript[N, p])[t]}, {t, 0, 
          80}];

    Plot[{(Subscript[N, s])[t] /. %, (Subscript[N, p])[t] /. %}, {t, 0, 80}, 
      PlotStyle -> {Hue[0], Hue[0.7]}, 
      PlotRange -> All]

    Unset[(Subscript[ϵ, s])]; Unset[(Subscript[ϵ, p])];

    Unset[(Subscript[β, s])]; Unset[(Subscript[β, p])];

    Unset[(Subscript[γ, s])]; Unset[(Subscript[γ, p])];

More complicated oscillators

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

For small initial displacements, the anharmonic oscillator behaves like the harmonic one, but the effect of the third order term can be seen by increasing s. With the third order term accounted for, the force takes on larger values and the oscillation speeds up.

deq := 
        m x''[t] + k x[t] + ϵ x[t]^3 == 
          0;

    k = 0.5; s = 2; m = 0.5; T = 20;

    ϵ = 0.0;

    sol1 = 
        NDSolve[{deq, x[0] == s, x'[0] == 0}, 
          x[t], {t, 0, T}];

    ϵ = 0.1;

    sol2 = 
        NDSolve[{deq, x[0] == s, x'[0] == 0}, 
          x[t], {t, 0, T}];

    Plot[{x[t] /. sol1, x[t] /. sol2}, {t, 0, T}, 
      PlotStyle -> {Hue[0], Hue[0.6]}, 
      PlotPoints -> 100]

Note that for ParametricPlot, one layer of the solution list has to be truncated (flattened) using the Part function.

deq1 = (Subscript[m, 1]) (Subscript[x, 1])''[
              t] == -(Subscript[k, 1]) 
              (Subscript[x, 1])[t] - (Subscript[k, 2]) ((Subscript[x, 1])[t] - 
                  (Subscript[x, 2])[t]);

    deq2 = (Subscript[m, 2]) (Subscript[x, 2])''[
              t] == -(Subscript[k, 2]) ((Subscript[x, 2])[t] - 
                (Subscript[x, 1])[t]);

    (Subscript[s, 1]) = 0.1; (Subscript[s, 2]) = -0.2; (Subscript[m, 1]) = 0.1; (Subscript[m, 2]) = 0.05; (Subscript[k, 1]) = 1; 
    (Subscript[k, 2]) = 0.5; T = 20;

    sol = 
        Part[NDSolve[{deq1, deq2, (Subscript[x, 1])[0] == (Subscript[s, 1]), 
              (Subscript[x, 2])[0] == (Subscript[s, 2]), (Subscript[x, 1])'[0] == 
                0, (Subscript[x, 2])'[0] == 0}, {(Subscript[x, 1])[t], (Subscript[x, 2])[t]}, {t, 0, 20}], 
          1];

    Plot[{(Subscript[x, 1])[t] /. sol, (Subscript[x, 2])[t] /. sol}, {t, 0, T}]

    ParametricPlot[{(Subscript[x, 1])[t] /. sol, (Subscript[x, 2])[t] /. sol}, {t, 0, T}, 
      PlotPoints -> 300]

Piecewise functions, conditional evaluation

1. Define a function f as
	f(x)= {({{0,  x ≤ 0}, {1,  x > 0}})
Use both the /; method and the If command.

2. Define a piecewise function that gets as its value (x^2) when x < 3, x + 5 when 3 ≤ x < 7 and 14 when x ≥ 7. Plot the function f. How could you redefine the function so that it would be continuous?

Solutions

I will define functions f and g differently only in order to tell them apart when plotting them.

Clear[f, g, 
        x];

    f[x_] := 0 /; x ≤ 0

    f[x_] := 1 /; x > 0

    g[x_] = If[x ≤ 0, 1, 2];

    Plot[{f[x], g[x]}, {x, -1, 1}, 
      PlotStyle -> {Hue[0], Hue[2/3]}]

In the next excercise, f could be made continuous in a number of ways. They all would call for the redefinition of f in at least two intervals. Figure them out if you wish.

Clear[f, x]

    f[x_] := x^2 /; x < 3

    f[x_] := x + 5 /; 3 ≤ x < 7

    f[x_] := 14 /; x ≥ 7

    Plot[f[x], {x, 1, 9}]

Programming

1. Write a program that defines a function f(x) = sin(10 cos(x)) and plots it.

2. Make a program that solves the equation of a damped harmonic oscillator
	\(TraditionalForm\`m\ \(x\& .. \) + γ\ x\& .  + k\ x = 0, \ 
      x(0) = s, \ \(x\& . \)(0) = 0\).
Write the program into a separate file. The program should set values for m, γ, k and s, then solve the equation and finally plot the solution.

Solutions

Files can be input into Mathematica using
	SetDirectory["c://path//to//files"]
	<< "file.m"
Note that you can't use the referring technique %n when giving input from a file. Also, the outputs of substitution expressions are automatically suppressed. Otherwise, you can pretty much write ordinary Mathematica code in your file.

Clear[f, 
        x];

    f[x_] = Sin[10  Cos[x]];

    Plot[f[x], {x, -Pi, Pi}]

Clear[m, x, t, g, k, 
        s];

    m = 1.0; g = 0.1; k = 1.5; s = 0.1;

    sol = 
        Part[DSolve[{m x''[t] + g x'[t] + k x[t] == 0, 
              x[0] == s, x'[0] == 0}, x[t], t], 
          1];

    Plot[x[t] /. sol, {t, 0, 20}]