ClearAll["Global`*"]

Roots

You may recall that we used the FindRoot function to determine the algebraic solutions to equations that are otherwise too hard for us. An equation that can’t be solved explicitly is called “transcendental”. Mathematica has many options for solving algebraic equations, whether they are transcendental or not. A “root” is a solution to a transcendental (or other algebraic) equation. Consider this formula:
	10 (sin^3)(θ) - (cos^3)(θ) +4 sin(θ)*cos(θ) + ((θ^2)/3)= 0.
Nobody can solve this for θ algebraically; we can only solve it “numerically”.
It’s a good idea to plot it before we solve it, because there are always subtle things that can go wrong that are visually easy to see.

f1[θ_]=10 ((Sin[θ])^3)-((Cos[θ])^3)+4 Sin[θ]*Cos[θ]+((θ^2)/3);
Plot[f1[θ],{θ,-10,10}]

Notice that the placement of the powers of sine and cosine don’t get placed in the same location that they would be if you were writing a properly formatted math textbook. In Mathematica, there’s no such thing as (sin^2)(x) or (Sin^2)[x], but there are both Sin[(x^2)], and ((Sin[x])^2). This first squares the angle, and the second squares the function (the "correct" version here!). Both of these “Mathematica” notations are ABSOLUTELY WRONG for any written work!

In the plot range shown, there are already four roots of this equation: near θ = -2.25 rad, near θ = +0.25 rad, near θ = +4.25 rad, and near θ = +5.1 rad. We can use FindRoot to solve for all of them by using different initial guesses. FindRoot always give answers that are numbers:

a=FindRoot[f1[θ],{θ,-2.25}]
b=FindRoot[f1[θ],{θ,0.25}]
c=FindRoot[f1[θ],{θ,4.25}]
d=FindRoot[f1[θ],{θ,5.1}]
3*b

Something odd to notice again is that all 4 of these answers appear inside curly brackets, which means that in some sense, Mathematica thinks of these answers as lists or arrays. The attempted math operation 3*b resulted in a non-numeric result, too! That’s problematic if you want to assign a variable to one of these answers, or to do more math using a numerical result. These answers are “lists” that contain several elements: one is a variable name, and another is a number. There’s also the arrow in there.

As we saw last time, this particular object is called a “rule” in Mathematica. When you have a “rule”, it is almost always convenient to “apply” the rule, which, in a sense, means to make it more permanent. 
The operator “/.” does this. Even though it looks a little like division, it’s really a “convert” operator.

bnew=θ/.b
3*bnew

You can see that these results are actual numbers. This operator is called “Replace All”. In the syntax, the first part is the name of the variable that was solved for in the rule (θ, in this case). Then, after the “/.” operator, you put the name of the expression that defines the rule for θ (in this case, we chose b, which was the root nearest to 0.25). You can also combine real math, the FindRoot, and the ReplaceAll operator into a single line:

3*θ/.FindRoot[f1[θ],{θ,0.25}]

Polynomials and lists of answers
While FindRoot is the most generic and powerful root finding tool in Mathematica, there are also a large number of other root-finding tools that are specifically designed to work with “easy” functions. Two of these are “Solve” and NSolve”. The first gives symbolic results, and the second gives numeric results. Let’s try NSolve first:

f2[x_]=(x^3)+2 (x^2) + 6;
NSolve[f2[x]==0,x] (* if "==" is omitted, NSolve assumes you mean "== 0" *)

You’ll notice that there are three answers, but two of them are complex. And, the answer contains multiple curly brackets within curly brackets. If you know you only want the real answer, you can use the “Reals” keyword:

NSolve[f2[x]==0,x,Reals]

As before, you’ll probably want to convert this into an actual number using the ReplaceAll operator.
Consider this quadratic equation with two real roots:

f2[x_]=(x^2)+6 x + 2;
NSolve[f2[x]==0,x]

How do we extract the numeric values for each of these?
The answer given is an array containing arrays of rules... note the double-layers of curly brackets! One option is to use the operators  “First”, “Last”, or “Part” to extract one part.
For example:

x/.First[NSolve[f2[x],x]]
x/.NSolve[f2[x],x][[2]] (* an alternative notation that gives the 2nd result only *)

In the second case, I am using the “Part” command, although you don’t see the word “Part” anywhere. To use it, you append a double-square bracket “[[2]]” to the end of any list to specify which elements you want. I asked for element #2. 
There is a similar function called “Take” which lets you choose specific parts of a list or array, but the result is still always a list, even if it has only one element in it. The second argument of “Take” is a list containing the start and end indices of the items you want:

mylist={1,3,5,9,22};
Take[mylist,{3,5}]  (* start with part 3, end with part 5 *)

We still haven’t tried the “Solve” command yet...

Solve[f2[x]==0,x]

We can also “apply” these two rules as before:

x1=x/.Solve[f2[x]==0,x][[1]]
x2=x/.Solve[f2[x]==0,x][[2]]

“Solve” works even if the coefficients are themselves variables... so "Solve" is an algebra tool! This might be the very best feature of Mathematica.

ClearAll["Global`*"]
Solve[c x + d == 2,x]

An Important Integral: The “Error Function”

Here’s an integral that you can’t do, and neither can your calculator:

Integrate[(E^(-(x^2))),x]

It is the most important integral in statistics, so it obviously relates closely to measurement sciences like physics and engineering.
Symbolically, this answer “Erf” is a function, just like “Sine” or “ArcSecant”. Like all functions, it takes an input number and generates an output number. You don’t know how ArcSecant works, and you don’t care. Why should Erf be different? For example:

Print["Erf of 0.2 is ",Erf[0.2]]
Print["Erf of 0.8 is ",Erf[0.8]]
Print["Erf of 2 is ",Erf[2]]
Print["Erf of 2 as a number is ",N[Erf[2]]]

You can also perform this integral as a definite integral, which is what you’ll usually need:

N[Integrate[(E^(-(x^2))),{x,0,3}]]
NIntegrate[(E^(-(x^2))),{x,0,3}]  (* since it's used so often, they pre-combined these 2 functions! *)

As you’ll see on your homework, this integral is related to probability calculations. Specifically, if you make a series of measurements of some quantity and the average is μ and the standard deviation is σ, then the probability that any subsequent individual measurement will lie between (Subscript[x, min]) and (Subscript[x, max]) is p = (1/(σ (Sqrt[2 π])))(∫_(Subscript[x, min])^(Subscript[x, max]))(e^(-(((x-μ))^2)/(2(σ^2))))dx. This assumes, of course, that the distribution of the kind of measurement errors you make is “normal” or “Gaussian”. In other words, a Gaussian distribution says that you’re more likely to make small mistakes than huge mistakes, and you’re equally likely to make a measurement that is too small as you are to make one that is too large. In fact, this integral defines what a “normal” distribution is. 

Most computational software packages (including Excel) know that this problem is done so often that they give you a direct function to compute it.  Mathematica instead requires you to put together all the steps on your own. Below are shown three ways to do the same problem. 
For this problem, suppose we have made some measurements of the length of a table, and discovered that the length is (3.00 + 0.22) meters (note that just having these numbers means that somebody has already done a lot of work!). If you wanted to know the probability that the next measurement will lie between 2.90 and 2.95, the answer is 8.54%:

μ=3.00;σ=0.22;xmin=2.90;xmax=2.95;
N[Integrate[PDF[NormalDistribution[μ, σ], x],{x,xmin,xmax}]] (* method 1 *)
(1/(σ (Sqrt[2 π])))(∫_xmin^xmax)(E^(-(((x-μ))^2)/(2(σ^2))))\[DifferentialD]x   (* method 2... Pogo's preferred way of thinking about it *)
(1/2)(Erf[((xmax-μ)/(σ (Sqrt[2])))]-Erf[((xmin-μ)/(σ (Sqrt[2])))])   (* method 3 *)

There are some obvious truths about probability:
1. The probability that the next measurement will be exactly some value is always zero, since any definite integral from x1 to x1 is always zero.
2. The probability that the next measurement will lie between -∞and +∞ is 1.00 (or, 100 percent).
3. If use the letter “p” for probability, and p(a,b) to mean the probability that a measurement will be between a and b, then:
     p(a,b) = p(a,c) - p(b,c).  With these subscripts, I'm sorta assuming a < b < c, but it doesn't really matter!

a=2.8;b=3.2;c=3.5; μ=3.00;σ=0.22;  (* recopied μ and σ for no good reason *)
paa=N[Integrate[PDF[NormalDistribution[μ, σ], x],{x,a,a}]];
pinfinf=N[Integrate[PDF[NormalDistribution[μ, σ], x],{x,-99999.9,99999.9}]]; (* close to ∞! *)
pab=N[Integrate[PDF[NormalDistribution[μ, σ], x],{x,a,b}]];
pbc=N[Integrate[PDF[NormalDistribution[μ, σ], x],{x,b,c}]];
pac=N[Integrate[PDF[NormalDistribution[μ, σ], x],{x,a,c}]];
Print["paa is ",paa]
Print["p-∞ to + ∞ is ",pinfinf]
Print["pab is ",pab]
Print["pbc is ",pbc]
Print["pac is ",pac]
Print["pac - pbc is ",pac-pbc, " which should match pab"]

Complex Numbers

Mathematica has good complex number utility. As usual, Mathematica generally doesn’t want to give decimals if it can avoid them:

Print["4\!\(\*SuperscriptBox[\(E\), \(\*FractionBox[\(π\), \(2\)]\\ I\)]\) reduces to ",4 (E^((π/2) I))]
Print["3\!\(\*SuperscriptBox[\(E\), \(π\\ I\)]\) reduces to ",3 (E^(π I))]
Print["As numbers, 5\!\(\*SuperscriptBox[\(E\), \(\*FractionBox[\(π\), \(4\)]\\ I\)]\) reduces to ",5 (E^((π/4) I))//N]

In this final case, I wanted the “number” form, so I used the N[ ] function. As you see, I used an alternative syntax: If you follow an input line with //, and then the name of a function (here, it’s “N”), it is equivalent to using the function. So, these two are the same:

Sin[0.6]
0.6//Sin

You can also easily do simplifications with complex numbers:

Print["\!\(\*SuperscriptBox[FractionBox[\((1 + 7  I)\), \(4 - I\)], \(3\(\\ \)\)]\) reduces to ",((((1+7I))/(4-I))^(3 ))//N]
Print["(5(E^(47  I))\!)^3 reduces to ",(((5(E^(47I))))^3)//N]

Series

In our first week, we reviewed some stuff about series. For example, we saw that sin(x) was:

f3[x_]=(∑_(i=0)^20)(((((-1))^i) (x^(2 i +1)))/((2*i+1)!));  (* but in class, the upper limit was actually infinity! *)

I got the Σ “summation” object from the Basic Math Assistant Palette, and you can see that I made an approximation that infinity is about the same as 20. It makes you wonder... how big is big enough to be virtually infinite? Let’s compare some plots of this function and the “true” sin(x):

plotrange=15;
plot0=Plot[Sin[x],{x,0,plotrange},PlotStyle->Directive[RGBColor[1,0,0],AbsoluteThickness[6]]];
plot1=Plot[f3[x],{x,0,plotrange},PlotStyle->Directive[RGBColor[0,0,1],AbsoluteThickness[2]]];
Show[plot0,plot1]  (* red is "true" sin (x), blue is series version *)

Here, I overlaid the two plots. The red one is the “true” sin(x), the blue one is the approximation, which actually looks pretty good. I made the red one fatter so that it would be visible beneath the blue one.
While the approximation looks great over this range, we see problems if we look at a larger range:

plotrange=25;
plot0=Plot[Sin[x],{x,0,plotrange},PlotStyle->Directive[RGBColor[1,0,0],AbsoluteThickness[6]]];
plot1=Plot[f3[x],{x,0,plotrange},PlotStyle->Directive[RGBColor[0,0,1],AbsoluteThickness[2]]];
Show[plot0,plot1]

Oops! That starts to look really bad around x = 16 or so. One of many possible ways to quantify the “quality” of our decision to approximate “infinity = 20” is to just look at the difference between the two functions directly:

diff[x_]=f3[x]-Sin[x];

plotrange=10;
Plot[diff[x],{x,0,plotrange},PlotStyle->Directive[RGBColor[1,0,0],AbsoluteThickness[2]],PlotRange-> Full,Frame->True]

It’s a shame that you have to go to the trouble to ask Mathematica to show you the whole plot using the PlotRange->Full command; I don’t think a single other math package in existence forces the user to do this. So, you have to watch out: by default, Mathematica will try to hide some answers from you. 
Look at the vertical axis scale... you notice that assuming infinity is about 20 results in errors that are always smaller than (10^(-10)), at least when x< 10. That’s not bad at all! 
We could make a similar plot over a larger range... specifically, looking at x greater than 16, where the obvious disconnect happened above.

plotrange=25;
Plot[diff[x],{x,0,plotrange},PlotStyle->Directive[RGBColor[1,0,0],AbsoluteThickness[2]],PlotRange-> Full,Frame->True]

Now look! The maximum errors are now around 10 million! And that’s for a function whose actual value never exceeds one!
Conclusion: for this function, using infinity = 20 is Ok if x<10, but is a bad idea if x is larger than 15. Or, we could just try a larger value for "infinity", like perhaps 40 :

f4[x_]=(∑_(i=0)^40)(((((-1))^i) (x^(2 i +1)))/((2*i+1)!));
diff[x_]=f4[x]-Sin[x];
plotrange=25;
Plot[diff[x],{x,0,plotrange},PlotStyle->Directive[RGBColor[1,0,0],AbsoluteThickness[2]],PlotRange-> Full,Frame->True]

OK, based on the vertical scale, that looks pretty good again.
So, for sin(x), using infinity = 40 is not too bad, at least until x exceeds 25.
All those error warning are telling you that using 40 terms is silly anyway... starting at around term #32, the little bits of the series that are being added are becoming smaller than the rounding error incurred because of the way computers represent numbers internally. These warnings are Mathematica’s way of telling you that you were foolish to use 40 terms in the first place, because your accuracy is actually going down as you use more terms!

Matrices

Matrices are a huge weakness of Mathematica! A matrix is the same thing as an array. Internally, Mathematica treats all arrays as one-dimensional, and wants you to think of them that way, too. However, that’s pretty inconvenient. You can make some matrices using the Basic Math Assistant Palette (which also shows you a very complex “hotkey” to use, if you want):

A=(({{2, 3, 4}, {1, -1, -3}, {1, 0, 2}}))
B=(({{4}, {2}, {3}}))
A  B  (* i.e., try to multiply these two matrices... *)

As you can see, these results are beyond ugly. Also, the matrix multiply gave me an error!
It turns out that Mathematica wants you to use a dot (a period) for multiplying matrices. 
You can force Mathematica to at least display results like a grown-up with the following command:

A//MatrixForm
B//MatrixForm
A.B//MatrixForm

Warning! Don’t actually use MatrixForm to convert a matrix into a new form, just use it for displays!
If you change what it is, then it won’t work right in subsequent algebra anymore!

A2 = MatrixForm[A];
Print["Here is true A times B: ", A.B]
Print["Here is MatrixForm A times B: ", A2.B]
Print["Here is true A times B displayed in MF: ", A.B//MatrixForm,"   (this one is nice!)"]
Print["Here is MatrixForm A times B displayed in MF: ", A2.B//MatrixForm]