. Differential Equations

Paul Abbott

Pre-requisites

To improve the appearance of graphics, TraditionalForm will be used as the default format type for text that appears in graphics,

$FormatType = TraditionalForm;

and the default style for text in graphics is set as follows.

$TextStyle = {FontSize -> 10};

Unprotect[Equal];

listableQ(f_) := 
      MemberQ(Attributes(f), Listable)

\(TraditionalForm\`Equal /: \ 
      lhs : f_Symbol?
            listableQ[___, \ _Equal, \ ___]\  := 
\ \ \ \ \ \ \ \ Thread[\ 
        Unevaluated[lhs], \ Equal\ ]\)

Protect[Equal];

Now, we can directly manipulate equations.

In this section you will explore the trajectories of a planet moving in the gravitational field of two fixed, equal mass suns. The suns are of mass M and are located at the points +D and -D along the x-axis. Rescaling the variables by x -> D x, y -> D y, t -> t (Sqrt[ (D^3)/(G M)]), the equations become dimensionless:

x''(t)==\(\(1 - x(t)\)\/\((\((x(t) - 1)\)\^2 + \(y(t)\)\^2)\)\^\(3/2\) - \(x(t) + 1\)\/\((\((x(t) + 1)\)\^2 + \(y(t)\)\^2)\)\^\(3/2\)\),
y''(t)==\(\(-\(\(y(
                    t)\)\/\((\((x(t) - 1)\)\^2 + \(y(t)\)\^2)\)\^\(3/2\)\)\) - \(\(\(y(t)\)\/\((\((x(t) + 1)\)\^2 + \(y(t)\)\^2)\)\^\(3/2\)\)\(.\)\)\)

These coupled differential equations are very difficult to solve analytically.  In fact, the general 3-body problem in classical mechanics (where all the masses can be different and none are fixed) is a famous problem. Series solution reveals that the solution is unstable to small perturbations and the system dynamics are chaotic.

eqns={x''(t)==((1 - x(t))/((((((x(t) - 1))^2) + ((y(t))^2)))^(3/2))) - ((x(t) + 1)/((((((x(t) + 1))^2) + ((y(t))^2)))^(3/2))),y''(t)==-((y(
                          t))/((((((x(t) - 1))^2) + ((y(t))^2)))^(3/2))) - ((y(t))/((((((x(t) + 1))^2) + ((y(t))^2)))^(3/2)))};

and initial conditions:

ics={x(0) == 0,y(0) == 0,x'(0)==1,y'(0)==1};

(i.e., starting from the origin with a specified velocity), here is a parametric plot of the solution for 8 time units:

nsol = 
      NDSolve[Join[eqns, ics], {x, y}, {t, 0, 8}];

SetOptions[ParametricPlot, 
        Epilog -> {PointSize[0.02], Hue[1], Point[{-1, 0}], 
            Point[{1, 0}]}];

ParametricPlot[
        Evaluate[{x(t), y(t)} /. nsol], {t, 0, 
          8}];

It is easy to turn these steps into a procedure which can be used to investigate the trajectory as a function of the initial velocity, {(Subscript[v, x]), (Subscript[v, y])}, and total time, τ:

TrajectoryPlot(vx_, vy_, τ_:  25):=Module[{ics, nsol},ics={x(0) == 0,y(0) == 0,x'(0)==vx,y'(0)==vy};nsol = 
                NDSolve[Join[eqns, ics], {x, y}, {t, 0, τ}];ParametricPlot[
                Evaluate[{x(t), y(t)} /. nsol], {t, 
                  0, τ}];]

. Heat Equation

It is well known that  diffusion  processes cause ‘smoothing’. For example in the thermal context one would expect any steepness in the temperature profile to be smoothed out as time goes on. In this section we will consider the scaled heat equation,

((∂Q)/(∂t))=(((∂^2) Q)/(∂(x^2)))·

. Normal Modes of Coupled Oscillators

In this section you will explore the natural frequencies and normal modes of a system of N coupled oscillators.  The case with N = 3 looks like:

The parameters (Subscript[k, i]) are the spring constants; (Subscript[m, i]) are the masses of the blocks; and (Subscript[x, i])(t) are the displacements from the equilibrium positions.

.. Equations of Motion

Fixed end-points

In this section we focus on the case where the end-points are fixed, i.e., the boundary conditions are (Subscript[x, 0])(t) == (Subscript[x, N + 1])(t) == 
        0. For a system of N coupled oscillators, the equation of motion for (Subscript[m, 1])reads:

(Subscript[m, 1]) (Subscript[x, 1]^(′′))(t)==(Subscript[k, 2]) ((Subscript[x, 2])(t) - (Subscript[x, 1])(t)) - 
            (Subscript[k, 1]) (Subscript[x, 1])(t);

Symbolic analysis of the case with general values for (Subscript[k, i]) and (Subscript[m, i]) is quite complicated.  We will firstly examine the case where all the (Subscript[k, i]) are equal to k and all the (Subscript[m, i]) equal m. This case can be obtained using a substitution rule, e.g.,

% /. {(Subscript[m, n])_ -> m, 
          (Subscript[k, n])_ -> k} // Simplify

In the case of equal masses and spring constants, introducing the frequency (time-scale), ω = (Sqrt[k/m]), (i.e., rescaling time by this factor) the equations can be cast into the dimensionless form:

(Subscript[x, i]^(′′))(t) ==(Subscript[x, i + 1])(t) - 
            2  (Subscript[x, i])(t) + (Subscript[x, i - 1])(t),1 ≤ i ≤ N.

For the following we set

𝒩 = 4;

i.e., 4 masses. The fixed end-point boundary conditions are

(Subscript[x, 0])(t_) = 0;

and

(Subscript[x, 𝒩 + 1])(t_) = 0;

and the equations of motion for the N masses are

(Subscript[ℰ, n])_:=(Subscript[x, n]^(′′))(t)==(Subscript[x, n - 1])(t) + (Subscript[x, n + 1])(t) - 
              2 (Subscript[x, n])(t)/;1 ≤ n ≤ 𝒩

Matrix form of the equations

Since the motion of each mass is expected to be simple-harmonic, one can propose that the solution is of the form (this is called an ansätz):

(Subscript[x, i])_(t_) = 
      (Subscript[c, i]) cos(Ω t);

(({{2, -1, 0, 0}, {\[Placeholder], \[Placeholder], \[Placeholder], \[Placeholder]}, {\[Placeholder], \[Placeholder], \[Placeholder], \[Placeholder]}, {\[Placeholder], \[Placeholder], \[Placeholder], \[Placeholder]}})) (({{(Subscript[c, 1])}, {(Subscript[c, 2])}, {\[Placeholder]}, {\[Placeholder]}}))==\[Placeholder](({{(Subscript[c, 1])}, {\[Placeholder]}, {\[Placeholder]}, {\[Placeholder]}}))

The matrix on the left-hand side is called the normal-mode matrix. Recall that the eigenvalues, λ, and eigenvectors, x, of a matrix 𝒜 satisfy 𝒜 . x = λ x.

Remember that the eigenvalues, λ, of the matrix 𝒜 satisfy the characteristic equation, Det[𝒜 - λ I] == 
        0.

You might want to visit http://library.wolfram.com/examples/quintic/ and http://mathworld.wolfram.com/QuinticEquation.html.

General normal-mode matrix

By looking at (Subscript[ℳ, 4]) you should be able to deduce the general normal-mode matrix for a system of 𝒩 masses.

(Subscript[a, i_, j_]) := 2 /; i == j

(Subscript[ℳ, n])_ := 
      Table[(Subscript[a, i, j]),  ... ]

In general, the output from Eigensystem is not ordered. However, in most physical systems the ordering is relevant.

Let us find a method for sorting a pair of related lists:

.. Graphical Display of the Eigenvectors and Eigenvalues

ModePlot

The following code uses ListPlot to display the n-th eigenvector, 𝓊[[n]]:

ModePlot(n_, opts___) := 
      ListPlot[Join[{0}, 𝓊[[
            n]], {0}], PlotJoined -> True, 
        PlotLabel -> "Ω=" <> 
            ToString[(Sqrt[λ])[[n]]], 
        opts];

Note that the endpoints, (Subscript[x, 0]) = 0 and (Subscript[x, N + 1]) = 0 have been joined to the n-th eigenvector using the Join command and also that the (normal-mode) oscillation frequencies are computed from the (square root) of the eigenvalues.

Frequency-level diagram

energylevel[y_]:=Line[(({{0, y}, {1, y}}))]

Show[Graphics[energylevel /@ ν], 
        Axes -> True, AspectRatio -> 3, 
        Ticks -> {None, Automatic}];

General formula

It can be shown that, for the special case of equal masses and spring constants with both ends fixed, the general formula for the frequencies (i.e., the square root of the eigenvalues) satisfy (Subscript[Ω, i]) = 
        2 sin(((i π)/(2 (𝒩 + 1)))) where (Subscript[Ω, i]) is the frequency of the i-th mode for a system of N masses.

.. One End Free

.. General Tri-diagonal Matrix

In the analysis of this problem, tri-diagonal matrices have appeared, i.e. matrices with only the diagonal, sub-diagonal and super-diagonal, non-vanishing. These matrices often arise in physical problems. It can be shown that any matrix can be reduced to tri-diagonal form (using a finite sequence of unitarily similar transformations).  Hence it is worthwhile to study these objects in their own right.

(Subscript[α, i_, j_]) := .../;\[LeftBracketingBar]i - j\[RightBracketingBar] ≤ 1

(Subscript[α, i_, j_]) :=  ...

(Subscript[𝒯, n])_ := Table[ ... ]

(Subscript[α, i_, j_]) := (Subscript[β, i, j]) /; \[LeftBracketingBar]i - j\[RightBracketingBar] ≤ 1

(Subscript[α, i_, j_]) := 0

(Subscript[𝒯, n])_ := 
      Table[(Subscript[α, i, j]), {i, n}, {j, n}]

(Subscript[𝒯, 5])

(Subscript[δ, n])_ := (Subscript[δ, n]) = 
            Det[(Subscript[𝒯, n])]

(Subscript[δ, 1])

(Subscript[δ, 2])

(Subscript[δ, 3])

PolynomialReduce[(Subscript[δ, 3]), {(Subscript[δ, 2]), (Subscript[δ, 1])}, 
      Flatten[Table[(Subscript[β, i, j]), {i, 2, 3}, {j, i + 1, 3}]]]

PolynomialReduce[(Subscript[δ, 5]), {(Subscript[δ, 4]), (Subscript[δ, 3])}, 
      Flatten[Table[(Subscript[β, i, j]), {i, 2, 5}, {j, i + 1, 5}]]]

(Subscript[δ, 4])

(Subscript[ℰ, n])_:=(Subscript[x, n]^(′′))(t)==(Subscript[x, n - 1])(t) + (Subscript[x, n + 1])(t) - 
            2 (Subscript[x, n])(t)

(Subscript[ℰ, n]) /. 
      (Subscript[x, n])_ -> 
        Function[
          t, (Subscript[a, n]) (E^(I ω t))]

\(TraditionalForm\`%\/E\^\(I\ t\ ω\)\)

Simplify[%]

% /. 
      (Subscript[a, j])_ -> 
        a (E^(I (2  j γ + δ)))

\(TraditionalForm\`Simplify[%\/E\^\(I\ \((2  j\ γ + δ)\)\)]\)

% /. (ω^2) -> ω2

ω2 /. First[Solve[%, ω2]]

FullSimplify[%] /. γ -> (γ/2)

\(TraditionalForm\`Simplify[\@%%%, γ > 0]\)

TrigExpand[%]