Corey Reed
2-6-98
Math 118A

Homework #4 (Mathematica section)

3)

MatrixForm[xvector  =  {x[t], y[t]}]

MatrixForm[v = {y[t], -c*y[t] + x[t]*(1 - x[t])}]

Solve[v == 0]

So there are two critical points, one at the origin and one at (1,0). This implies that in order to meet the special properties of x[t], the origin must be unstable, the point (1,0) must be stable, and that the phase flow from the origin to the point (1,0) cannot cross the x = 0 axis (if it did, a particle travelling along the phase flow would have a negative speed at some time, which is not what we want).

MatrixForm[
      a = {{D[v[[1]], x[t]], D[v[[1]], y[t]]}, {D[v[[2]], x[t]], 
            D[v[[2]], y[t]]}}]

At the origin,

Eigenvalues[a /. x[t] -> 0]

Note that the eigenvalues can never be imaginary. If c < 0, the both eigenvalues will be positive. If c ≥ 0, then there will be one positive and one negative eigenvalue. In either case, the origin will be unstable in at least one direction.

At the point (1,0),

Eigenvalues[a /. x[t] -> 1]

If c is less than 2 in absolute value, the the eigenvalues here will be a pair of complex conjugates. But in that case, the phase flow would spiral around the point (1,0), which would cause the phase flow to cross the x axis eventually, which should not happen. If c ≥ 2, then both eigenvalues will be negative, which will be stable. If c ≤ -2, then both eigenvalues will be positive, which will be unstable. Thus, c must be greater than or equal to 2 in order to make (1,0) a stable critical point satisfying the desired conditions.

<< Graphics`PlotField`

PlotVectorField[
      (a /. c -> 2) . xvector, {x[t], - .2, 1}, {y[t], - .2,  .2}]

So the phase flow for c ≥ 2 fits the requirements.

PlotVectorField[
      (a /. c -> 0.5) . xvector, {x[t], - .2, 1}, {y[t], 
        - .2,  .2}]

When c is between -2 and 2, the phase flow spirals, causing some orbits to cross the x axis, which does not satisfy the requirement that the speed of the particle be positive for all time.

PlotVectorField[
      (a /. c -> -2) . xvector, {x[t], - .2, 1}, {y[t], 
        - .2,  .2}]

If c is less than or equal to -2, then the point (1,0) is not stable, which does not fit the requirement.

4.a)

MatrixForm[v = {y[t], -x[t] - k*y[t]}]

Solve[v == 0]

So the only critical point is at the origin.

MatrixForm[
      a = {{D[v[[1]], x[t]], D[v[[1]], y[t]]}, {D[v[[2]], x[t]], 
            D[v[[2]], y[t]]}}]

Eigenvalues[a]

So for values of k between -2 and 2, the eigenvalues will be a pair of complex conjugates. For values of k which are greater than 2 in absolute value, the eigenvalues will both be real. When k is equal to 2 or -2, the eigenvalues will be the same.

4.b)

MatrixForm[aeq = a /. k -> -2]

Eigenvalues[aeq]

MatrixForm/@({s, j} = JordanDecomposition[aeq])

MatrixForm[si = Inverse[s]]

By splitting [ j ] up into the sum of a diagonal matrix and an upper triangle matrix, we can see that raising e to the [ j ] * t  power will give: (see written portion)

MatrixForm[ej = {{Exp[t], t*Exp[t]}, {0, Exp[t]}}]

soln = (s . ej . si)[[1]] // MatrixForm

So in this case, the solution is a linear combination of exp(t) - t * exp(t) and t * exp(t) . The "position" (x[t]) grows exponentially in time becuase the "damping" term is negative, which means it will push the particle away from the origin, instead of pushing it toward the origin.

Phase plane for one real eigenvalue. k=-2 => λ=1,1:

PlotVectorField[aeq . xvector, {x[t], -1, 1}, {y[t], -1, 1}]

MatrixForm[aeq = a /. k -> 1]

eigs = Eigenvalues[aeq]

MatrixForm[c = Eigenvectors[aeq]]

MatrixForm[ci = Inverse[c]]

MatrixForm[l = {{eigs[[1]], 0}, {0, eigs[[2]]}}]

Because the matrix [ l ] is a diagonal matrix, we can read off  e to the [ l ] * t  power:

MatrixForm[el = {{Exp[l[[1, 1]]*t], 0}, {0, Exp[l[[2, 2]]*t]}}]

MatrixForm[soln = Simplify[(c . el . ci)[[1]]]]

So the solution is a linear combination of the entries of this matrix, which is an underdamped oscillator.

Phase plane for a pair of complex conjugate eigenvalues with real and imaginary parts. Here, k=1:

PlotVectorField[
      (a /. k -> 1) . xvector, {x[t], -1, 1}, {y[t], -1, 1}]

Here is a sample solution (note that it starts at one):

w[t_] = 
      ComplexExpand[
        Re[Simplify[
            ExpToTrig[
              (((((-1))^(1/3)) 
                    (-(E^(-(((-1))^(1/3)) t)) + 
                        (E^((((-1))^(2/3)) t))))/(1 + 
                    (((-1))^(1/3)))) + 
                (((((-1))^(1/3)) 
                      (E^(-(((-1))^(1/3)) t)) + 
                    (E^((((-1))^(2/3)) t)))/(1 + 
                    (((-1))^(1/3))))]]]]

Plot[w[t], {t, 0, 10}]

MatrixForm[aeq = a /. k -> 3]

eigs = Eigenvalues[aeq]

MatrixForm[c = Eigenvectors[aeq]]

MatrixForm[ci = Inverse[c]]

MatrixForm[l = {{eigs[[1]], 0}, {0, eigs[[2]]}}]

Because the matrix [ l ] is a diagonal matrix, we can read off  e to the [ l ] * t  power:

MatrixForm[el = {{Exp[l[[1, 1]]*t], 0}, {0, Exp[l[[2, 2]]*t]}}]

MatrixForm[soln = Simplify[(c . el . ci)[[1]]]]

So the solution is a linear combination of the entries of this matrix, which is an overdamped oscillator, and the particles will not oscillate, but will fall into the origin.

Phase plane for two distinct real eigenvalues. Here, k=3:

PlotVectorField[
      (a /. k -> 3) . xvector, {x[t], -1, 1}, {y[t], -1, 1}]