UC SCI 111, vM+
Notebook 3

Needs["Graphics`Colors`"];
Off[General::spell1];

Mathematica: finding zeroes.

Mathematica contains a whole list of preprogrammed procedures that can be used for finding zeroes of functions. We present a few of them in an order that corresponds with the difficulty of the problems that are solved.

Factorizing polynomials.

A class of functions for which it is more or less easy to find zeroes is the class of polynomials in one variable. When the coefficients of a polynomial p(x) are rational, that is, every coefficient is an integer or a fraction, the procedure Solve can find the roots of p(x). When the degree, n, of p(x) is two, Mathematica will just apply the famous abc-formula.

Solve[a x^2 + b x + c  == 0, x]

Similar formulas exist for degree three and four polynomials. (Try it!) It has been proved that no such formula exists for n bigger than four. Still there exist efficient methods to find the exact roots of polynomials of higher degree.

Consider for example the function f(x) = (x^5)-3 (x^3)- (x^2)+ 3.

f[x_] := (x^5)-3 (x^3)- (x^2)+ 3
Plot[
        f[x], {x, -3, 3}];

The picture gives the feeling that there should be three zeroes. Let's try.

Solve[f[x] == 0, x]

Mathematica finds five! But two of them are not real but so called complex numbers. Here you can recognize them because they contain some root of a negative number, whatever that means. What it means will be explained in one of the next chapters. But anyhow, Mathematica does find the three solutions that we expected.

Exercise 4.
Find the exact position and values of the extrema of the function

	f(x) = (x^6)-3 (x^4)+6 (x^3)-36 x+2

Approximating zeroes.

When we want an answer fast or if we do not need an exact answer, we might be satisfied with an approximation of the answer. That is, we can use a numerical method. The procedure NSolve will try to find all solutions of the equation we feed it.

Clear["@"];
f[x_] := (x^5)-3 (x^3)- (x^2)+ 3;
NSolve[f[x] == 0, x]

Now consider the function

	g(x) = (x^6)-5 (x^3)+ sin x
	
Again we want to know the extrema. Thus we try to solve g'(x)=0.

\(g[x_]\  := \ x^6 - 5\ x^3\  + \ Sin[x]; 
NSolve[\(g'\)[x] == 0, x]\)

Hmm. NSolve doesn't seem to be very happy about our function. Let's have a look at the graph of g.

Plot[g[x], {x, -1, 2}];

The graph clearly suggests the existence of local extrema. Now we could use Newton iteration. Fortunately, Newton's method is implemented in the function FindRoot. But in order tot apply it successfully, we need a first approximation of the zeroes. So let's plot the graph of g'.

Plot[g'[x], {x, -1, 2}];

The graph suggest as starting points x=-.2, x =.2 and x =1.3, for example, and this approach works!

FindRoot[g'[x] == 0, {x, -.2}]FindRoot[g'[x] == 0, {x, .2}]FindRoot[g'[x] == 0, {x, 1.3}]

Exercise 5.
Let f(x) = (sin^3)x +sin x +1. Determine the zeroes of f(x) on the interval [0,10].

Maxima and Minima.

One variable functions.

We start with the simplest optimization problem, the determination of (local) maxima and minima of a one variable function. Let us consider the following example.

f[x_] =  (x/((x^2) + 1))

We are interested in the maximum and minimum. To get an idea of what they are let us plot its graph.

Plot[f[x], {x, -20, 20}];

Can you see where the extrema are? To get a better look, play around with the plot range. It helps you to get a global picture of the graph. Once you have a rough idea of where the extrema are located we determine them by solving f'(x)=0.

NSolve[D[f[x], x] == 0]

Exercise 6.
Repeat the above steps for the functions 

	f(x)=(x^4)-(x^3)-4(x^2)+1 	and		 f(x)=(((x^2)+1)/((x - 1)(x - 2)))
	
Note that some of the functions may not have an absolute maximum or minimum, but local extrema instead. For each of the two functions give a sketch (on paper) of the graph and give the positions of the (local) extrema to two decimal places.

The functions we have considered are so-called  rational functions. That is, they are polynomials or quotients of polynomials.  They have at most finitely many zeroes and that is why Mathematica can usually detect them. It becomes another matter if we consider the so-called  transcendental functions such as the Sine and Cosine.

Repeat the above steps for the function x sin x.  We will find that NSolve does not have the answer for us.  From the graph we note that there are infinitely many local extrema, which explains why Mathematica has some problems with it. In that case we must help Mathematica a bit. In the case of the function x sin x  you have probably seen that the first positive local maximum is near x=2. To determine the value up to 6-digit precision we use the function FindRoot.  First we determine the derivative,

df = D[x Sin[x], x]

And graph this derivative,

Plot[df, {x, -10, 10}];

There is clearly one zero of the derivative in the vicinity of 2. In such a case FindRoot can locate this zero,

FindRoot[df == 0, {x, 2}]

As an aside we mention that Mathematica can compute the zero up to any desired precision. To that end we introduce some options in the command, which you do not need to remember.

FindRoot[df == 0, {x, 2}, WorkingPrecision -> 40, 
      AccuracyGoal -> 35]

Exercise 7.
Consider the function ((sin x)/x). Plot its graph so that you get a good impression of its shape. Then determine the three local extrema at smallest positive values of x to two decimal places.

The law of reflection.

Clear[a]; Clear[b];

First we deal with the law of reflection of lightrays. Draw a sketch (on paper) of the following situation. Consider the point A=(0,1) and a moving point B=(a,b) with a,b positive. Consider the X-axis as a reflecting line. From the point A a light ray leaves, reflects at the X-axis and reaches B. Fermat's principle tells us that the total length of this ray should be minimal. Let us determine at which point of the X-axis the ray should reflect. Call this point (x,0). The total length of the ray now equals

path = (Sqrt[1 + (x^2)]) + (Sqrt[(((a - x))^2) + (b^2)])

To minimize this length we take its derivative and set it to zero. Here is the derivative

D[path, x]

And here we determine its zeros,

Solve[D[path, x] == 0, x]

Wow, we get two solutions! Did you expect that? Which one do you think we should have? After you have made your decision observe that the tangent of the angle of the incoming ray is x and the angle of the reflected ray is (a- x)/b.  Angles are measured with respect to the normal of the reflecting axis.

Exercise 8.
Compute the ratio of these tangents for our optimal ray and verify that it is one. What is your conclusion?

In order to become acqainted with the idiosyncrasies of computer algebra systems we digress briefly on the fact that we seem to have found two local extrema. We consider the special case when a=2 and b=2. In this case we should have the solutions x=-2, 2/3.

specialpath = path  /. {a -> 2, b -> 2}

Plot the function,

Plot[specialpath, {x, -4, 4}];

We see only one local extremum, when we expected two... or did we expect one? Let us solve our derivative=0 problem in this case,

NSolve[D[specialpath, x] == 0, x]

As you see, we have only one extremum, which is as it should be. But then why did we get two local extrema in the general case?

Snell's law.

Here we shall verify Snell's law (in special cases)  using Mathematica . Sketch (on a piece of paper) the following situation. We have a point A=(0,- 1) below the X-axis  and the point B=(2,2) above the X-axis.  We consider a light ray travelling from A to B. In the media above and below the X-axis the speed of light has different values, say speed v- below the X-axis  and speed v+ above the X-axis.  According to Fermat's principle the light ray will travel along the trajectory which takes minimal time. So below the X-axis  it will travel from A to a point (x,0) in a straight line, and then from (x,0) to B in another straight line.  The ratio v+/v- is known as the index of refraction between the two media. Let us denote it by n. The time it takes to travel from A to (x,0) and then from (x,0) to B is proportional to

Clear[path]; 
    path[n_] := 
      n (Sqrt[1 + (x^2)]) + (Sqrt[(((x - 2))^2) + (((2))^2)])

Let us now find the value of x for which this is minimal,

Solve[D[path[n], x] == 0, x]

Perhaps you are now impressed by the computing powers of Mathematica, but unfortunately this answer is not going to do us any good. You may also understand now why we did not choose a general point like (a,b) for B. Instead of finding general expressions it is better to make a function that computes the optimal value of x for each n individually. Here it is,

optimalx[n_] :=  NSolve[D[path[n], x] == 0, x]

Exercise 9.
As we let n vary from zero to infinity, between which values do you think will optimalx[n] vary? You can check this by computing a few values of optimalx[n]. You could even try to make a plot of optimalx[n].

Note that the sine of the incoming ray at (x,0) is equal to

insine = x/Sqrt[x^2 + 1]

The sine of the outgoing ray is equal to

outsine = (2 - x)/Sqrt[4 + (2 - x)^2]

Now observe the value of

outsine/insine  /.  optimalx[1.4]

Put in some other positive values of n as well. The result clearly illustrates Snell's law. After all we have done it must seem amazing that someone in the 16-th  century, like Snell, could solve this problem completely using paper and goose feather. It shows the power of clever thinking as opposed to dogged computation.

Maxima and minima on segments.

The previous three subsections shared the common feature that the range of variables where to look for extrema was unbounded. The search range covered the whole line or the whole plane. In applications however, this is usually not the case. Often one requires the variables to positive. Or one is interested in a function of several variables restricted to a very small range of such variables.

In dimension one constraints are easy. We then deal with functions which are defined on a segment with one or two boundaries. In order to determine extrema we determine the function values in the stationary points and in the extremal points of the segment. The largest or smallest value then gives us the desired extremal function value.

As an example consider the function (x^4)-2 (x^3)+(x^2)+.1 x+1 on the interval [0,1].

f[x_] := x^4 - 2 x^3 + x^2 + 0.1 x + 1;

Here is its graph.

Plot[f[x], {x, -0.2, 1}];

To determine its extrema we determine the stationary points

solutions = NSolve[D[f[x], x] == 0, x]

And the function values in the stationary points and the endpoints of the interval.

{f[0], f[0.604574],  f[0.939443],  f[1]}

The extremal values are now clear.

Exercise 10.
Construct a function f on the interval [-1,1] that has
	• 5 internal maxima
	• 5 internal minima
	• A maximum at x=-1
	• A minimum at x=1