Code for Electronic Supplement

FindRoots

(* This version accepts options for FindRoot *) 

FindRoots[f_,{x_,a_,b_},opts___]:=Module[{p,seeds},p=Plot[f,{x,a-((b-a)/100),b+((b-a)/100)},Mesh->{{0}},MeshFunctions->(f/.x->#1&),Evaluate[Sequence@@FilterRules[{opts},Options[Graphics]]]];seeds=Cases[Normal[p],Point[{z_,_}]:>z,∞]⋃(SameTest->(Abs[#1-#2]<(1/(10^6))&));Select[If[seeds=={},{},Table[Re[x]/.FindRoot[f==0,{x,seed},Evaluate[Sequence@@FilterRules[{opts},Options[FindRoot]]]],{seed,seeds}]⋃(SameTest->(Abs[#1-#2]<(1/(10^6))&))],a≤#1≤b&]]

FindRoots2D

FindRoots2D::usage=
"FindRoots2D[funcs,{x,a,b},{y,c,d}] finds all nontangential solutions to {f=0, g=0} in the given rectangle.";

Options[FindRoots2D]={PlotPoints->Automatic, MaxRecursion->Automatic};

FindRoots2D[funcs_,{x_,a_,b_},{y_,c_,d_},opts___]:=Module[{fZero,seeds,signs,fy},
fy=Compile[{x,y},Evaluate[funcs[[2]]]];
fZero=Cases[Normal[ContourPlot[funcs[[1]]==0,{x,a-((b-a)/97),b+((b-a)/103)},{y,c-((d-c)/98),d+((d-c)/102)},Evaluate[FilterRules[{opts}, Options[ContourPlot]]]]], Line[z_] :> z,∞];
seeds=Flatten[(signs=
Sign[Apply[fy,#,{1}]];#1[[1+Flatten[Position[Rest[signs * RotateRight[signs]],-1]]]])&/@fZero,1];
If[seeds=={},{},Select[Union[
({x,y}/.FindRoot[{funcs[[1]],funcs[[2]]},{x,#[[1]]},{y,#[[2]]},Evaluate[FilterRules[{opts}, Options[FindRoot]]]]&)/@seeds, SameTest -> (Norm[#1-#2]<10^-6&)],a≤#[[1]]≤b&&c≤#[[2]]≤d&]]]

Solving Equations

Code for Chapter Opening Art

Clear[f,g];
f[x_,y_]:=-Cos[y]+2y Cos[(y^2)]Cos[2x];
g[x_,y_]:=-Sin[x]+2Sin[(y^2)]Sin[2x];
roots=FindRoots2D[{f[x,y], g[x,y]},  {x, -3.5, 4}, {y, -1.8, 4.2}];
Show[ContourPlot[{f[x,y]==0,g[x,y]==0},{x,-3.5,4},{y,-1.8,4.2},ContourStyle->{{Thick,Black,Red},{Thickness[0.004],Blue}}],Graphics[{EdgeForm[Thickness[.003]],Yellow,(Disk[#1,0.07]&)/@roots}],FrameTicks->{Range[-3,4],Range[-1,4],None,None},AspectRatio->Automatic,ImageSize->400]

.  Solve

Solve[(x^3)+2 (x^2)-1==0,x]

x /. Solve[(x^6)+2(x^5)+2(x^3)+3(x^2)-4x-4==0]

Solve[(x^7)+x-1==0]

x /. N[%]

Solve[x Log[x]==2]

Reduce[x Log[x]==2,x]

Plot[LambertW[x],{x,-(1/E),3}]

Series[LambertW[x],{x,0,7}]

D[LambertW[x],x]

∫LambertW[x]\[DifferentialD]x

ContourPlot[{(x^2)+(y^3)-2==0,(x^2)+((y^2)/4)-2==0},{x,-2,2},{y,-3,3},ContourStyle->{{Red,Thick},{Blue,Thick}}]

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

ContourPlot[{(x^2)+(y^3)-2,(x^2)+((y^2)/4)-2},{x,1.4,1.42},{y,-0.2,0.3},ContourStyle->{{Red,Thick},{Blue,Thick}}]

Reduce[13 (x^2)<10&&5 x+2>0]

sol=Reduce[(x^3)+13(x^2)<10&&5 x+2>-0]

ToRadicals[sol[[5]]]

N[sol]

Reduce[Sin[x]==Cos[x],x]

s=Reduce[Sin[x]==Cos[x]+Tan[x],x, GeneratedParameters->(n&)]

InputForm[Root[1-2 ((#1)^2)+4 ((#1)^3)+((#1)^4)&,4]]

x/.{ToRules[s[[2]]]}/.Root[u1_,u2_,u3_]:>N[Root[u1,u2,u3]]

.  Diophantine Equations

Solve[{7x==13, Modulus==23}, x]

x/.Solve[{7 (x^5)+13 x+13==0,Modulus==541},x]

Reduce[31 a+7 b==91,{a,b},Integers,GeneratedParameters->(n&)]

Reduce[(x^2)-7 (y^2)==1,{x,y},Integers,GeneratedParameters->(n&)]

{x,y}/. {ToRules[Reduce[{(x^2)-7 (y^2)==1,0≤x≤(10^6),0≤y≤(10^6)},{x,y},Integers]]}

n/. {ToRules[Reduce[{EulerPhi[n]==1000,n≥0},n, Integers]]}

PhiMultiplicity[d_]:=Length[Reduce[{EulerPhi[n]==d,n≥0},n,Integers]]
Union[PhiMultiplicity/@Range[0,5000]]

Reduce[13x+17y+21z==234 && And @@ Thread[{x,y,z}≥0], {x,y,z}, Integers]

FrobeniusSolve[{13,17,21},234]

FrobeniusSolve[{13,17,21},234,1]

FrobeniusSolve[{132345,172345235,2234523451},2323452345234523454,1]

FrobeniusNumber[{6,9,20}]

FrobeniusNumber[(10^10)+Prime[Range[5]]]

FrobeniusNumber[10000+Prime[Range[20]]]

PowersRepresentations[1729,2,3]

PowersRepresentations[1729,3,2]

SquaresR[20,500]//Timing

SetSystemOptions["ReduceOptions"->{"SieveMaxPoints"->10000000}];
FindInstance[31==(x^3)+(y^3)+(z^3)+(w^3),{x,y,z,w},Integers]

FindInstance[(10^100)+39==(x^2)+(y^2)+(z^2)+(w^2),{x,y,z,w},Integers]

.  LinearSolve

A={{1,2,3},{4,5,6},{2,3,1}};b={1,2,1};
Solve[A.{x,y,z}==b]

LinearSolve[A, b]

A=SparseArray[{{1,2}->5}, 5]

Normal[A]

MatrixForm[A]

n=10;
A=(SparseArray[{{i_,i_}->Prime[i]},n]+(#+Transpose[#])&)[SparseArray[Flatten[Table[{i,i+(2^j)}->1,{i,n-1},{j,0,Log[2,n-i]}]],n]]

MatrixForm[A]

b=ConstantArray[0,n];  b[[1]]=1;
LinearSolve[A,b][[1]]

n=20000;  b=ConstantArray[0,n];  b[[1]]=1;
A=(SparseArray[{{i_,i_}->Prime[i]},n]+(#+Transpose[#])&)[SparseArray[Flatten[Table[{i,i+(2^j)}->1,{i,n-1},{j,0,Log[2,n-i]}]],n]];

LinearSolve[A,N[b],Method->{"Krylov",Method->"ConjugateGradient"}][[1]]

diagonal=Table[A[[i,i]],{i,n}];

LinearSolve[A,N[b, 40],Method->{"Krylov",Method->"ConjugateGradient","Preconditioner"->(((#)/diagonal)&)}][[1]]

.  NSolve

x/.NSolve[(x^11)+(x^7)-3 (x^2)==0,x]

NSolve[(E^x)==(x^10),x]

Reduce[(E^x)==(x^10),x, GeneratedParameters->(n&)]

N[Reduce[(E^x)==(x^10),x,Reals]]

AllRoots[deg_,cmax_]:= (
coeffs=Tuples[Range[-cmax,cmax], deg+1];
polys=Table[c . (z^(Range[0,deg])),{c,coeffs}];
Quiet[Cases[Flatten[Table[{Re[z],Im[z]}/.NSolve[p,z],{p,polys}],1],_?VectorQ]]);

AllRootsUpToDegree[dmax_,cmax_,opts___]:= (roots=Join @@ Table[AllRoots[d,cmax], {d, 1,dmax}];
Graphics[{PointSize[0.002],Point@roots}, Sequence@@FilterRules[{opts},Options[Graphics]],Frame->True, PlotRange->{{-4,4}, {-2.5, 2.5}},FrameTicks->{Automatic,{#, #I} & /@ Range[-2,2],None,None}]);

AllRootsUpToDegree[3,4]

.  FindRoot

θ/.FindRoot[((Sin[θ])/θ)==(5280/5281),{θ,0.2}]

obj[y_Real] :=NIntegrate[Sin[Sin[x]], {x,0,y}]
FindRoot[obj[y]==0.2, {y, .1}]

x/.FindRoot[(x^2)==2,{x,1},AccuracyGoal->20,WorkingPrecision->30]

sol={x,y}/.FindRoot[{Sin[x y]==0.5,Cos[x+(y^2)]==0.6},{x,-1},{y,2.5}]

ContourPlot[{Sin[x y]==0.5,Cos[x+(y^2)]==0.6},{x,-2,0},{y,0,3},ContourStyle->{{Red,Thickness[0.01]},{Blue,Thickness[0.01]}},Epilog->{PointSize[Large],Point[sol]}]

.  Finding All Roots in an Interval

Clear[f];f[x_]:=Sin[(x^2)]-(E^((Cos[x])/2))+1;
p=Plot[f[x],{x,-10,0},Mesh->{{0}},MeshFunctions->(f[#1]&)];
seeds=Cases[Normal[p],Point[z_]:>z[[1]],∞]

roots=Table[x/.FindRoot[f[x],{x,s}],{s,seeds}]⋃(SameTest->(Abs[#1-#2]<(1/(10^6))&))

Plot[f[x],{x,-10,0},Epilog->{PointSize[0.013],Red,Point[Transpose[{roots,f[roots]}]]},Frame->True,Axes->False]

FindRoots[f_,{x_,a_,b_},opts___]:=Module[{p,seeds},p=Plot[f,{x,a-((b-a)/100),b+((b-a)/100)},Mesh->{{0}},MeshFunctions->(f/.x->#1&),Evaluate[Sequence@@FilterRules[{opts},Options[Graphics]]]];seeds=Union[Cases[Normal[p],Point[{z_,_}]:>z,∞],
              SameTest->(Abs[#1-#2]<(10^(-6))&)];
Select[If[seeds=={},{},
Union[Table[Re[x]/.FindRoot[f==0,{x,seed},Evaluate[
Sequence@@FilterRules[{opts},Options[FindRoot]]]],
                {seed,seeds}],
SameTest->(Abs[#1-#2]<(10^(-6))&)]],a≤#≤b&]]

f[x_]:=(x/5)  Cos[x]+Sin[(x^3)];
solns=FindRoots[f[x],{x,-3,3}]

Plot[f[x],{x, -3, 3},Frame->True,AxesStyle->GrayLevel[0.6],
	 Epilog->{PointSize[0.02],Point[{#,0}&/@solns]}]

Reduce[f[x]==0&&-3<x<3,x][[1]]

N[x/.{ToRules[Reduce[f[x]==0&&-3<x<3,x]]}]

FindRoots[(E^x)-(x^10),{x,-40,40}]

FindRoots[Log[(E^x)]-Log[(x^10)],{x,-40,40}]

FindRoots[RiemannSiegelZ[t],{t,0,100}]//Timing

N[Im[ZetaZero[Range[29]]]]//Timing

.  FindRoots2D

f[x_,y_]:= x - (y^2) Cos[y];
g[x_,y_]:=-y +x Sin[x];
fcnVec[{x_,y_}]:= {f[x,y], g[x,y]};

Short[fZero=Cases[Normal[ContourPlot[f[x,y]==0,{x,-10,10},{y,-10,10}]],Line[z_]:>z,∞],15]

Show[Graphics[Point/@fZero],ContourPlot[g[x,y]==0,{x,-10,10},{y,-10,10},ContourStyle->{Thick,Red}],Frame->True]

(signs=Sign[Apply[g,fZero[[3]],{1}]])//Short

(changes=Rest[signs RotateRight[signs]])//Short

posns=Flatten[Position[changes,-1]]

seeds=fZero[[3,Flatten[posns]]]

roots=Table[v/.FindRoot[fcnVec[v],{v,s}],{s,seeds}]

Show[ContourPlot[{f[x,y]==0,g[x,y]==0}, {x,-10,10}, {y,-10,10}, ContourStyle->{{Thick,Red},{Thick,Black}}],Graphics[{Thickness[0.005],Blue,Table[Circle[r, .3], {r, roots}]}], Frame->True]

FindRoots2D::usage=
"FindRoots2D[funcs,{x,a,b},{y,c,d}] finds all nontangential solutions to {f=0, g=0} in the given rectangle.";

Options[FindRoots2D]={PlotPoints->Automatic, MaxRecursion->Automatic};

FindRoots2D[funcs_,{x_,a_,b_},{y_,c_,d_},opts___]:=Module[{fZero,seeds,signs,fy},
fy=Compile[{x,y},Evaluate[funcs[[2]]]];
fZero=Cases[Normal[ContourPlot[funcs[[1]]==0,{x,a-((b-a)/97),b+((b-a)/103)},{y,c-((d-c)/98),d+((d-c)/102)},Evaluate[FilterRules[{opts}, Options[ContourPlot]]]]], Line[z_] :> z,∞];
seeds=Flatten[(signs=
Sign[Apply[fy,#,{1}]];#1[[1+Flatten[Position[Rest[signs * RotateRight[signs]],-1]]]])&/@fZero,1];
If[seeds=={},{},Select[Union[
({x,y}/.FindRoot[{funcs[[1]],funcs[[2]]},{x,#[[1]]},{y,#[[2]]},Evaluate[FilterRules[{opts}, Options[FindRoot]]]]&)/@seeds],a≤#[[1]]≤b&&c≤#[[2]]≤d&]]]

Clear[f,g];
f[x_,y_]:=-Cos[y]+2y Cos[(y^2)]Cos[2x];
g[x_,y_]:=-Sin[x]+2Sin[(y^2)]Sin[2x];
roots=FindRoots2D[{f[x,y], g[x,y]},  {x, -3.5, 4}, {y, -1.8, 4.2}];
Length[roots]

Show[ContourPlot[{f[x,y]==0,g[x,y]==0},{x,-3.5,4},{y,-1.8,4.2},ContourStyle->{{Thick,Red},{Thick,Blue}}],Graphics[{EdgeForm[Thickness[0.003]],Yellow,(Disk[#1,0.07]&)/@roots}],AspectRatio->Automatic]

.  Two Applications

A Challenging Minimization

f[x_,y_]:=(E^(Sin[50 x]))+Sin[60 (E^y)]+Sin[70 Sin[x]]+Sin[Sin[80 y]]-Sin[10 (x+y)]+(1/4) ((x^2)+(y^2));
Plot3D[f[x,y],{x,0,1},{y,0,1}]

Plot3D[f[x,y], {x, -100, 100}, {y,-100, 100}, RegionFunction->(Norm[{#1,#2}]<100&),MeshFunctions->(#3&)]

Min[Table[f[x,y],{x,-1,1,0.01},{y,-1,1,0.01}]]

N[IntervalUnion[f[Interval[{1,∞}],Interval[{-∞,∞}]],f[Interval[{-∞,-1}],Interval[{-∞,∞}]],f[Interval[{-∞,∞}],Interval[{1,∞}]],f[Interval[{-∞,∞}],Interval[{-∞,-1}]]]]

gradf[x_, y_] :=Evaluate[D[f[x,y], {{x,y}}]];
CPs=Flatten[Table[FindRoots2D[gradf[x,y],{x,i,i+0.5},{y,j,j+0.5}],{i,-1,0.5, 0.5},{j,-1,0.5, 0.5}],2];
Length[CPs]

data=Table[{pp,
CPs=Flatten[Table[FindRoots2D[gradf[x,y],{x,i,i+0.5},{y,j,j+0.5},PlotPoints->pp],{i,-1,0.5, 0.5},{j,-1,0.5, 0.5}],2];
Length[CPs]},{pp,20,70,5}]

Show[ContourPlot[Evaluate[gradf[x,y]],{x,0,0.5},{y,0,0.5},ContourShading->False,PlotPoints->50,ContourStyle->{{Thickness[0.005],Red},{Thickness[0.005],Blue}}],Graphics[({EdgeForm[Black],Yellow,Disk[#,0.0045]}&)/@Select[CPs,0<#[[1]]<0.5&&0<#[[2]]<0.5&]]]

SortBy[Table[{cp,f[cp[[1]],cp[[2]]]},{cp,CPs}],Last][[1]]

Classification and Morse Theory

ClassifyCP2D[f_,{x_,a_}, {y_,b_}]:= Module[{
derxx=(∂_(x,x))f/.{x->a,y->b},disc},
disc=derxx (∂_(y,y))f-((((∂_(x,y))f))^2)/.{x->a, y->b};
Which[disc>0 && derxx<0, "Maximum",
disc>0 && derxx>0, "Minimum",
disc<0,"Saddle",
disc==0,"Unknown"]]

classData=Table[ClassifyCP2D[f[x,y],{x,cp[[1]]},{y,cp[[2]]}],{cp,CPs}];
Table[{type, Count[classData, type]}, {type, {"Minimum", "Maximum", "Saddle"}}]

693+667

d=0.1;
Show[Plot3D[f[x,y],{x,0,.4},{y,0,0.4},MeshFunctions->(#3&), Mesh->{Range[-5,5,10/15]}, BoundaryStyle->BoundaryStyle->{Thick, Black}],
Plot3D[f[0,y] ((((x+d))/d))-5(1-(((x+d))/d)),{x,-d,0},{y,0,0.4},MeshFunctions->(#3&), Mesh->{Range[-5,5,10/15]},BoundaryStyle->None],

 Plot3D[f[x,0] ((((y+d))/d))-5(1-(((y+d))/d)),{x,0,0.4},{y,-d,0},MeshFunctions->(#3&), Mesh->{Range[-5,5,10/15]},BoundaryStyle->None], 
ParametricPlot3D[{r Cos[θ],r Sin[θ],(1-(r/d))f[0,0]-(r/d)5},{θ,π,3π/2},{r,0,d},MeshFunctions->(#3&), Mesh->{Range[-5,5,10/15]},
BoundaryStyle->None],PlotRange->All]

{Sign[gradf[Interval[{-1,1}],-1][[2]]],Sign[gradf[1,Interval[{-1,1}]][[1]]],Sign[gradf[Interval[{-1,1}],1][[2]]],Sign[gradf[-1,Interval[{-1,1}]][[1]]]}

der[1]=(∂_x)f[x,-1];
der[2]=(∂_y)f[1,y];
der[3]=(∂_x)f[x,1];
c1=FindRoots[der[1],{x,-1,1}];
c2=FindRoots[der[2],{y,-1,1}];
c3=FindRoots[der[3],{x,-1,1}];
Length[c1]+Length[c2]+Length[c3]

ClassifyCP1D[f_,{x_,a_}]:=Sign[(∂_(x,x))f/.x->a]/.{-1->"Maximum",1->"Minimum",0->"Unknown"}

Count[c1, xx_ /; ClassifyCP1D[f[x,-1],{x,xx}]=="Maximum"]+
Count[c2, yy_ /;ClassifyCP1D[f[1,y],{y,yy}]=="Maximum"]+
Count[c3, xx_ /;ClassifyCP1D[f[x,1],{x,xx}]=="Maximum"]
Count[c1, xx_ /; ClassifyCP1D[f[x,-1],{x,xx}]=="Minimum"]+
Count[c2, yy_ /;ClassifyCP1D[f[1,y],{y,yy}]=="Minimum"]+
Count[c3, xx_ /;ClassifyCP1D[f[x,1],{x,xx}]=="Minimum"]

N[Table[(gradf@@cor).cor,{cor,Tuples[{-1,1},2]}]]

GraphicsGrid[{{Show[Plot[f[-1,-1+y],{y,-0.01,0},PlotStyle->{Black,Thick}],Plot[f[x-1,-1],{x,0,0.01},PlotStyle->{Red,Thick}],PlotRange->All,Frame->True,Axes->False,PlotLabel->"Southwest"],
Show[Plot[f[x+1,-1],{x,-0.0001,0},PlotStyle->{Black,Thick}],Plot[f[1,y-1],{y,0,0.00001},PlotStyle->{Red,Thick}],PlotRange->All,Frame->True,Axes->False,PlotLabel->"Southeast"]},{Show[Plot[f[1,y+1],{y,-0.0001,0},PlotStyle->{Black,Thick}],Plot[f[1-x,1],{x,0,0.001},PlotStyle->{Red,Thick}],PlotRange->All,Frame->True,Axes->False,PlotLabel->"Northeast"],
Show[Plot[f[-1-x,1],{x,-0.01,0},PlotStyle->{Black,Thick}],Plot[f[-1,1-y],{y,0,0.01},PlotStyle->{Red,Thick}],PlotRange->All,Frame->True,Axes->False,PlotLabel->"Northwest"]}}]