Algebra

from Essential Mathematica for Students of Science
© James J. Kelly, 1998
Last revised: January, 2006

In this notebook we introduce the techniques for manipulating algebraic expressions, solving equations, and finding eigenvalues and eigenvectors.

Initialization

ClearAll["Global`*"];
Off[General::spell,  General::spell1]

Needs["Utilities`Notation`"]

Manipulation of algebraic expressions

Rational expressions

Together collects terms over a common denominator and cancels common factors.

(a/(b + c)) -  (d/(e - f)) // Together

Apart separates an expression into a sum of terms with simpler denominators.

% // Apart

Cancel cancels common factors between the numerator and denominator of a rational expression.

((2 (a^2) b)/((b + c) a))  // Cancel

The Numerator or Denominator of a rational function of polynomials can be extracted using functions of the same name.  Note that the numerator is considered the part of the expression which does not have a "superficially" negative power.

expr  =  ((a (b^(-2)) (g^(-n)))/((c - d) (((e - f))^(-1))));
{Numerator[expr], Denominator[expr]}

Note that the option Trig->True is needed for use with trigonometric expressions to avoid automatic conversions of expressions of the type 1/Cos[x]->Sec[x].

expr = Sin[(x^2)]/Cos[(x^3)];
{Numerator[expr], Denominator[expr], 
      Denominator[expr, Trig -> True]}

Combine: (x/(1 - (x^2))) - (1/(1 + x))

Expanding products

ExpandNumerator expands numerators while ExpandDenominator expands denominators.  (Duh!)

(((a + b) (b - 2  a))/((a - b) a))  // 
      ExpandNumerator

(((a + b) (b - 2  a))/((a - b) a))  // 
      ExpandDenominator

Expand expands numerators leaving denominators in factored form, cancelling common factors where possible.

(((a + b) (b - 2  a))/((a - b) a))  // 
      Expand

ExpandAll expands both numerators and denominators.  Note that common factors are not cancelled.

(((a + b) (b - 2  a))/((a - b) a))  // 
      ExpandAll

Cancel common factors from the expression above.

Factoring sums

Factor reduces an expression to a product of factors.

(a + b) (2  a - b) (c - a) (d + b)  // 
      Expand

% // Factor

(((a + b) (2  a - b))/((c - a) (d + b)))  // 
      ExpandAll

% // Factor

By default, Factor is limited to integers and thus will not factor

(x^2) + 1  // Factor

The option GaussianIntegers->True can be used to extend Factor to complex numbers.

Factor[(x^2) + 1, GaussianIntegers -> True]

Options[Factor]

Similarly, with the option Trig->True multiple-angle formulae are used to factor expressions involving trigonometric and hyperbolic functions.

Factor[Sin[2  x] Cos[2  x], Trig -> True]

Factor[Sinh[2  x] Cosh[2  x], Trig -> True]

Factor : -24 x + 12 (x^2) - 24 y + 4 x y + 
        10 (x^2) y - 3 (x^3) y - 8 (y^2) + 10 x (y^2) - (x^2) (y^2) - 
        (x^3) (y^2) + 2 x (y^3) - (x^2) (y^3)

Factor: 
a) Cos[x]+Cos[y]	b) Cosh[x]+Cosh[y]

Collecting terms

Collect[expr,x] collects terms involving the same powers of x.

(a  +  b x) (c x  - d (x^2))  // 
      Collect[#, x] &

Collect[expr,x,f] applies the function f  to each coefficient separately. Thus, Collect[expr,x,Simplify] can be used to simplify the coefficients or Collect[expr,x,Factor] to factor them.

(a  +  b x) (a c x  - d (x^2))  // 
      Collect[#, x, Factor] &

Collect[expr,{x1,x2}] collects terms involving the same powers from a list of objects.

(a y  +  b x + f) (c (y^2)  - d (x^2)  - g)  // 
      Collect[#, {x, y}] &

(a y  +  b x + f) (c (y^2)  - d (x^2)  - g)  // 
      Collect[#, {y, x}] &

Note that collection priority is ordered according to the listed order of the variables.

FactorTerms[expr,x] extracts any that are independent of x.

-a d x - b d y + a c x y + b c (y^2)  // 
      FactorTerms[#, x] &

-a d x - b d y + a c x y + b c (y^2)  // 
      FactorTerms[#, {x, y}] &

With only one argument an overall rational factor is pulled out of a polynomial.

2 x  +  (1/2)  +  ((x^2)/3)  // FactorTerms

Coefficients are obtained using Coefficient or CoefficientList.  Coefficient[expr,x] reports the coefficient of x found in expression expr; note that for this purpose (x^2) is considered distinct from x, such that the coefficient returned will not contain any power of x.  Coefficient[expr,x,n] returns the coefficient of (x^n).

expr  =  (((a  + 3  x y  + (b/x)))^15);
{Coefficient[expr, (x^7)], Coefficient[expr, (x^2), 7], 
      Coefficient[expr, x, -3]}

CoefficientList[expr,x] returns a list of coefficients beginning with (x^0).  Coefficient[expr,{x1,x2,...}] returns a rectangular array.  Note that since both Coefficient and CoefficientList are both intended for use primarily with polynomials, functions other than powers may appear in the coefficients.  Also, CoefficientList does not handle negative powers.

expr  =  (((x  +  y Cos[y]))^5);
CoefficientList[expr, x]

CoefficientList[(((x + y))^5), {x, y}]  // MatrixForm

Let expr = Nest[μ #(1-#)&,x,4].
a) Use Collect to produce an explicit polynomial in μ and to factor the coefficients.
b) Alternatively, use CoefficientList to list the coefficients of (μ^n) and then factor those coefficients.

Deduce the linear transformation x -> a z + b that takes an arbitrary quadratic function α (x^2) + β x + γ to the simpler form (z^2) + c by determining {a, b, c} in terms of {α, β, γ}.

Using the AlgebraicManipulation palette

There exists an AlgebraicManipulation palette which can be used to perform many of the operations by point-and-click rather than by entering commands directly.  If you are unsure of whether that palette is active, use the Window submenu of the main toolbar to display a list of the active windows.  If AlgebraicManipulation is present, clicking on it will bring it to the fore; otherwise use File ↦ Palettes to activate it.  The palette contains a list of operations that can be applied to a highlighted expression or subexpression.  Highlight the following expression or a portion thereof. Then apply one of the transformations, such as Expand, by clicking on the palette.  That operation is performed in situ, leaving the result highlighted.  Clicking on Factor undoes the previous the operation.  Expand followed by Simplify produces another variation.  Experimenting with various transformation sequences can produce the desired result more quickly than typing and editing commands.  Of course, if you expect to repeat these evaluations later, you should record the best sequence of transformations in the notebook itself.

-(2 a - b) (a + b) (a - c) (b + d)

This method can be applied to parts of expressions also.  To get some familiarity with this technique, perform the following operations on a copy of the expression below.
	1) put the first term in the numerator over a common denominator
	2) factor the second term of the numerator
	3) simplify the numerator
	4) try simplifying the denominator — the result is not likely to be satisfactory
	5) try other sequences on the whole or parts of the denominator

(((a  + b (x/(c - d))) ((c^2) - (d^2)))/(a + (b/(c + (d/(x - d)))) + (1/(c x  - d c + d))));

Manipulation of trigonometric expressions

Many functions are also available for manipulation of expressions involving trigonometric functions, often appearing as counterparts to functions that transform algebraic expressions in similar ways.  Both circular and hyperbolic trigonometry are handled by these functions.  In addition, there are functions which transform between trigonometric and exponential representations.

TrigExpand expands products of trigonometric functions into a sum of terms and applies multiple-angle formulas to reduce the arguments to their most basic form.

Cos[2  x] Sin[2  x]  // TrigExpand

Cosh[2  x] Sinh[2  x]  // TrigExpand

TrigFactor factors expressions into products of terms.

2 Cosh[x]^3 Sinh[x] + 2 Cosh[x] Sinh[x]^3 // 
      TrigFactor

TrigReduce uses trigonometric identities to simplify expressions, generally attempting to reduce the number of trigonometric functions involved, often using multiple-angle formulas.

Sin[x]^2 + Cos[x]^(2 )

Sin[x]^2 + Cos[x]^2 // TrigReduce

2 Cos[x]^3 Sin[x] - 2 Cos[x] Sin[x]^3 // TrigReduce

TrigToExp converts trigonometric functions to exponential form.

stuff = Cosh[2  x]^2  -  Sinh[x]  // TrigToExp

ExpToTrig converts from exponential to trigonometric functions.

stuff // ExpToTrig

Often several transformations can be profitably combined.

(∑_(m = -j)^j) Exp[((m y)/j)] // ExpToTrig // 
      TrigFactor

Identities can also be checked, but it will often be necessary to use Simplify to obtain comparable expressions in both arguments of Equal.

Cos[x + y]^2 == (1/2) (1 + Cos[2 x + 2 y]) // 
      Simplify

Expand:
a) Cos[4x]	b) Tanh[2x]

Simplify using TrigReduce: \(TraditionalForm
      \`2 \(\( cosh(x)\)\^3\) \(sinh(x)\) + 
        2 \( cosh(x)\) \(sinh(x)\)\^3\).

Expand and place over common denominator: 
a) cos(x + y)		b) tanh(x - y)

Expand the following product and use TrigReduce to simplify the result: \(TraditionalForm
      \`\(\((\(cos(x)\) \(cos(y)\) - \(sin(x)\) \(sin(y)\))\)\^2 . \)\)

Use ExpToTrig and TrigFactor to simplify \(TraditionalForm\`\(1 + e\^x\)\/\(1 - e\^x\)\).  Compare with the effect of Simplify.

Watch out for that branch!

Mathematica tends to be more careful about its algebra than the average user.  For example, many new users are puzzled by the "failure" of Simplify to reduce the following expression

(Sqrt[(a^2)]) // Simplify

or Expand to alter

Log[a b]  // Expand

but Mathematica does not automatically assume that the positive branch of the square root is desired or that variables are positive or even that they are real.  If you are confident that all relevant quantities are positive, PowerExpand can be used to simplify powers and products, expand logarithms, and perform related transformations that are valid for positive quantities.  Thus,

(Sqrt[(a^2)]) // PowerExpand

gives the positive branch of the square root and

Log[(a^2) (b^3)] // PowerExpand

separates the log of a product into the sum of their logs assuming that both factors are positive and extracts their exponents properly.  However, although PowerExpand can be very useful, you must ensure that all affected expressions satisfy its assumptions.  Thus, if a > b and c > 0 and all are real, then

Log[(a - b) c] // PowerExpand

is handled correctly, but otherwise you must determine the appropriate branch of the complex logarithm function.

Similarly, ComplexExpand assumes that all variables are real, but not necessarily positive.  Thus,

(Sqrt[(a^2)]) // ComplexExpand

still cannot alter this expression because (Sqrt[(a^2)]), though real, could still be ±a.  Similarly,

Log[a b]  // ComplexExpand

involves the argument function Arg because the sign of a b is not known a priori.  Nevertheless, ComplexExpand is helpful for expressions like

Sin[x + I y]  // ComplexExpand

where explicit real and imaginary parts are expressed in terms of real variables.  ComplexExpand[expr,{x1,x2,⋯}] assumes that {x1,x2,⋯} are complex but that all other variables are real.

Exp[(z^2) (x  +  I y)]  // 
      ComplexExpand[#, {z}] &

More information about the algebra of complex variables can found in complex.nb.

Expand and simplify coth(x + i y) assuming that x and y are real.

Simplicity is in the eye of the beholder

Perhaps the most useful, but often the most frustrating, Mathematica function for manipulation of symbolic expressions is Simplify, which attempts to reduce an expression to a simpler form.  Simplify includes expansion, factorization, and many other algebraic transformations.  With the option Trig->True, which is the default, it applies trigonometric identities also.  Thus, because Simplify includes all of the functions described above, plus others, the first step in simplifying an expression is usually to append a simplification command in postfix notation, namely expr//Simplify, and seeing what you get.

In the contrived example below this simple operation produces immediate gratification

(4 (Cosh[x]^6 Sinh[x]^2 + 2 Cosh[x]^4 Sinh[x]^4 + 
                Cosh[x]^2 Sinh[x]^6))/(8 + 12 x + 6 (x^2) + 
            (x^3))  // Simplify

whereas for this next expression it is useful to convert from exponential to trig before simplifying:

(((y^2) (E^(y  )))/(((1 + (E^y)))^2))  // ExpToTrig // Simplify

Sometimes when Simplify is stumped, FullSimplify can help.  FullSimplify tries a much wider range of transformations and includes rules for many special functions.  Thus, FullSimplify is useful for the next expression where Simplify is not.

2  n  BesselJ[n, x]  - x BesselJ[n + 1, x] // 
      FullSimplify

However, because FullSimplify must test a broader range of possibilities, it can become extremely time consuming and can exhaust the memory of your computer.  If you hear the hard disk churning, your expression is probably too complicated to handle without human guidance.  It is a good idea to save your work before attempting to apply FullSimplify to an unwieldy expression in case it becomes necessary to abort the evaluation, which is not always without risk!

Furthermore, one often has an aesthetic sense of the form that is desired and there is no guarantee that Simplify shares your notion of simplicity or aesthetics.  You can then attempt to guide the simplification process by applying transformation functions, or explicit replacement rules, in the order you believe will lead to the desired form.  Such a calculation is generally developed by accretion, adding transformations sequentially in postfix form until the goal is reached.  Sometimes the same transformation has to be applied several times at different stages of the calculation, perhaps interspersed with simplification commands, with intermediate stages organized using parentheses.  Unfortunately, the appropriate order for these steps can be difficult to work out and  sometimes changes from one release of Mathematica to the next.  There are few general rules we can impart — simplification guided by aesthetics is a fine art acquired only by experience.  On the other hand, unless you have a good reason to insist upon the simplest possible expressions, simplification is often superfluous.  One can calculate and plot functions or investigate many of their properties without getting bogged down in unnecessary manipulations.

Simplify: 
a) 3  sin(x) - sin(3  x)		b) cos((sin^(-1))(x)) 	c) \(TraditionalForm\`\(3 \( tanh(x)\) + \(tanh\^3\)(x)\)\/\(1 + 3 \(\( tanh\^2\)(x)\)\)\)

For the expression 2 (cosh^3)(x) sinh(x) + 2 cosh(x) (sinh^3)(x), compare the following transformations:
a) Factor	b) Factor[expr,Trig->True]	c) TrigFactor	
d) TrigReduce	e) Simplify	f) Simplify[expr,Trig->False]

Express ln(((x + 1)/(x - 1))) in terms of hyperbolic trigonometric functions.  [Hint: you might need the substitution x -> (E^(2  y)) and several steps.]   I leave it to you to judge which form is simplest — most statistical physics textbooks express this formula in terms of hyperbolic trigonometric functions, but rational expressions are pretty simple also.  Often one's preference depends upon context.

Example

The entropy for a system is proportional to the logarithm of the total number of states with the same energy.  Suppose that a system is described by a multiplicity function of the form

multiplicity =  (((((n! )/(m! (n - m)!))))^2);

where n and m are very large positive integers, of order (10^23), with m≪n.  Under those circumstances it is useful to apply the Stirling approximation for Log[n!] in the form

StirlingApprox = {Log[x_!] -> x Log[x] - x};

First, note that Simplify is ineffective because it does not recognize that the variables in this expression are positive and that the logarithm can be reduced without specifying the appropriate branch, but PowerExpand does help.

Log[multiplicity]  // Simplify

Log[multiplicity]  // PowerExpand

Next we apply the Stirling approximation.  Note that parentheses are included to ensure that PowerExpand is executed first.  (Try to do without the parentheses.)

(Log[multiplicity]  // PowerExpand)  /. StirlingApprox

Because we know that m < n, it is convenient to make the substitution m -> x n where 0 ≤ x ≤ 1 is the ratio between m and n.  (The physics would also be discussed in terms of x if this were a course in statistical physics.)  Therefore, we append this replacement rule and another simplification step.

(Log[multiplicity]  // PowerExpand)  /. StirlingApprox  /. 
        m -> x n  // Simplify

Finally, although we probably could use built-in functions, it is easier to simplify these logarithms with explicit rules, using cut-and-paste editing.  It is again necessary to use parentheses to ensure proper order of operations.

((Log[multiplicity]  // PowerExpand)  /. StirlingApprox  /. 
              m -> x n  // Simplify)  /. {Log[n x] -> 
            Log[n]  +  Log[x], 
          Log[n - n x] -> Log[n]  +  Log[1  - x]}  // 
      Simplify

We have now obtained a tidy expression, and could transform it in other ways if desired.  Rather than retaining each step of such a calculation, usually only the last pair of input/output cells would be kept in a notebook.  Although the fairly complicated input cell above would then appear to have sprung from the author's forehead fully armed in Athenian splendor, it would actually have been assembled like an oyster's shell one layer at a time.

FunctionExpand

Expressions involving special functions with complicated arguments can often be simplified using FunctionExpand.

FunctionExpand[Sinh[ArcCosh[x]/2]]

FunctionExpand[Cos[4  ArcTan[x]]] // Simplify

Simplify: Tan[ArcCos[x]/2]

Simplification with assumptions

A powerful new feature of simplification became available in Version 4.0.  It is often possible to simplify symbolic expressions much further when the properties of the various symbols, such as their numerical ranges, are known in advance.  The syntax Simplify[expression,assumptions] permits a list of assumptions to be employed during the simplification of expression.  The assumptions are specified either as a list, {assumption1,assumption2,⋯} or as a logical expression, such as assumption1&&assumption2&&assumption3.  Assumptions can specify the type (domain) for various variables, allowed ranges for values, or relationships between variables.  Assumptions can be employed with Simplify, FullSimplify, FunctionExpand, Refine, Limit, or Integrate.

Consider the following expression.

expr = Log[a - b] + Log[a + b];

Simplify[expr]

Simplification is ineffective without further information about a and b, but if we know that a > b

Simplify[expr, {a > b}]

then the two functions can be combined.  Note that the assumption a > b implicitly contains the additional assumption that both a and b are real numbers, so that their magnitudes can be compared directly.  Similarly, the following expressions can be simplified by employing inequalities to define ranges.

{Simplify[(Sqrt[(((x - y))^2)])], 
      Simplify[(Sqrt[(((x - y))^2)]), x > y]}

Simplify[
      2 (Sqrt[ x y]) ≤ x + y, {x ≥ 0, 
        y ≥ 0}]

Limit[Exp[x y], y -> ∞, 
      Assumptions -> {x < 0}]

Often it is sufficient to specify a domain (Integers, Rationals, Reals, Algebraics, Complexes, Booleans, Primes) for one or more variables.  The function Element[x,domain] specifies that x∈domain is an element of the specified domain.  The element operator ∈ can be entered with the keystroke sequence EscelemEsc.  The variations {x,y,z}∈domain or (x|y|z)∈domain assign each member of a list to the specified domain.  Domains can also be specified for patterns according to patt∈domain.

{Simplify[Cos[n π]], 
      Simplify[Cos[n π], {n ∈ Integers}]}

Simplify[Cos[n π], {((n + 1)/2) ∈ Integers}]

A related function is Refine, which uses assumptions to produce output that more closely approximates the form that you might use when its symbols satisfy explicit numerical assumptions.  The result is often better than produced by Simplify, or at least more specific.

{Simplify[Abs[x], x < 0], Refine[Abs[x], x < 0]}

{Simplify[Log[x], x < 0], Refine[Log[x], x < 0]}

{Simplify[Cos[x + π  y], y ∈ Integers], 
      Refine[Cos[x + π  y], y ∈ Integers]}

If one is performing a long derivation involving many symbols, it can become cumbersome to specify a long list of assumptions many times.  The following example illustrates a useful technique for performing such calculations.  First we create a list a assumptions that apply to the variables of interest.  Then we define our own versions of the simplification functions using these assumptions.  Two examples using these functions follow.

MyAssumptions = {0 < (Subscript[θ, e]) < π, ω > 0, q > Q, Q > 0, 
          0 < ϵ < 1, (Subscript[ε, f]) > 
            0, (Subscript[ε, i]) > (Subscript[ε, f])};
MySimplify = Simplify[#, MyAssumptions] &;
MyFullSimplify = FullSimplify[#, MyAssumptions] &;

sol1 = Solve[ϵ == (((1 + 2 ((q^2)/(Q^2)) Tan[((Subscript[θ, e])/2)]^2))^(-1)), (Subscript[θ, e])][[2]] // MySimplify

Solve[
              (Q^2) == 
                4 (Subscript[ε, i]) ((Subscript[ε, i]) - ω) Sin[((Subscript[θ, e])/2)]^2, (Subscript[ε, i])][[2]] /. ω -> (Sqrt[(q^2) - (Q^2)]) /. sol1 // 
      MyFullSimplify

Although the formatting procedure used by Mathematica still has some peculiarities (e.g., too many minus signs), these results were simplified fairly completely with relatively little effort, whereas much more complicated expressions would have been produced had we not told Mathematica what assumptions it could employ  (Try it!).  In a real problem we would probably use our simplification routines many more times, perhaps adding additional assumptions as the solution is developed.

Under what assumptions will ((((x^m)))^n) reduce to (x^(m n))?  Verify.  List a few examples which show that uncritical use of the proposed replacement rule leads to incorrect results.

Prove (a^n) - 
          (b^n) = (a - 
              b) ((a^(n - 1)) + (a^(n - 2)) 
                b  +  ⋯ +  (b^(n - 1))) for any positive integer n and real a, b.

Symbolic solution of equations

Basic syntax

The basic tools for symbolic solution of algebraic equations are Solve, Reduce, and Eliminate.  Solve[eqs,vars] attempts to solve an equation or list of equations, eqs, for the variable or list of variables, vars.  Equations must be expressed in the form lhs==rhs.

eq1  =  a (x^2) +  b x  +  c  == 0;

sol1 =  Solve[eq1, x]

Solutions are returned as sets of replacement rules.  In this case there are two possible solutions, each containing a single replacement rule for the only independent variable.  A particular solution can then be selected by using the replacement rule of your choice to make an assignment

x1  =  x /. sol1[[1]]

Caution: if you use assignments of the form x=x/.solution you will contaminate the symbol x and may have difficulty revising your equation later.

Sometimes it is useful to verify that solutions are correct by substitution into the original equations.  It is usually necessary to perform some simplification steps before the veracity of the equations becomes manifest.

eq1 /. sol1 // Simplify

A system of equations is expressed in the form of a list of equations, {lhs1==rhs1, lhs2==rhs2, ⋯} or as an expression in which equations are combined using And (&&) in the form lhs1== rhs1 && lhs2== rhs2 ⋯.

Solve[{3  (x^2)  - 2  (y^2)  ==  1,  
        (x^2)  +  4  (y^2) == 3}, {x, y}]

Now Solve returns four solutions, each consisting of two replacement rules.  In other cases, degenerate solutions are duplicated.

Solve[36 + 12 x - 11 (x^2) - 2 (x^3) + (x^4) == 0, x]

Solve does not check for special cases that might arise for some choices of parameters.  For example, the solution proposed for the following equation

Solve[a (x^2) + b x  == 1, x]

does not accommodate the special case a -> 0.  More complete solutions are provided by Reduce, which checks special cases and returns solutions in the form of logical conditions upon both variables and parameters.

Reduce[a (x^2) + b x  == 1, x]

Thus, we recover the original solutions when a ≠ 0, but in addition find a solution for the special case of a = 0 also.  (Actually, the solution for a -> 0 can be obtained by Taylor expansion of the general case but Reduce requires less work.)

Eliminate eliminates a variable or set of variables from a system of equations and returns a smaller system combined with logical operators.  Often this is useful when it is convenient to formulate a system of equations in terms of one or more intermediate quantities that help clarify the relationships between the independent variables.  Sometimes it is useful for transformation of variables.  For example, the following expression transforms a hyperbola and a line from Cartesian to polar coordinates.

Eliminate[{(x^2) - (y^2) == (a^2),  x  + y  == b, 
        x == r Cos[θ], y == r Sin[θ]}, {x, 
        y}]

Solve[eqs,vars,elims] first eliminates variables in elims and then solves the resulting equations for vars, thereby combining the functionality of Solve and Eliminate in a single expression.

Solve the preceding system of equations for {r,θ} and simplify assuming that {a,b,r,θ} are all positive.

Find the roots of LegendreP[6,x] symbolically and verify that these roots are, in fact, real.

A parabola is the locus of all points in a plane that are equidistant from a fixed point, called the focus, and a fixed line, called the directrix.  Construct the parabola whose focus is given by focus={a,b} and whose directrix is the x-axis.  Then determine the coordinates of the vertex, which is the point closest to the directrix.

High-order polynomial equations

Solve produces symbolic or exact numerical solutions to polynomial equations of degree 4 or less, but will also return symbolic roots to higher-order polynomials as well.  Suppose that we wish to determine the fixed points for the (4^th) iteration of the logistic map, which are defined by the following equation.

eq = 
      x == (Nest[μ # (1 - #) &, x, 4] // 
            Expand)

Although there exist no general methods for solving equations of order 16, Solve will nonetheless return Root objects which can be evaluated numerically.  (Open to view)

fixedpoints = Solve[eq, x]

Length[fixedpoints]

The list begins with roots that can be obtained by factoring and is followed by Root objects representing the roots of the remaining irreducible polynomial.  Each Root object has two arguments: a pure function represents the irreducible part of the polynomial and an index identifies a particular root.  These roots can be displayed for a range of μ using Plot.

Plot[Evaluate[x /. fixedpoints], {μ, 2.5, 5}];

The error messages are generated for roots that are complex for some range of μ and are no cause for alarm.  As μ increases more of these roots become real and gradually assemble the pitchfork bifurcation pattern for the logistic map.  More detail about the nature of period doubling and the onset of chaos within deterministic systems can be found in cobweb.nb.

Solve the following equation for the quantity r=y/x and plot the real solutions for {a,-1,1}:
		(2/(y^4)) == (1/(x^4)) + (a/(((x + y))^4))
[Hint: substitute y->r x first and then simplify the equation as much as possible.  (You may need to use Map to perform the same transformations to both sides of the equation.)]

Algebraic equations in transcendental disguise

Solve is designed to handle polynomial equations, but will often attempt to find solutions for more general equations.  For example, in the trivial example

Solve[Cos[x] == Sin[x], x]

two solutions are reported, but the 2π periodicity is overlooked when only the principal branches of inverse trigonometric functions are used.  However, in other cases inverse functions do not any provide solutions at all and numerical solutions must be sought (see below).

Solve[x == Tan[x], x]

Examples

Example: y-scaling for quasielastic scattering

Suppose that a high-energy electron scatters from one of the nucleons in a nucleus with mass (Subscript[m, A]), delivering energy transfer ω and momentum transfer q, such that a nucleon of mass (Subscript[m, N]) is ejected leaving a residual nucleus of mass (Subscript[m, B]).  In the spectator model, the target nucleus is represented at the time of the interaction as the struck nucleon with initial momentum y plus spectator (Subscript[m, B]) with momentum -y, such that the ejectile emerges with final momentum y + q.  Conservation of energy (relativistic) then requires

eq1 = ω + 
            (Subscript[m, A]) == (Sqrt[((Subscript[m, N])^2) + (((y + q))^2)]) + (Sqrt[((Subscript[m, B])^2) + (y^2)]);

Don't worry if you are unfamiliar with the physics behind this equation — we only pose it as an example of a nontrivial algebraic equation which we would rather solve using Mathematica than by hand.

sol1 = Solve[eq1, y] // Simplify

Mathematica produces two solutions, but the ordering of these solutions is to some degree arbitrary and varies from version to version.  We must select the appropriate solution using additional information, such as a physical interpretation or constraint.  For this problem, it is useful to specialize to a target which consists only of a single nucleon without any residual spectator.  Furthermore, we know that the momentum transfer must be greater than or equal to the energy transfer.  Hence, it is useful to solve for ω in terms of q for this special case.

We use Thread to produce new equations for this special case from the list of possible solutions for the general case.

newEqs  = 
      Thread[0 == y /. sol1 /. {(Subscript[m, A]) -> (Subscript[m, N]), 
            (Subscript[m, B]) -> 0}]

Nested expansion and simplification steps are needed to print attractive versions of these equations, but that does not affect their solution.

newEqs  = (newEqs // Simplify // PowerExpand) // 
        Simplify

Solve[newEqs, ω]

It is now clear that the second solution is the one that we want.

Example: nonrelativistic binary collisions

A nonrelativistic binary collision is described by a mass-balance equation of the form m1 + m2==m3 + m4, where we consider m1 to be the projectile, m2 the target, m3 the scattered particle, and m4 the scattered target, but it is also possible that the reactants exchange mass during the collision.  The kinematics of the reaction are governed by conservation of momentum and of energy, which are expressed by the following equations.

collisionEqs = {m1 + m2 == m3 + m4, p1 + p2 == p3 + p4, 
          p1 . (p1/(2 m1)) + p2 . (p2/(2 m2)) + qvalue == 
            p3 . (p3/(2 m3)) + p4 . (p4/(2 m4))};

Here the qvalue measures the inelasticity of the reaction and vanishes for elastic scattering.  Note that we use Dot, with operator notation ., to indicate the scalar product of two vectors.  The momenta are defined using a set of replacement rules.  Note that the since the reaction is confined to a plane, we can omit the z-components and work with vectors in two dimensions.

momentumComponentRules = {p1 -> {p1x, p1y}, 
          p2 -> {p2x, p2y}, p3 -> {p3x, p3y}, 
          p4 -> {p4x, p4y}, p1x -> p1r Cos[theta1], 
          p1y -> p1r Sin[theta1], p2x -> p2r Cos[theta2], 
          p2y -> p2r Sin[theta2], p3x -> p3r Cos[theta3], 
          p3y -> p3r Sin[theta3]};

Because there are 4 independent equations, we can solve this system for 4 variables.  It is convenient to specify three masses {m1,m2,m3}, the projectile momentum {p1x,p1y} in terms of magnitude p1r and angle theta1, the target momentum {p2x,p2y} in terms of magnitude p2r and angle theta2, and the projectile scattering angle theta3.  Mathematica will then happily report a general solution for the remaining mass m4, the momentum p3r for the scattered particle, and the momentum {p4x,p4y} of the recoil particle.  Although we could use Simplify, the minimal return is not worth the considerable investment of time required to simplify expressions of this complexity!

gensol=Solve[
        collisionEqs //. momentumComponentRules, {m4, p3r, p4x, 
          p4y}]

Familiar special cases can be examined by choosing appropriate parameters.  For example, consider the elastic scattering of two particles of equal mass assuming that the target is at rest.

(gensol  /. {m1 -> m, m2 -> m, m3 -> m, 
                qvalue -> 0, theta1 -> 0, p2r -> 0}  // 
            Simplify // PowerExpand)  // Simplify

Clearly the first solution is the one we seek while the second solution is degenerate, describing only the special case of a head-on collision in which the projectile transfers all of its momentum to the struck particle.  More general situations can now be investigated by varying the parameters in the first solution.  A detailed exploration of these possibilities is beyond the scope of this notebook, but can be found in collisions.nb, which also develops functions which produce a variety of useful vector diagrams.

Example: entropy for one-dimensional Ising model

The following equations arise in the one-dimensional Ising model of magnetic systems.  (Don't worry about the physics, we just need a nontrivial system of equations for illustration of Mathematica functions.)

Needs["Utilities`Notation`"]

Symbolize[(Subscript[N, ↑])]; Symbolize[(Subscript[N, ↓])]; Symbolize[(Subscript[N, ↕])];
Symbolize[(Subscript[N, T])]; Symbolize[(Subscript[λ, 1])]; Symbolize[(Subscript[λ, 2])]

eqs  =  {0 == (Subscript[λ, 2]) - B (Subscript[λ, 1]) μ + 
              Log[(Subscript[N, ↑])] - Log[(Subscript[N, ↑]) - ((Subscript[N, ↕])/2)], 
          0 == (Subscript[λ, 2]) + B (Subscript[λ, 1]) μ + 
              Log[(Subscript[N, ↓])] - 
              Log[(Subscript[N, ↓]) - ((Subscript[N, ↕])/2)], 
          0 == 
            2 J (Subscript[λ, 1]) + 
              (1/2) (Log[4] + Log[(Subscript[N, ↓]) - ((Subscript[N, ↕])/2)] + 
                    Log[(Subscript[N, ↑]) - ((Subscript[N, ↕])/2)] - 
                    2 Log[(Subscript[N, ↕])]) , 
          (Subscript[N, T])  ==  
            (Subscript[N, ↑])  +  (Subscript[N, ↓])} ;

Obviously, we can eliminate (Subscript[λ, 2]) and (Subscript[N, ↓]) immediately.

newEqs = 
      Eliminate[eqs, {(Subscript[λ, 2]), (Subscript[N, ↓])}] // Simplify

We find that there are two basic quantities for which substitutions can be used to help simplify the algebra.

newEqs = 
      newEqs  //. {Exp[
                a_.  (B μ + 
                      b_. )] -> (γ^(-a/(2 (Subscript[λ, 1])))) Exp[a b], 
            Exp[a_.  (J  +  
                      b_. ) ] -> (φ^(-a/(4  (Subscript[λ, 1])))) Exp[a b]} // Simplify

Mathematica can now produce relatively simple solutions to these equations. (I preferred the expressions produced by version  4, but these are equivalent.)

sol = 
      Solve[newEqs, {(Subscript[N, ↑]), (Subscript[N, ↕])}] // 
        FullSimplify[#, {γ > 0, φ > 0, 
              (Subscript[N, T]) > 0}] &

At this point we would need to use physical arguments to select the appropriate solution, but that is beyond the scope of this course.  Nevertheless, it should now be more clear how some of these Mathematica functions can be used to simplify and solve systems of equations.

Numerical solution of equations

Polynomial equations

An equation or system of equations involving only polynomials can be solved numerically using NSolve[eqs,vars,n] where eqs is an equation or list of equations, vars is the variable or list of variables, and n is an optional parameter specifying the precision sought in terms of number of decimal digits.  The result is returned as a list of sets of replacement rules.

NSolve[(x^7) - 3  (x^3) == x - 2, x]

NSolve[{(x^2) y - x == 1, (y^2) + x y == -2}, {x, 
        y}]

NSolve will also handle equations which are really polynomials in disguise

NSolve[{(Sqrt[(x^2) y  -  (y^2) x]) == ((1 + x y)^(1/4)), 
        x (Sqrt[ x  -  y]) == 2}, {x, y}]

and many others for which inverse functions are available

NSolve[Cos[x]^2 - Sin[x]^2 == 0.5, x]

NSolve[(E^(2  x)) - (E^x) == 1, 
      x]

The warning message about inverse functions is usually nothing to worry about, and can be disabled using Off, but after solving a problem one should go back and determine whether alternative solutions have been overlooked.

However, NSolve is unable to solve equations which involve transcendental functions in an essentially nonalgebraic manner that makes transformation to polynomial form impossible.

lhs  =  (ξ/(Sqrt[1 -  (ξ^2)]));   rhs  =  
      Cot[ξ];
eq1 =  lhs == rhs;

NSolve[eq1, ξ]

Write a function using NSolve and Cases which returns the real solutions to  the following equation in terms of the parameter a :
		2 == (r^4) + ((a (r^4))/(((1 + r))^4))
Plot these solutions for {a,-1,1}.  [Hint: it is probably easiest to use ListPlot for this purpose.]

Transcendental equations

The more general function FindRoot can be used to find numerical solutions to equations or systems of equations given appropriate starting conditions.  For a single equation, FindRoot[eq,{x, (Subscript[x, 0])}] uses a variant of the secant method to search for a solution in the vicinity of (Subscript[x, 0]) while FindRoot[eq,{x, {(Subscript[x, 0]), (Subscript[x, 1])}}] uses Newton's method based upon two starting values {(Subscript[x, 0]), (Subscript[x, 1])}.  The secant method can be used when symbolic derivatives are possible, but Newton's method must be used otherwise.  Also note that FindRoot[expr,{x, (Subscript[x, 0])}] where expr is an expression rather than an equation is equivalent to FindRoot[expr==0,{x, (Subscript[x, 0])}].

For the simple equation below we recognize by inspection that there is only one root for y > 0 and, hence, need not be too careful in selecting the starting value for the secant method.

eq1  =  y Tanh[y] == 1;

sol1  =  FindRoot[eq1, {y, 0.1}]

To test the accuracy of the solution, we change the equation to a difference between its two sides and then substitute the solution.

eq1 /. Equal -> Subtract /. sol1

Even though our starting value is poor, the accuracy of the solution is quite good; note that the solution contains more significant figures than are routinely printed.

For more complicated equations with several roots, we should compare the left-hand and right-hand sides graphically to choose the appropriate starting conditions.

lhs  =  (ξ/(Sqrt[16 -  (ξ^2)]));   rhs  =  
      Cot[ξ];
eq1 =  lhs == rhs;

Plot[Evaluate[{lhs, rhs}], {ξ, 0, 4}, 
        PlotRange -> {Automatic, {-5, 5}}];

Here we find two positive roots by supplying two sets of starting conditions based upon the figure.

FindRoot[eq1, {ξ, #}] &  /@ {1, 3.5}

Find the first 3 positive roots of tan(x) = 1/x.

Find the first 3 positive roots of BesselJ[2,x].  Then explore the sensitivity to starting values; for example, compare results for starting points of 6.8, 7.0, 7.2.

Find all real solutions to the pair of equations (x^4) + (y^4) = 1 and (e^x) - (e^y) =  1.  [Hint: use ImplicitPlot to display the two equations and to locate appropriate starting values.]

Linear algebra

Matrix manipulation

A matrix is represented as a list of sublists in which each sublist has the same size.  Mathematica provides many built-in functions with obvious names for manipulating and analyzing matrices.  The most important matrix functions are summarized in the accompanying table.

({{function, output}, {Array[mat,{m, n}], m×n matrix with elements mat[i, j]}, {DiagonalMatrix[list], diagonal matrix using elements from list}, {IdentityMatrix[m], m-dimensional unit matrix}, {MatrixForm, print matrix in traditional form}, {Transpose, interchange rows and columns}, {Inverse[A], inverse matrix (A^(-1))}, {Dot[A, B], product A . B for compatible matrices}, {Det, determinant of square matrix}, {Eigenvalues[A] , list of eigenvalues λ,  
              such that A . x = λ x}, {Eigenvectors[A], list of eigenvectors x,  such that A . x = λ x}, {Eigensystem, \(list\ containing\ both\ list\ of\ eigenvalues\ and\ 

              list\ of\ eigenvectors\)}, {LinearSolve[A, rhs], solution vector x to A . x==rhs}, {NullSpace[A], list of vectors x for which A . x==0 }, {MatrixPower[A, b], (A^b)=matrix A raised to power b}, {MatrixExp[A], Exp[A] as defined by power series}})

A few examples of some of these functions are given below.

matrix  = Array[m, {3, 3}]

matrix  // MatrixForm

matrix // Transpose // MatrixForm

The product of two matrices with compatible dimensions is formed by Dot with operator form (.).

var  =  Array[v, {3}];

matrix . var

Evaluate MatrixExp[{{0,θ},{-θ,0}}] and interpret the result.

Matrix form of linear equations

A system of equations of the form A . x == b can be formed using Thread to move the function Equal through the head List.  (What happens if Thread is omitted?)

matrix  = Array[A, {3, 3}];
var = Array[x, {3}];
rhs =  Array[b, {3}];

Thread[matrix . var == rhs]

The result is then a form suitable for Solve.  Alternatively, LogicalExpand

LogicalExpand[matrix . var == rhs]

also produces an expression suitable for input to Solve.  Although these methods are most often used for systems of linear equations, they can be employed for more general problems if the coefficient matrix and variable vector are constructed properly.

Alternatively, LinearSolve[coeff,rhs] is a specialized version for linear equations of the form coeff.x==rhs.  The solution is a vector x which is returned as a simple list, rather than as a set of replacement rules; note that it is not necessary to supply variable names.  LinearSolve works with either symbolic or numerical expressions.  For underdetermined systems LinearSolve returns only one of the possible solutions, whereas Solve returns the general solution.  Also note that for sparse matrices it is usually more efficient to use Solve.

coeff  =  { {1, -2,  (3/8)}, {2, 2, -3}, {(1/2), -1, 
            4}};
rhs = {-1, 3, 0};
LinearSolve[coeff, rhs]

The null space of a matrix is defined by a linear combination of basis vectors satisfying matrix.m==0 and is obtained using NullSpace[matrix].

a  =  {{1, 1, 0}, {-1, 0, 0}, {1, 1, 0}};
b  = NullSpace[a]

a . b[[1]]

Note that LinearSolve and NullSpace can be used for rectangular as well as square matrices, as illustrated by the following exercise.

Let a={{1,2,1},{-1,-1,2}} and b={2,3}.  Determine vectors x and c such that a.(x+λ c)==b for any value of λ.  [Hint: use both LinearSolve and NullSpace.]

Eigensystems

The determinant of a square matrix is given by Det.  Eigenvalues finds the eigenvalues and Eigenvectors the eigenvectors of a matrix, either symbolically or numerically.  Repeated eigenvalues appear with the requisite multiplicity.  Eigensystem returns both in the form of a list with two sublists, the first being the eigenvalues and the second the corresponding eigenvectors.  Eigenvectors are not normalized automatically.  The package LinearAlgebra`Orthogonalization` provides tools to create orthonormal bases.  Rather than attempt to describe all available functions in general terms we prefer to use an example.

Example: linear triatomic molecule

Suppose that three masses move along a line subject to pairwise forces represented by springs.  The figure above displays the interactions between neighboring masses, but omits a possible interaction between the outer masses with spring constant (Subscript[k, 1, 3]).  One can assume, without loss of generality, that the resting lengths of both interior springs are equal to b.  For simplicity we begin by assuming that the two outer masses are equal, (Subscript[m, 3]) = (Subscript[m, 1]), but not necessarily equal to the central mass (Subscript[m, 2]), and that (Subscript[k, 1, 2]) = (Subscript[k, 2, 3]) while (Subscript[k, 1, 3]) = 0.  With these simplifications we can obtain simple but nontrivial symbolic expressions for the frequencies and normal modes of vibration; two different methods are used to illustrate several Mathematica functions.  Numerical results can also be obtained easily for more general systems, but symbolic expressions become unwieldy.

ClearAll["Global`*"];
Needs["Utilities`Notation`"];

Symbolize[(Subscript[m, 1])]; Symbolize[(Subscript[m, 2])];Symbolize[(Subscript[m, 3])];
Symbolize[(Subscript[(Overscript[x, ¨]), 1])]; Symbolize[(Subscript[(Overscript[x, ¨]), 2])];Symbolize[(Subscript[(Overscript[x, ¨]), 3])]

Symbolize[(Subscript[k, 1])]; Symbolize[(Subscript[k, 2])];
Symbolize[(Subscript[k, 1, 2])];Symbolize[(Subscript[k, 2, 3])]; Symbolize[(Subscript[k, 1, 3])]

The equations of motion for the general system, assuming that the spring force is proportional to its length change, are obtained from Newton's second law as follows.

\(\(equations = \ {
		\(m\_1\) x\&¨\_1 == 
            k\_\(1, 2\)\ \((x\_2 - x\_1 - b)\) + \(k\_\(1, 3\)\) \((x\_3 - 
                    x\_1 - 2  b)\), 
		\(m\_2\) 
              x\&¨\_2 == \(k\_\(2, 3\)\) \((x\_3 - x\_2 - 
                    b)\) - \(k\_\(1, 2\)\) \((x\_2 - x\_1 - 
                    b)\), 
		\(m\_3\) 
              x\&¨\_3 == \(k\_\(2, 3\)\) \((x\_2 - x\_3 + 
                    b)\) - \(k\_\(1, 3\)\) \((x\_3 - x\_1 - 
                    2  b)\)\ };\)\)

For our initial special case, these equations reduce to:

eq1 = 
      equations /. {(Subscript[m, 3]) -> (Subscript[m, 1]), (Subscript[k, 1, 2]) -> k, 
          (Subscript[k, 2, 3]) -> k, (Subscript[k, 1, 3]) -> 0}

It is useful to eliminate the resting length with the following change of variables.

Symbolize[(Subscript[(Overscript[η, ¨]), 1])]; Symbolize[(Subscript[(Overscript[η, ¨]), 2])];Symbolize[(Subscript[(Overscript[η, ¨]), 3])]

eq2  = 
      eq1  /. {(Subscript[x, 1]) -> (Subscript[η, 1]) - b, (Subscript[x, 2]) -> (Subscript[η, 2]), 
          (Subscript[x, 3]) -> (Subscript[η, 3]) + b, 
          (Subscript[(Overscript[x, ¨]), 1]) -> (Subscript[(Overscript[η, ¨]), 1]), 
          (Subscript[(Overscript[x, ¨]), 2]) -> (Subscript[(Overscript[η, ¨]), 2]), 
          (Subscript[(Overscript[x, ¨]), 3]) -> (Subscript[(Overscript[η, ¨]), 3])}

Normal modes are defined to be solution vectors of the form η(t) = η(0) 
          Exp[-i ω t].  Thus, if we make the substitution \(TraditionalForm\`η\&¨ -> \(-ω\^2\) η\) for each coordinate, the time dependence divides out.  It is also convenient to replace the equations by differences so that we can extract the coefficient matrix using Coefficient.

eq3 = 
      eq2  /. {(Subscript[(Overscript[η, ¨]), 1]) ->  -(ω^2) (Subscript[η, 1]), (Subscript[(Overscript[η, ¨]), 2]) ->  -(ω^2) (Subscript[η, 2]), (Subscript[(Overscript[η, ¨]), 3]) ->  -(ω^2) (Subscript[η, 3])}

matrix  = Coefficient[#, {(Subscript[η, 1]), (Subscript[η, 2]), (Subscript[η, 3])}] &  /@  (eq3 /. Equal -> Subtract);
matrix // MatrixForm

Therefore, the vibrational frequencies for the normal modes are solutions to the equations matrix.η==0, which is equivalent to the condition Det[matrix]==0.  For an n-dimensional system, this condition yields an (n^th) degree polynomial equation known as the secular or characteristic equation for the system.  To avoid negative frequencies in the present problem, it is easiest to solve for (ω^2) instead of ω itself and to take the positive root later.

freqRule2  = 
      Solve[Det[matrix /. (ω^2) -> ω2] == 
            0, ω2]  // Simplify

freqRule  =  ω -> (Sqrt[ω])2 /. freqRule2

The modes are then found by substituting the frequencies back into the equations of motion and solving for the amplitudes.  Because the overall amplitude is arbitrary for linear systems, it is useful to set one of the amplitudes to unity (but not the central one because of the symmetry of the problem).

Table[
      Solve[eq3 /. 
            freqRule[[
              i]] /. (Subscript[η, 1]) -> 
            1, {(Subscript[η, 2]), (Subscript[η, 3])}], {i, 1, 3}]

The nature of these solutions should now be obvious.  The zero-frequency mode represents the motion of the entire system together without any internal vibration.  In the second mode the central mass remains stationary while the two outer mass vibrate in opposite directions symmetrically.  For the third mode the two outer masses move together while the central mass moves in the opposite direction, with the relative amplitudes proportional to the ratio of masses.

There is a more elegant way to obtain the eigenvalues and eigenvectors.  We can transform the system of equations to the form h . η == -(ω^2) η and then use Eigensystem to obtain both the eigenvalues λ = -(ω^2) and the corresponding eigenvectors directly.  First we extract the right-hand sides of eq2 and divide out the masses.

rhs = #[[2]] &  /@ 
            eq2 /{(Subscript[m, 1]), (Subscript[m, 2]), (Subscript[m, 1])} // Simplify

Then we construct the coefficient matrix.

h = Coefficient[#, {(Subscript[η, 1]), (Subscript[η, 2]), (Subscript[η, 3])}] & /@ 
          rhs;

h // MatrixForm

Finally, we evaluate the eigenvalues and eigenvectors

{λ, vectors} = Eigensystem[h]

and compute the frequencies.

frequencies = Sqrt[-λ]

Note that Eigensystem does not normalize the eigenvectors.  It is often useful to construct a matrix containing the normalized eigenvectors in column form.

modes  =  vectors/(Sqrt[# . #] &  /@ vectors)

Transpose[modes] // MatrixForm

Of course, this procedure can be generalized and streamlined considerably if we formulate the problem in terms of the appropriate matrices from the beginning!

Develop the coupling matrix for arbitrary parameters in the linear triatomic model.  It is more convenient for this and later exercises not to divide out the masses.

Write a function which returns the eigenfrequencies and normal modes given a coupling matrix, h, for which h . η == 
        m . (Overscript[η, ¨]) == -(ω^2) m . η  where η is a coordinate vector and m is a diagonal matrix carrying the masses.

Find the eigenfrequencies and normal modes assuming {(Subscript[k, 1, 2]) = (Subscript[k, 2, 3]) = k, (Subscript[k, 1, 3 = 0]), (Subscript[m, 1]) = (Subscript[m, 3])}.  Verify that the previous results are recovered.  Plot the frequency and ratio of amplitudes for the mode in which the central mass vibrates against the outer masses as a function of the ratio of masses.

Find the eigenfrequencies and normal modes assuming {(Subscript[k, 1, 2]) = (Subscript[k, 2, 3]) = (Subscript[k, 1]), (Subscript[m, 1]) = (Subscript[m, 3])} and including the interaction (Subscript[k, 1, 3]) = (Subscript[k, 2]) between the outer masses.  Verify that the previous results are recovered when the spring constant for the new force is set to zero.  What is the effect of (Subscript[k, 2])?

Generalize to the case of three arbitrary masses.  You will probably find that symbolic expressions for the eigenvalues are too lengthy to display comfortably and are difficult to simplify.  Nevertheless, you can still plot interesting quantities.  For example, plot the two nonzero frequencies together as functions of the mass ratios (either in 3D or in 2D as functions of one ratio for several values of the other).

More examples

Paschen.nb — investigates the energy levels of the hydrogen atom in external magnetic fields, including both weak and strong-field limits.

Table[Factor[(x^n) - 1], {n, 1, 10}] // ColumnForm

hilbert[n_] :=  
      Table[(1/(i + j - 1)), {i, 1, n}, {j, 1, n}]