Physics 229& 100 Homework #2
Name: Nasser Abbasi

1. Find all real solutions of  (e^x)=5 (x^5)-30 (x^3)+4.  Provide convincing commentary in a Text cell which explains why you have found all of them.

Since this is not a polynomial equation we need to use FindRoot[], which requires an initial guess. So we do a plot first to see where can give the guess at.

Remove["Global`*"];
eq1=(E^x);
eq2=5(x^5)-30(x^3)+4;
Plot[eq1-eq2,{x,-3,3}]
Plot[eq1-eq2,{x,-3,16}]

Now looking at the plots above, we can determine the roots.

r1= FindRoot[eq1==eq2,{x,-3}]
r2=FindRoot[eq1==eq2,{x,0}]
r3=FindRoot[eq1==eq2,{x,3}]
r4=FindRoot[eq1==eq2,{x,15}]

The above are all the real roots found. I used 2 methods to convince myself of this. 

Method 1: First by looking at the plot, we see that there can be no more real roots after x=15 since for larger x the exponential terms dominates and increases much faster than the eq2 function which has only x^5 as the largest term. So they not cross again. This can also be seen by plotting for larger and larger x. So the function must continue to increases, so it can not cross the x-axis again. For larger and larger x in the negative direction, the exponential term (eq1) decays rapidly to zero, while x^5 term become the dominant term and it grows larger and larger. Hence eq1 and eq2 curves can never meet again.

Method 2: expend Exp[x] in Taylor series for very large number of terms. This made the terms eq1-eq2 a polynomial equation. Then use Solve. Then count the number of real roots found. I got 4. This below is an example.

nTermsInTaylor=50; 
eq1=Normal[Series[Exp[x],{x,0,nTermsInTaylor}]];
sol=N[x/.Solve[eq1-eq2==0,x],16];
Select[sol,Element[#,Reals]&]

2.  a) Consider the function f(z)= (∫_( 0)^z)(x^x)dx. Construct a plot of this function for x in the range (0,5).

Remove["Global`*"];
f[z_?NumericQ]:=NIntegrate[(x^x),{x,0,z}]
Plot[f[z],{z,0,5},PlotLabel->"Area under (x^x) from 0 to 5",AxesLabel->{"x","Area"},PlotRange->All]

b) For what value of w is f(w)=250?

sol=z/.FindRoot[f[z]==250,{z,4}];
area=N[f[sol],16];
Print["At x=",sol," The area is ",area];

c) If values of f(w) were required repeatedly in a longer calculation, it would be wasteful to recompute the value of the integral each time. It would be more efficient to make a table (using Table) of values, and then to use Fit to find a smooth representation for the values. Find such a fit ( a polynomial in (e^x) works reasonably well) and show that the fit is within ±0.5 of the exact value by constructing an error plot .

del=0.005;
xcoord=Range[0,5,del];
xycoord=N[Table[{w,f[w]},{w,0,5,del}],16];

g=Fit[xycoord,{1,Exp[x],Exp[2 x],Exp[3x],Exp[4x],Exp[5x]},x];

Print["Fitting function is ",g];

error=Table[{xcoord[[i]],Abs[xycoord[[i,2]]-(g/.x->xcoord[[i]])]},{i,Length[xcoord]}];
maxError=Max[error[[All,2]]];

ListPlot[error,PlotRange->All,
PlotJoined->True,
PlotLabel->"Error in fitting. Max error ="<>ToString[maxError],
AxesLabel->{"x","abs error"}
]

3. The sensitivity of the roots of a polynomial to  the exact value of its coefficients is illustrated by a famous pathological polynomial. This polynomial has -1,-2,-3...-20  as roots. 
a) Construct this polynomial using Product and Expand. Verify that the coefficient of the (x^(n-1))term is the negative of the sum of the n=20 roots. (use Coefficient and Sum).  Verify that the constant term is the product of the roots. These are general properties of polynomials.

Remove["Global`*"];
roots=Table[r,{r,-1,-20,-1}];
p=Expand[Product[x-r,{r,-1,-20,-1}]]

c=CoefficientList[p,x];
n=20;
Print["Coefficient of (x^19) term is ",c[[n]]]
Print["negative Sum of the first 20 roots is ",-Total[roots]]

Print["Constant term is ",c[[1]]];
Print["Product of all roots is ",Product[roots[[i]],{i,1,Length[roots]}]]

b) Verify that NSolve finds all the roots of the polynomial. Now construct a new polynomial by adding 1.0 *(10^(-6 )) to the coefficient of (x^19). Use NSolve to find the roots of the new polynomial. Make a plot of their location in the complex plane by constructing a separate list of the real and imaginary parts, and then combining them into {x,y} pairs using Transpose. Make a plot using ListPlot.

solBefore=NSolve[p==0,x]

(*slightly change the coefficient for the (x^19) term*)

delta=(10^(-6));
p[[20,1]]=p[[20,1]]+delta;
solAfter=x/.NSolve[p==0,x]

Notice how a very small perturbation in one term coefficient made very large change in solution

x=Re[solAfter];
y=Im[solAfter];
pts=Transpose[{x,y}];
ListPlot[pts,PlotStyle->PointSize[.02],PlotLabel->"roots of perturbated polynomial"]

4. Consider a uniformly charged disk of radius R. Find the electrostatic potential at a point z on the axis; call it potz (set up the appropriate integral using Integrate). Find a series expansion for potz valid for large z. Do this in two different ways:

a) First, make a substitution z->1/ϵ; in this new formula, ϵ->0 is equivalent to z->∞. Find the series expansion in ϵ about zero, and then substitute ϵ->1/z.

Note: a differential area on the disk is da=r dθ dr

Remove["Global`*"]

d=(Sqrt[(r^2)+(z^2)]);
potz[z_]:=(σ/(4π ϵ))Assuming[Element[{R0,r},Reals]&&R0> 0&&z> 0,Integrate[Integrate[(r/d),{theta,0,2 Pi}],{r,0,R0}]];
sol=potz[z]

Print["Integral is " ,solutionInZ=sol];

solutionInZeta=solutionInZ/.z->(1/ζ);
seriesInZetaAboutZero=Series[solutionInZeta,{ζ,0,10}]

seriesInZAboutZero=Normal[seriesInZetaAboutZero]/.ζ->(1/z)
Print["Static potential at large z is ",Limit[Normal[seriesInZAboutZero],z->∞]]

b) Second, expand potz about the point ∞ (see Series). Confirm that you get the same answer both ways. (the second way is obviously easier).

seriesInZAboutInfinity=Series[sol,{z,Infinity,10}]

(*Now verify the same answer as part a*)
Normal[seriesInZAboutZero]-Normal[seriesInZAboutInfinity]

Limit[seriesInZAboutInfinity,z->∞]

c) If we measure z in units of the disk radius R, we can set R->1; Plot potz and the first two terms of the series expansion on the same graph (set σ and any other constants->1) . Show that the series expansion is a terrible approximation at z=1/2, but is excellent at z=2. The first two terms are called the monopole and the quadrupole terms in the potential. Check your calculation by making sure that the monopole term is the potential of the total charge concentrated at the origin.

(*Find potential due to total charge concentrated at center*)
pe=(1/(4 π ϵ))((π (R^2)σ)/z);
Print["Potential due to total change at center=",pe];
Print["Monopole term in series is ",(1/z)Coefficient[Normal[seriesInZAboutZero],1/z]];

plotit[numberOfTerms_]:=Module[{p},
p=Plot[ Evaluate[{potz[z],
Normal[Series[potz[z],{z,Infinity,numberOfTerms}]]
}/.{σ->1,R0->1,ϵ->1}
],
{z,0,2.5},
PlotLabel->"Compare with series expanded at infinity using "<>ToString[numberOfTerms]<>" terms",
PlotStyle->{Red,Blue},
AxesLabel->{"z",""}
]
];

numberOfTerms=10;
plotit[numberOfTerms]

We see from the plots above the following: The BLUE curve represents potz expressed in series expansion around infinity, while the RED curve is potz function. We see that for z less than around .5, the approximation is very bad. The approximation improves quickly after about z=1. At z=2 for example, one sees no difference between the series expansion and potz itself.

5.

Find an expression for the voltage at point d in the circuit above. Go here for a reminder about rules for circuits. Clearly describe what variables you are using and what principles you use to generate your equations.  A good way to do this is with a drawing.You should cut and paste the drawing above into a graphics package  such as Paint that works with bit map graphics and add current labels and arrows.

The above diagram shows the method to use for generating the equations. We have 3 loops. Using Kirchoff laws, we know that the sum of voltages in each loop must add to zero. Using this we generate 3 equations (one per loop), which we solve for I1,I2,I3. Next we use Ohm's law to find Voltage at point d.

Remove["Global`*"];
loop1=V-I1 R1+I2 R1 -I1 R3+ I3 R3==0;   
loop2=-I2 R2 -I2 R5 +I3 R5 -I2 R1 +I1 R1==0;
loop3=-I3 R5 +I2 R5 -I3 R4 -I3 R3 +I1 R3==0;
sol=Solve[{loop1,loop2,loop3},{I1,I2,I3}]//Simplify//Flatten

Now we can find the voltage at point d

Vd=(V-I2 R2)/.sol//Simplify

6. Graphics exercises. 
a) Use Mathematica to make a drawing  of unit vectors emanating from the origin that looks like

Use Arrow, which  requires loading a package, as described here.  Arrow is similar to Line; the arguments are Arrow[ {tail},{head}] , where head and tail are vectors that specify the location of the head and tail of the arrow.  Use Table to make a list of Arrows.

Remove["Global`*"];

nArrows=32; 
delAngle=((2π)/nArrows);

endCoordinates[r_,angle_]:=Module[{c},
c=CoordinatesToCartesian[{r,angle,0},Cylindrical];(*No polar? used cyl with z=0*)
{c[[1]],c[[2]]}
]

arrows=Table[Arrow[{0,0},endCoordinates[1,n delAngle]],{n,0,nArrows-1}]
Show[Graphics[arrows],AspectRatio->Automatic]

b) Make another drawing which adds a small red circle at the end of each arrow . Information about graphics directives to change properties such as color are described here.

circles=Table[Circle[endCoordinates[1.1,n delAngle],.1],{n,0,nArrows-1}];
Show[Graphics[arrows],Graphics[{Thickness[.01],RGBColor[1,0,0],circles}],AspectRatio->Automatic];

c) The equation (x^3)-(y^3)= 6 x y  defines a curve in the plane. Plot this curve using either ContourPlot or ImplicitPlot as described here.

Needs["Graphics`ImplicitPlot`"]
Remove[x,y,eq];
eq=(x^3)-(y^3)==6x y;

ImplicitPlot[eq,{x,-5,5},
FrameLabel->eq,
AxesLabel->{"x","y"},
Frame->True,
PlotStyle->RGBColor[1,0,0]
];

d) Use ContourPlot3D to make a plot of the unit sphere.

Remove["Global`*"];
Needs["Graphics`ContourPlot3D`"];

eq=(x^2)+(y^2)+(z^2)-1;
ContourPlot3D[eq,{x,-1,1},{y,-1,1},{z,-1,1},PlotLabel->"Unit sphere in 3D"];

e) Use ParametricPlot3D to make a plot of the unit sphere. The {x,y,z} coordinates of the points on a sphere are described here.

Remove["Global`*"];
Needs["Graphics`ParametricPlot3D`"];

(*use the parameters ϕ,θ to parameterize the coordinates x,y,z based on spherical coordinates mapping*)

r=1;
x=r Cos[ϕ]Sin[θ]; y=r Sin[θ]Sin[ϕ]; z=r Cos[θ];
ParametricPlot3D[{x,y,z},{ϕ,0,2π},{θ,0,π}];

f) Construct a drawing of 4 mutually tangent spheres centered on the vertices of a tetrahedron.  Use ParametricPlot3D to draw the spheres. Turn on RealTime3D so the drawing can be rotated with the mouse.

Remove["Global`*"];
Needs["Graphics`ParametricPlot3D`"];
Needs["Graphics`Polyhedra`"];
<<Default3D`;

nVertices=4;  v=Vertices[Tetrahedron];

plotAt[{xx_,yy_,zz_}]:=Module[{r,x,y,z,θ,ϕ},
r=1.4;(*I make r little smaller so spheres do not touch each others*)
                  (*So that we can see the Tetrahedron. Make r=1.4 otherwise*)
x=xx+r Cos[ϕ]Sin[θ]; y=yy+r Sin[θ]Sin[ϕ]; z=zz+r Cos[θ];
p=ParametricPlot3D[{x,y,z},{ϕ,0,2π},{θ,0,π},DisplayFunction->Identity]
];

p=Table[plotAt[v[[i]]],{i,1,nVertices}];

(*First show it in default 3D*)

Show[{Polyhedron[Tetrahedron,{0,0,0},1],Map[p[[#]]&,Range[1,nVertices]]},DisplayFunction->$DisplayFunction];

<<RealTime3D`;

Show[{Polyhedron[Tetrahedron,{0,0,0},1],Map[p[[#]]&,Range[1,nVertices]]},DisplayFunction->$DisplayFunction];

<<Default3D`

7. Computational algorithms are said to work in polynomial time if the execution time is proportional to (N^n), where N is a characteristic scale of the length of the input data, and n is a fixed constant that characterizes the efficiency of the algorithm. Using the discussion of the computing time for Det as a guide, estimate n for the Mathematica function Eigenvalues. How long would it take to find the eigenvalues of a 3000 x 3000 matrix on your computer? ( Don't actually try to do it).

Remove["Global`*"];

makematrix[n_]:=Table[Random[Real,{-10,10}],{i,n},{j,n}];
TimevsN=Table[{n,First[Timing[Eigenvalues[makematrix[n]]]]/Second},{n,100,800,20}];

ListPlot[TimevsN,
PlotJoined->True,
PlotLabel->"Time in seconds vs N",
AxesLabel->{"N","t (sec)"}];

{times,matsize}=Transpose[TimevsN];

LogtvsLogn=Transpose[{Log[times],Log[matsize]}];

ListPlot[LogtvsLogn,
PlotJoined->True,
PlotLabel->"Log[t]vs Log[N]",
AxesLabel->{"Log[N]","Log[t]"}
];

y=Fit[LogtvsLogn,{1,x},x];
Print["Fitted line is ",y];

slope=Coefficient[y,x];

Print["Execution time propertional to N^",slope]

nElements=3000.;
timeWillTake=y/.x->Log[nElements];
Print["It will take ",Exp[timeWillTake]," seconds to compute eigenvalues of 3000x3000 matrix"]

8 a) Calculate the exact value of the sum 1-(1/4)+(1/7)-(1/10)+(1/13)...  ( Hint: a convenient way to find an expression for the nth term is to use Fit and a linear function of n  on the terms you know)

By inspection the general term is 1/(((-1))^n)(4n-(n-1)) for n=0..∞

Remove["Global`*"];
term[n_]:=(((-1))^n)(1/(4n-(n-1)))
Table[term[n],{n,0,10}]

(*lets make sure the sequence converges*)
Limit[term[n],n->Infinity]

Print["Exact sum is =", exactSum=Sum[term[n],{n,0,Infinity}]]
Print["Which is equal to ",N[exactSum,16]]

write a function that generates the partial sum of the first Nt terms of the series. Make a plot of the value of the partial sums as a function of Nt  from 1 to 50.

partSum[upTo_]:=Sum[term[n],{n,0,upTo}]
pts=Table[{n,partSum[n]},{n,1,50}];
ListPlot[pts,PlotJoined->True,PlotRange->All,PlotLabel->"Sum as function of n",AxesLabel->{"N","Sum"}];

How many terms of the sum are required to ensure that the error is less than 0.001 ?

minError=0.001;
nTerms=0;
keepLooking=True;

While[keepLooking,
nTerms=nTerms+1;
error=Abs[ N[ exactSum-partSum[nTerms],16]]; 
If[error<minError,keepLooking=False]
];

Print["Number of terms needed is ",nTerms];

9

Use MonteCarlo integration using Method->MonteCarlo  in NIntegrate to compute the volume of high dimensional spheres.
a) first verify that you get the right answer for a 3 dimensional sphere of radius =1

First I will use spherical coordinates for 3D. Using spherical coordinates, the volume is found by integrating a unit volume. We obtain the formula that V= 8(∫_0^(π/2))(∫_0^(π/2))(∫_0^R)((r^2)Sin[ϕ])drdϕdθ

Remove["Global`*"];
R=1;
numericalAnswer=8*NIntegrate[(r^2)Sin[ϕ],{r,0,R},{ϕ,0,(π/2)},{θ,0,(π/2)},Method->MonteCarlo];
exactAnswer=(4/3)π (R^2);
Print["Using NIntegrate witg MonteCarlo we obtain V= ",numericalAnswer,"\nUsing the standard formula the answer is V= ",N[exactAnswer,16]];

Using the normal Cartesian coordinates systems,we get

numericalAnswer=NIntegrate[1,
{x1,-R,R},
{x2,-Sqrt[(R^2)-(x1^2)],Sqrt[(R^2)-(x1^2)]},
{x3,-Sqrt[(R^2)-(x1^2)-(x2^2)],Sqrt[(R^2)-(x1^2)-(x2^2)]},
Method->MonteCarlo
];
Print["Using NIntegrate with MonteCarlo we obtain V= ",numericalAnswer];

NIntegrate[If[x1^2+x2^2+x3^2<1,1,0],{x1,-1,1},{x2,-1,1},{x3,-1,1},Method->MonteCarlo,PrecisionGoal->3,MaxPoints->1000000]

b) Next, calculate the volume of a sphere in 7 dimensions with radius =1

To do this in higher dimensions, I will use the normal Cartesian coordinates systems.

(*Now generialize to n=7*)
numericalAnswer=NIntegrate[1,
{x1,-R,R},
{x2,-Sqrt[(R^2)-(x1^2)],Sqrt[(R^2)-(x1^2)]},
{x3,-Sqrt[(R^2)-(x1^2)-(x2^2)],Sqrt[(R^2)-(x1^2)-(x2^2)]},
{x4,-Sqrt[(R^2)-(x1^2)-(x2^2)-(x3^2)],Sqrt[(R^2)-(x1^2)-(x2^2)-(x3^2)]},
{x5,-Sqrt[(R^2)-(x1^2)-(x2^2)-(x3^2)-(x4^2)],Sqrt[(R^2)-(x1^2)-(x2^2)-(x3^2)-(x4^2)]},
{x6,-Sqrt[(R^2)-(x1^2)-(x2^2)-(x3^2)-(x4^2)-(x5^2)],Sqrt[(R^2)-(x1^2)-(x2^2)-(x3^2)-(x4^2)-(x5^2)]},
{x7,-Sqrt[(R^2)-(x1^2)-(x2^2)-(x3^2)-(x4^2)-(x5^2)-(x6^2)],Sqrt[(R^2)-(x1^2)-(x2^2)-(x3^2)-(x4^2)-(x5^2)-(x6^2)]},
Method->MonteCarlo
];

Print["Using NIntegrate with MonteCarlo we obtain V for n=7 as  ",numericalAnswer];

c)( grads)  write a function Vsphere[ndim_] that will use montecarlo integration to calculate the volume of the unit n-sphere in ndim dimensions. Make a table of values and a plot for dnim={2,12}, which shows that  as a function of n, a 5-sphere has the maximum volume and that ( counterintuitively) the volumes of spheres for dimensions >5 decreases rapidly with dimension. You may have to experiment with PrecisionGoal and MaxPoints.
To automatically generate the limits for your function, you will want to make a list of limits called lims, and then use Evaluate[Apply[Sequence,lims]] to splice the arguments into NIntegrate. See Further Examples in the Help for Sequence.

NOTE: I am attempting this as extra, I am not signed as a Grad

Using a formula based on Gamma function for finding sphere Volume in higher n, we can generate the above function.

v[n_]:=(((π^(n/2))(R^n))/(Gamma[(n/2)+1]))

pts=Table[{n,v[n]},{n,2,12}];

p1=ListPlot[pts,
PlotStyle->PointSize[.03],
PlotRange->{{0,14},{0,6}},
AxesOrigin->{0,0},
DisplayFunction->Identity
];

p2=Graphics[Line[pts]];

Show[{p1,p2},
DisplayFunction->$DisplayFunction,
PlotLabel->"Volume of sphere as N changes",
AxesLabel->{"ndim","Volume"}
];