University College Utrecht  SCI 311: Advanced Mathematics

Initialization

The cells below are initialization cells. These cells are automatically evaluated by Mathematica as soon as the notebook is opened.

(* Hmm..  It doesn't seem like these are automatically evaluated.  You can verify this by leaving this as it is, and producing the very first Plot of the notebook, a linear fit for the first 20 prime numbers.  Then come back here, evaluate these cells, and produce the Plot again.  If the color of the axis and the background of the plot changes, these cells weren't automatically evaluated when the notebook was opened. *) 

(* This turns of the warning that is generated when the name of a new symbol is similar to the name of an already existing symbol. *)
Off[General::"spell"]

(* This defines some specific colors that we will use often. *)
lightbluecolor=RGBColor[220/256,248/256,248/256];
darkbluecolor=RGBColor[34/256,44/256,162/256];
darkgreencolor=RGBColor[0,0.5,0];
redcolor=RGBColor[1,0,0];

(* This makes it the default option for both Plot and ListPlot to draw the axes in dark green on a light blue background and a width of 400 pixels. *)
defaultplotoptions=Sequence[AxesStyle->RGBColor[0,0.5,0], Background->lightbluecolor,ImageSize->400];
SetOptions[ListPlot, Evaluate[defaultplotoptions]];
SetOptions[Plot, Evaluate[defaultplotoptions]];

LU Decomposition

In the Mathematica Help Browser, see 
     Advanced Documentation->Linear Algebra->Linear Algebra in Mathematica, Matrix Decompositions, LU Decomposition
for information on how Mathematica can be used to find this decomposition.

Least Squares

Remark on the notation for vectors in Mathematica

Mathematica does not make a distinction between row and column vectors. This makes it easier to multiply all kinds of vectors and matrices using the . (dot operator). It is very inconvenient, however, if one wants to multiply a column vector with a row vector or a column vector with a matrix. This happens several times in the book, for example in the decomposition of a positive definite matrix on page 25.

We define two vectors a1 and a2.

a1={1,2,3}
a2={a,b,c}

If we multiply a1 and a2 using the dot operator we get the inner product of a1 and a2

a1.a2

This is the reason for writing below:
	P[xList_] := (A.x – b) . (A.x – b) 
and not: 
            P[xList_] := Transpose[A.x – b] . (A.x – b) 
Mathematica does not know how to take the transpose of a vector.
One method to deal with row and column vectors is to define the vectors a1 and a2 as matrices. For example:

a1={{1},{2},{3}}

a2={{a,b,c}}

MatrixForm[a1]

MatrixForm[a2]

MatrixForm[a2.a1]

MatrixForm[a1.a2]

(* It would have been a good idea to Clear[a1, a2] here *)

Remark on MatrixForm

The purpose of the wrapper MatrixForm (see the Mathematica Help Browser for details) is not to make A into a matrix. Its only purpose is to make the output pretty. That is why the assignment  A={{1,0},{1,1},{1,3},{1,4}}  below is inside the argument of MatrixForm. If we would put it outside, as in A =MatrixForm[{{1,0},{1,1},{1,3},{1,4}}] , then the result would be disastrous. Mathematica would not treat A as a matrix, but as a form of output.

Treatment by Mathematica of Example 1 in Section 1.4 in the textbook

As a first illustration we verify the calculations in the example by means of Mathematica. The data in the example are the 4 by 2 matrix A and the vector b in (R^4).

MatrixForm[A={{1,0},{1,1},{1,3},{1,4}}]
MatrixForm[b={0,8,8,20}]
(* Note that A has been defined as a matrix, but the vector b as a list.  A product like c = Transpose[A].b below will thus again be a list *)

The minimizing vector xmin is obtained by solving the following normal equations; we do this using the built-in procedure LinearSolve in Mathematica (see the Help Browser for details).

(* We denote the follow matrix product by K, inspired by the stiffness matrix K = Transpose[A].C.A of Example 2. *)
K=Transpose[A].A;
c=Transpose[A].b;
MatrixForm[xmin=LinearSolve[K,c]]

Next we compute the orthogonal projection p of b onto the column space of A, the error e and its norm, the fact that e is perpendicular to the column space of A and therefore to p in particular, and furthermore that p + e = b.

MatrixForm[p=A.xmin]
MatrixForm[e=b-A.xmin]
min1val=Norm[e]
(* min1val = first value of (first way of calculating) the minimum; a second value will be calculated a bit further on and these two values will be checked to be equal. *)
MatrixForm[Transpose[A].e]
p.e
(* Note that on p. 38 the inner product of p and e is expressed the other way around as (e^T).p *)
p+e==b

The example arises as a minimization problem for a quadratic form, which is given by the inner product of the vector A x - b with itself. Mathematica is unable to take the transpose of a vector, unless it is defined to be a matrix with one row.

P[x_List]:=(A.x-b).(A.x-b)
(* Note that P[x] is defined here to accept only a list x (and not a matrix) as an argument *)
x={x1,x2};
Expand[P[x]]
(* Without Expand the polynomial P[x] would directly be expressed as a sum of squares. *)

We use FindMinimum to minimize this function. The output to FindMinimum consists of two parts, the value of the minimum and a point where it is attained. 
First, we check that the two methods yield the same vector where the minimum is attained.

(* Search for a local minimum of P[x] in a neighborhood of the origin {x1, x2} = {0, 0}. *)
{minimum,location}=FindMinimum[P[x],{{x1,0},{x2,0}}];
(* The output of FindMinimum is of the form {minimum, x -> xxmin}.  The following command returns xxmin *)
xxmin=x/.location
(* Note that here xxmin = {1., 4.} is displayed with a decimal point while xmin = {1, 4} was not.  This is Mathematica's way of saying that xmin was an exact result, whereas xxmin is a approximate numerical result.  This is why we now check whether the two values indeed coincide (see the Help Browser for details on how == evaluates equality of approximate numbers) *)
xmin==xxmin

Furthermore, we verify that both methods yield the same value for the minimum.

(* Remember that P[x] = ((\[LeftDoubleBracketingBar]e\[RightDoubleBracketingBar])^2), so we need to take the square root to compare it to our other calculation min1val of the error. *)
min2val=(Sqrt[minimum]);
(* Here too, mimimum = 6.6332495807108005` is an approximate numerical result, whereas min1val = 2 (Sqrt[11]) was an exact result, which is why we want to compare them. *)
min1val==min2val

We use Clear to remove the values assigned to the variables we just introduced.

Clear[A,b,K,c,xmin,p,e,min1val,P,x,x1,x2,minimum,location,xxmin,min2val]

A final remark: all built-in Mathematica expressions begin with a capital letter, while expressions like E and I have a fixed meaning already. Therefore it is good practice to use lower case to denote your own expressions in Mathematica code. Furthermore, it is advisable to make them longer than the single letter ones used above. For instance, as it happens, A was not yet in use, but look at

?K

?C

?D

The built-in way

Now we consider a new application of the method of least squares. Here is a table of the first 20 primes in the natural numbers.

tableprimes=Table[Prime[n],{n,20}]

The following is a plot of these data (in correct English a plural).

plotprimes=ListPlot[tableprimes,PlotStyle->PointSize[0.015`]]
(* The ` in 0.015` is a so-called NumberMark, denoting the effective number of digits in the precision of 0.015. *)

From the Mathematica Help Browser we learn about the built-in procedure Fit. 
This gives a linear fit to the data.

lfit=Fit[tableprimes,{1,t},t]

Below is a plot of the linear fit.

plotlfit=Plot[lfit,{t,0,20},PlotStyle->Blue]
(* What would happen if we remove this PlotStyle specification? Blue is already Mathematica's default color *)

Next we seek a quadratic fit to the data.

qfit=Fit[tableprimes,{1,t,(t^2)},t]

Here is a plot of the quadratic fit.

plotqfit=Plot[qfit,{t,0,20},PlotStyle->Red]

We show both fits superimposed on the original data. Visually, the quadratic fit is better than the linear one.

(* We remove the tick marks on the axis to remove all distraction for a visual comparison of the two fits *)
Show[plotqfit,plotlfit,plotprimes,Ticks->None]

Comparing the errors

We compare the errors in the linear and quadratic approximation.

lfittable=Table[lfit,{t,20}];
qfittable=Table[qfit,{t,20}];
lerror=Table[ Abs[tableprimes[[n]]-lfittable[[n]]],{n,20}];
qerror=Table[ Abs[tableprimes[[n]]-qfittable[[n]]],{n,20}];

pl1=ListPlot[lerror,PlotStyle->{Blue,PointSize[0.02`]},Joined->True];
pl2=ListPlot[qerror,PlotStyle->{Red,PointSize[0.02`]},Joined->True];
Show[pl1,pl2]

With linear algebra

Now we want to approach this problem by means of linear algebra. Instead of leaving everything to the built-in function Fit, we set up the equations ourselves.  
We wish to determine the quadratic polynomial whose values at the natural numbers approximate as much as possible those of the corresponding primes. The column of values of the polynomial at the points 1 through 20 is a linear combination of the columns of the following 20 by 3 matrix amatrix

MatrixForm[amatrix=Table[{1,n,(n^2)},{n,20}]]

In order to find the coefficients of this polynomial we must solve the normal equations; we do this using LinearSolve. Geometrically, this comes down to projecting the vector tableprimes belonging to (R^20) orthogonally onto the three-dimensional linear subspace in (R^20) spanned by the columns of amatrix.

coeff=LinearSolve[Transpose[amatrix].amatrix,Transpose[amatrix].tableprimes]
myfit=coeff.{1,t,(t^2)}

plotmyfit=Plot[myfit,{t,0,20},PlotStyle->Magenta];
Show[plotmyfit,plotprimes,Ticks->None]

Exercise. Show that myfit evaluated at the numbers from 1 through 20 equals the orthogonal projection onto the appropriate space. Then check with Mathematica  that the error is indeed perpendicular to the approximation.

Now we check that our fit agrees with the quadratic fit obtained with Fit.

qfit-myfit
Chop[qfit-myfit]

Exercise. Find a third-degree fit with the linear algebra method.

On larger intervals there are much better ways to approximate the distribution of primes. The Prime Number Theorem tells to use logarithms.

(Optional) Approximation of the kth prime number and the number of primes

The k–th prime number(Subscript[p, k])(k≥6)canbeboundedinthefollowingway:k(log(k)+log(log(k)))≥(Subscript[p, k])≥k(log(k)+log(log(k))-1)
Here are these bounds for the first 1000primes visualized.

With[{n = 1000},
Print[Show[Graphics[{Thickness[0.002],
 (* upper bound *)
 {Red, Line[Table[{k, k (Log[k] + Log[Log[k]])}, {k, 2, n}]]},
 (* lower bound *)
 {Green, Line[Table[{k, k (Log[k] + Log[Log[k]] - 1)}, {k, 2, n}]]},      
 (* the primes *)             
 {Blue, Line[Table[{k, Prime[k]}, {k, 2, n}]]}}],AspectRatio->.6,
      PlotRange -> All, Axes->True, Background->lightbluecolor, ImageSize->800] // N]]

Next we give a more precise approximation of (Subscript[p, k]). The picture  on the left shows a plot of the approximation in black (thick line) and the primes in red.  The right picture shows the relative error.

pApprox[k_] = k (Log[k] + Log[Log[k]] - 1 + (Log[Log[k]] - 2)/Log[k] -(Log[Log[k]]^2 - 6 Log[Log[k]] + 11)/Log[k]^2);

With[{n = 10^4},
Print[Show[GraphicsArray[
Block[{$DisplayFunction = Identity},
  {Plot[pApprox[k], {k, 2, n}, PlotStyle->{Thickness[0.02]},
        (* show primes as points *)
       Epilog -> Table[{Red, Point[{k, Prime[k]}]},
                       {k, 2, n}], AxesStyle->RGBColor[0,0.5,0], PlotRange -> All],
   ListPlot[Table[N[{k, (pApprox[k] - Prime[k])/Prime[k]}], 
                    {k, 2, 10000}], PlotStyle->RGBColor[1,0,1]
		    ]}]], ImageSize->800]]]

Determining a good approximation of π(n), i.e. the number of primes smaller or equal to n, is one of the most interesting problems in mathematics. Here are some estimates:

(* modeling step curve *)
exactπ[n_] := Line[(Sequence @@ {# - {0, 1}, #})& /@
          Thread[{Table[Prime[i], {i, n}], Table[i, {i, n}]}]];

(* Hadamard’s approximation *)
Hadamardπ[n_] := Line[Table[N[{x, x/Log[x]}],
                            {x, 2, Prime[n], Prime[n]/120}]];
(* Gauss’s approximation *)
Gaussπ[n_] := Line[Table[N[{x, LogIntegral[x] - LogIntegral[2]}],
                            {x, 2, Prime[n], Prime[n]/120}]];
(* Guiasu’s approximation *)
Guiasuπ[n_] := Line[Table[N[{x,
Function[h, 1/h (LogIntegral[x^(1./h)] - LogIntegral[2.]) +
            Plus @@ ((LogIntegral[x^(1/#)] - LogIntegral[
            x^(1/(# + 1))])/#& /@ Reverse[N @ Range[h - 1]])][
                                     Floor[Log[2, x] // N]]}],
          {x, 2, Prime[n], Prime[n]/120}]];

The exact curve exhibits a characteristic step shape; the other three are good approximations for large numbers.

PrimeNumberSet[n_,opts___]:=Show[Graphics[(*the four curves*){exactπ[n],{AbsoluteDashing[{7,7}],Hadamardπ[n]},{AbsoluteDashing[{5,5}],Gaussπ[n]},{AbsoluteDashing[{3,3}],Guiasuπ[n]}}],opts,Axes->True]

Needs["PlotLegends`"]

styledText[text_]:=Text[text,Background->Blue];

ShowLegend[PrimeNumberSet[50,ImageSize->800,"Background"->lightbluecolor],(*the rest is all for the legend*){{{Graphics[{Black,Line[{{0,1},{4,1}}]}],styledText["exact"]},{Graphics[{Black,AbsoluteDashing[{5,5}],Line[{{0,0},{4,0}}]}],styledText["Hadamard"]},{Graphics[{Black,AbsoluteDashing[{3,3}],Line[{{0,-1},{4,-1}}]}],styledText["Gauss"]},{Graphics[{Black,AbsoluteDashing[{2,2}],Line[{{0,-2},{8,-2}}]}],styledText["Guiasu"]}},LegendPosition->{0.2,-0.4},LegendSize->{0.65,0.4}}]

Clear[tableprimes,plotprimes,lfit,plotlfit,qfit,plotqfit,lfittable,qfittable,lerror,qerror,pl1,pl2,amatrix,coeff,myfit,plotmyfit,myfittable,error,amatrix3,
myfit3,plotmyfit3,pApprox,exactπ,Hadamardπ,Gaussπ,Guiasuπ,PrimeNumberSet,styledText]