Numerical and Symbolic
Solution of Equations

by
K. D.  Stroyan
University of Iowa

copyright 1997 by Academic Press, Inc. - All rights reserved.

Special Mathematica Functions

Needs["Graphics`Colors`"];

Notebook Overview

This NoteBook shows various ways to solve equations.  Exact symbolic solutions are not always possible and then we need to use approximate methods.  The approximate methods are sensitive to scales.

References to the Text

This NoteBook goes with Chapter 11 of the text in particular because we need to solve various equations in MAX-min problems.

Sample Solution of Equations

The Quadratic Formula

Mathematica has good algorithms to find solutions of algebraic equations.  The Simplest case is the quadratic formula that you learned in high school.  <Enter> the following computations:

Clear[x, y, a, b, c]y := a (x^2) + b x + cy

Solve[y == 0, x]

Now solve the equation with numbers :  y = x2 + 3 x + 2.

a = 1b = 3c = 2Solve[y == 0, x]

y

Next use  y = 2 x2 + 2 x - 3.  Don't forget  c  is negative in this case.

Cell[BoxData[{
    \(a = \(?;\)\), 
    \(b = 
      \(?;\ \ (*Change\ \ y\ \ and\ solve\ the\ new\ equation.*)",LongForm\
->False]\)\), 
    \(c = Solve[y == 0, x]\)}], "Input",
  AspectRatioFixed->True,
  FontFamily->"Times"]

Last try  y = 2 x2 + 2 x + 3.

\(a\  = \ 
      \(?
b\ =\ ?	(*Change\ \ y\ \ again\ and\ solve\ the\ equation.*)
        \)\)\(c\  = \ \(?
Solve[y\ ==\ 0\ ,\ x]\)\)

Complicated Algebraic Solutions "ToRules"

Mathematica's solve can lead to solutions that are not expressed in radicals.  Here is what they look like

f := (x^5) - 3 (x^4) + 3 x + 2Solve[f == 0]

These solutions "ToRules" can still be useful and seen as numerical approximations as follows:

N[Solve[ f == 0 ]]

NSolve[ f == 0 , x]

"Essentially" Non-Algebraic Equations:

Mathematica's NSolve[.] can not deal with equations that involve "too many" transcendental funcrtions like Sin[.], Cos[.], Exp[.] and  Log[.].  We saw one such example in Exercise 2.?? where we want to solve the equation

s - Log[s]/c = k

for constants  c  and  k.  (This is the invariant equation for an S-I-R epidemic.)

Mathematica's NSolve[.] can not find this solution because of the Log[.]:

c = 6.8; k = (1/3.0) + (1/30.0) - Log[(1/3.0)]/cNSolve[s - Log[s]/c == k, s]

Mathematica's  FindRoot[ expr1 == expr2, {var, guess}]  will give us an accurate approximate based on a guess at the solution:

c = 6.8; k = (1/3.0) + (1/30.0) - Log[(1/3.0)]/cFindRoot[s - Log[s]/c == k, {s,  .2}]

Numerical Scaling and FindRoot[.]

Numerical solutions cannot tolerate vast differences in scale.  Planck's formula (from the Scientific Projects) for intensity in terms of wavelength is:

Clear[a, b, T, l, dl, i]; i[l_] := (a/((l^5) Exp[(b/(l T))] - 1)); i[l] dl

which has the specific physical parameters (in mks units)

\(h = 6.6255*10\^\(-27\); 
c = 2.9979\ *10\^10; 

    k = 1.3805*10\^\(-16\); 
a = 8\ N[π]\ h\ c\^2; 
\(b = \(h\ c\)\/k; 
    \)\)i[l]

Suppose we fix the temperature T = 1000 and try to Find the Root of  di/dl = 0 without any attention to the sizes of the constants.  This leads to an incorrect solution of  l = 0.01.

T = 1000; FindRoot[i'[l]==0,{l,  .01}]

While  l = 0.007  is NOT correct, the value of the derivative actually is small (relative to 1) at this value of  l.  Let's see why:

i'[0.00742649]

To obtain a useful accurate solution, we make the change of variables,  x = b/(l T), so that the critical point equation is 5 = (5 - x) Exp[ x ], with all parameters on a scale around  5.  We find this root and convert back to  l  as follows:

FindRoot[5 == (5 - x) Exp[x], {x, 5.0}]

x = 4.96511; l = (b/(x T))i'[l]

Notice that the value of the derivative is six orders of magnitude smaller at this value of  l.

The Moral: 
Numerical root finding requires that we scale the problem to about unit size.

Solutions of Systems of Equations

The parabola  y = (x^2)  crosses the circle  (x^2) + (y^2) = 1 in two places.  Mathematica can find these by applying the Solve[.] function to a list of equations as follows:

equns = {y == (x^2), (x^2) + (y^2) == 1}Solve[equns]N[Solve[equns]]NSolve[equns]

Notice that there are only  2  real solutions.

Solutions of Systems of Equations with Parameters

The parabola  y = (x^2)  crosses the general ellipse  ((((x/a)))^2) + ((((y/b)))^2) = 1 in two places.  Mathematica can find these in terms of the parameters by applying the Solve[.] function  as follows:

equns = {y == (x^2), ((((x/a)))^2) + ((((y/b)))^2) == 1}Solve[equns, {x, y}]

MAX-min Examples

The following examples are explained in detail in the text.

Example 11.1 Max & min of f[x] = (x^3) - 6 (x^2) + 9 x + 1,  for  0  ≤ x ≤ 5

First we find all critical points:

\(Clear[f, x]; 
f[x_]\  := \ x\^3 - 6\ x\^2 + 9\ x + 1; 

    Solve[\(f'\)[x]\  == \ 0]\)

Since we want a max and min of a smooth function over a compact interval  [0,5], the max and min must occur at one of the endpoints or critical points:

f[0]f[1]f[3]f[5]

\(curve\  = \ 
      Plot[{f[x], Sign[\(f'\)[x]]}, {x, 0, 5}, PlotRange -> {\(-1\), 22}, 
        PlotStyle -> {{Thickness[0.02], Red}, Blue}, 
        DisplayFunction -> Identity]; 

    Show[curve, 
      Graphics[{PointSize[0.03], Brown, Point[{1, f[1]}], Point[{3, f[3]}], 
          Point[{0, f[0]}], Point[{5, f[5]}]}], 
      DisplayFunction -> $DisplayFunction]; 

    Print["Slope ± in Blue"]; \)

Example 11.2 Max of  V[x]  = 4 (x^3) - 48 (x^2) + 144 x for 0 ≤ x ≤ 6

First find the critical points:

Clear[V, x]; V[x_] := 4 (x^3) - 48 (x^2) + 144 xSolve[V'[x]==0]

Next check the candidate points:

V[0]V[2]V[6]

\(curve\  = \ 
      Plot[{V[x], 10  Sign[\(V'\)[x]]}, {x, 0, 6}, 
        PlotRange -> {\(-10\), 130}, 
        PlotStyle -> {{Thickness[0.02], Red}, Blue}, 
        DisplayFunction -> Identity]; 

    Show[curve, 
      Graphics[{PointSize[0.03], Brown, Point[{0, V[0]}], Point[{2, V[2]}], 
          Point[{6, V[6]}]}], DisplayFunction -> $DisplayFunction]; 

    Print["Slope ± in Blue"]; \)

Example 11.3: MAX & min of  f[x] = x + (625/x)  for x > 0, if they exist

Find critical points:

\(Clear[f, x]; 
f[x_] := x + 625\/x; 
Solve[\(f'\)[x]\  == 0]\)

Compute signs of the derivative:

f'[10]f'[30]

This shows that the graph goes down over up.

Find Limits:

Limit[f[x], x -> 0]Limit[f[x], x -> ∞]

Finally, take a look:

\(curve\  = \ 
      Plot[{f[x], 10  Sign[\(f'\)[x]]}, {x, 0.01, 50}, 
        PlotRange -> {\(-10\), 100}, 
        PlotStyle -> {{Thickness[0.02], Red}, Blue}, 
        DisplayFunction -> Identity]; 

    Show[curve, Graphics[{PointSize[0.03], Brown, Point[{25, f[25]}]}], 
      DisplayFunction -> $DisplayFunction]; 

    Print["Slope ± in Blue"]; \)

Example 11.4 MAX & min  T[p]  =  (p - 0.50) ((1000/(1 + (p^2))))for 0.5 ≤ p < ∞

Clear[T, p];T[p_] := (p - (1/2)) (1000/(1 + (p^2)));
Solve[T'[p]==0]N[%]

Check signs of the slope on each side of the critical value:

T'[1]T'[2]

The graph goes up over down.

Limiting values:

T[0.50]Limit[T[p],  p -> ∞]

The zero limiting values means  T[(1/2) (1 + (Sqrt[5]))] is maximal and zero at p = 0.50 is minimal.

N[T[(1/2) (1 + (Sqrt[5]))]]

\(curve\  = \ 
      Plot[{T[p], 10  Sign[\(T'\)[p]]}, {p, 0.5, 10}, 
        PlotRange -> {\(-10\), 315}, 
        PlotStyle -> {{Thickness[0.02], Red}, Blue}, 
        DisplayFunction -> Identity]; 

    Show[curve, 
      Graphics[{PointSize[0.03], Brown, 
          Point[{1\/2\ \((1 + \@5)\), T[1\/2\ \((1 + \@5)\)]}]}], 
      DisplayFunction -> $DisplayFunction]; 

    Print["Slope ± in Blue"]; \)

Exercises on Solution of Equations

Exercise 2.6.3

(a)  In the rubella epidemic with c=6.8, s[0]=1/3, i[0]=1/30, find the invariant  k = s + i - Log[s]/c.

(b)  We know from the text Section 2.5 that i[t] tends to zero as t tends to infinity, so the limiting susceptible fraction also satisfies

\(\(		0 + s[∞] - \(1\/c\) Log[s[0]] = k\)\)

with the same value of k as you computed in part (a) of the exercise.  We could view this as the following question:  Find all solutions in s of the equation

s - (1/c) Log[s] = k

Exercises on Solution of Equations with MAX-min

Exercise 11.3.2

Find the maximum and minimum of each of the following functions over the specified interval.
	
	(a)  f(x) = x3 - 3x2 + 3x + 2,	x ∈ [-2,2]
	
	(b)  f(x) = x3 - 6x + 4,	x ∈ [-3,3]
	
	(c)  f(x) = x3 + 7x + 4,	x ∈ [-10,10]
	
	(d)  f(x) = x3 + 7x + 4,	x ∈ [0,10]
	
	(e)  f(x) = x / (x2 + 1),	x ∈ [-5,5]

Exercise CD11.6.3

Find the point closest to (x,y) = (1,1) on each branch of the hyperbola

((((x/2)))^2) - ((((y/3)))^2) = 1