Physics 229& 100 Final Exam 2006
Due Tuesday, Dec 12 at noon 
Open book, but no discussion with anyone except me. You may use any packages that seem helpful. Provide clear commentary so that I can follow your reasoning.
grads do all problems
undergrads do 1-13 ( 14&15 are extra credit if you do them ) 

Name: Nasser Abbasi

1.  Find a set of vectors that span the hyperplane w - 2.3 x- 2 y- 3.4z==6

Solution outline: 

First find a normal to the hyperplane. 

Use the fact that given a general equation of a hyperplane when written in the form

                                 (Subscript[c, 1])(Subscript[x, 1])+(Subscript[c, 2])(Subscript[x, 2])+...+(Subscript[c, n])(Subscript[x, n])==d
            
Then a normal is just the vector {(Subscript[c, 1]),(Subscript[c, 2]),...,(Subscript[c, n])}.  For example, given the equation of a line, written in the above form, it is a x+b y=c Hence a normal to this line is the vector {a,b}.

Once the normal is found, then we can obtain basis for the null space of this normal. The nullspace of this normal gives us the basis vectors which span the space that this normal can not reach. 

But the space that this normal vector can not reach is the hyperplane itself. 

Hence all what we need to do is to use NullSpace on a Matrix which contains one row, which is the normal vector.
                        
step 1. Find the normal vector:

Remove["Global`*"];
normalVectorToHyperplane={1,-2.3,-2,-3.4}

Step 2. Find the NullSpace:

A=NullSpace[{normalVectorToHyperplane}]

2. Find the solution to the equations 
2x+3y+2z=10
5x-2y+7z=0
which has minimum length.

Solution outline:
Write the system as Ax=b and use  PseudoInverse to find the minimum solution

Remove["Global`*"];
Needs["LinearAlgebra`MatrixManipulation`"];
eq1=2x+3y+2z==10;
eq2=x+y-z==11;

Convert equations to matrix form

{A,b}=LinearEquationsToMatrices[{eq1,eq2},{x,y,z}];
MatrixForm[A]

MatrixForm[b]

And now use PseudoInverse to find the minimum solution

v=PseudoInverse[A].b//N

The above is the solution which gives the minimum length which is

Sqrt[v.v]

3. Calculate (∫_(-∞)^∞)((x^2)/(((1+(x^2)))^2))dx using residues.  Show that it is the real part of a complex integral . Specify the contour you are using by drawing a sketch. Explain your reasoning. Check your answer using Integrate.

Solution outline:

Find the poles of the function that we are integrating. (the integrand in the above). Then Find a laurent series expanded around this pole, and find the coeff. of the term (1/(z-pole)). This is the residue. The solution is then, by the residue theorem, is the sum of all the residues * 2πi. In this problem there is one pole, hence one residue is found.

Step 1. Start by finding the pole of the integrand

Remove["Global`*"];
f[z_]:=((z^2)/((((1+(z^2)))^2)))
poles=Solve[Denominator[f[z]]==0,z]

Step 2. I need to consider only poles that are inside the contour, which will be in the upper half plan. Hence there are pole are z=I of order 2. Find Laurent series around the pole, and pick up the coefficient of the term (1/(z-pole))

s=Series[f[z],{z,I,6}]

step3 . Now find the correct coeff.

residue=Coefficient[s,(1/(z-I))]

Step 4. Now use the residue theorem to find the value of the integral

result=Abs[2π I  residue]

Now I need to show it is the real part of a complex integral. To do this I consider an integral of f[z]=((z^2)/(((1-(z^2)))^2))and write the integral as a contour integral that extends over the real x axis from -r to +r and then goes anticlock wise in a circle of radius r from θ=0 to θ=π as shown in this diagram


Now since z=((re)^(iθ)) then dz=((rie)^(iθ))dθ and so the contour integral can be written as

integral=(∫_(-r)^r)((x^2)/(((1+(x^2)))^2))dx+(∫_0^π)((((r (e^(i θ))))^2)/(((1+(((r (e^(i θ))))^2)))^2))dθ

Now since we want to integrate from x=-∞ to x=∞, then we start by making r larger and larger.

Then we can take r to go to infinity. But notice in the second integral above, in the denominator, r is to the power of 4, while in the numerator r is to the power of 2. 

Hence the second integral above (the integral over the circular part) will vanish as r becomes larger and larger. 

Hence the integral becomes the real part only which is what we are asked to evaluate which is

integral=(∫_(-∞)^∞)((x^2)/(((1+(x^2)))^2))dx

Hence we showed that (∫_(-∞)^∞)((x^2)/(((1+(x^2)))^2))dx is the real part of contour complex integral (∫_(-r)^r)((z^2)/(((1+(z^2)))^2))dz where z is a complex variable z=(re^(iθ))

Now I check my answer using Integrate

Integrate[f[z],{z,-Infinity,Infinity}]

The above matches the result using residue method as shown above.

4.  Find all solutions to (z^7)+2(z^5)-3(z^2)+5z+1==0 Make a plot which shows the location of the roots in the complex plane.

Solution outline:
Since this is a polynomial, we can use Solve. Since this is a polynomial of degree 7, this by the fundamental theory of algebra, we will have 7 roots. Use the command ListPlot to make a plot of the roots in the complex plane.

Remove["Global`*"];
f[z_]:=(z^7)+2(z^5)-3(z^2)+5z+1
sol=z/.Solve[f[z]==0,z]//N

ListPlot[Transpose[{Re[sol],Im[sol]}],PlotStyle->PointSize[.03],AxesLabel->{"x","I y"}];

5.  Find the 4 smallest positive real  roots of x (e^(-x)) = Cot[x+1]

Solution outline:
Make a plot (try different ranges until you see one which contains at least 4 real roots). Then use FindRoot to locate exact roots.

Find make a plot to see the shape

eq=x (E^(-x))-Cot[x+1];
Plot[eq,{x,0,12}];

Now by looking at the above plot, I see the smallest 4 positive real roots around 0,4,6,10. So I can use FindRoot with an initial guess at these values. Since my starting guess values are all real, then FindRoot will only find real roots. I also supply a max/min brackets for the range of search by FindRoot for each call.

getRoots[f_,var_,at_]:=x/.FindRoot[f==0,{var,at,at-2,at+2}];

Hence the 4 smallest positive roots are

getRoots[eq,x,{0,4,8,12}]//TableForm

6.  A drunken sailor (or a molecule adsorbed on a surface, or a polymer chain) takes a step of unit length in either the +x,-x,+y, or -y directions every unit of time. Write a procedure which will draw a typical trajectory of nstep steps. Use a Do loop, AppendTo, Random[Integer] and ListPlot with the PlotJoined and AspectRatio-> Automatic option. Draw a trajectory of 2000 steps.

Solution outline:
Make an empty list. Generate random integer from 1 to 4. depending on the value make the appropriate step which will be an new entry, add it to the list. At the end use ListPlot to display it.

Declare data and do initialization

Remove["Global`*"];
Needs["Graphics`Colors`"]
nSteps=2000;
history={};
x=0 ;y=0;

Start the main loop.

Do[
m=Random[Integer,{1,4}];

Which[m==1,x=x+1,
m==2,x=x-1,
m==3,y=y+1,
m==4,y=y-1
];

AppendTo[history,{x,y}],

{n,nSteps}
];

ListPlot[history,
PlotJoined->True,
AspectRatio->Automatic,
PlotStyle->{Red,Thickness[.005]}
];

7. Two identical masses are attached to three springs between two fixed walls, as shown in the figure below. The spring constant of the middle spring is different from the outside springs. Neglect gravity.

Solution Outline:
Generate the lagrangian, then the equation of motion, from this we find the mass matrix and the stiffness matrix. Solve Det([k]-w^2 [M])=0 to find the eigenfrequencies.

First Update the above diagram as shown below to make it clear where displacement are

a) Construct the Lagrangian , and set up the equations of motion for this system

Set up symbols using notation for ease of reading.

Remove["Global`*"];
Needs["Utilities`Notation`"];
DynamicSymbolize[f_,k_]:=Symbolize[NotationBoxTag[SubscriptBox[f,k]]]
DynamicSymbolize["u","1"];
DynamicSymbolize["u","2"];
DynamicSymbolize["k","1"];
DynamicSymbolize["k","2"];

Now define the kinetic energy KE and the potential energy PE

KE=(1/2)m ((((Subscript[u, 1])'[t]))^2)+(1/2)m ((((Subscript[u, 2])'[t]))^2)

PE=(1/2)(Subscript[k, 1])((((Subscript[u, 1])[t]))^2)+(1/2)(Subscript[k, 2])((((Subscript[u, 1])[t]-(Subscript[u, 2])[t]))^2)+(1/2)(Subscript[k, 1])((((Subscript[u, 2])[t]))^2)

Find the lagrangian L

L=KE-PE//Simplify

Here I can use EulerEquations[] in ["Calculus`VariationalMethods`"] to get the equations, or do it by hand. This is easy enough so I just write it down from the definition (this way I also do not forget it)

eq1=D[D[L,(Subscript[u, 1])'[t]],t]-D[L,(Subscript[u, 1])[t]]==0//Simplify

eq2=D[D[L,(Subscript[u, 2])'[t]],t]-D[L,(Subscript[u, 2])[t]]==0//Simplify

(eqs={eq1,eq2})//TableForm

I can also display the above EQM is a more standard format as follows (even though this is not needed to proceed)

Needs["LinearAlgebra`MatrixManipulation`"]
vars={(Subscript[u, 1])''[t],(Subscript[u, 2])''[t]};
deps={(Subscript[u, 1])[t],(Subscript[u, 2])[t]};

{massMatrix,b}=LinearEquationsToMatrices[eqs,vars];
{stiffnessMatrix,doNotCare}=LinearEquationsToMatrices[b,deps];

Print[MatrixForm[massMatrix],MatrixForm[Transpose[List[vars]]],"+",MatrixForm[stiffnessMatrix],MatrixForm[Transpose[List[deps]]],"== 0"]

Multiply the LHS and RHS of the second equation above by minus sign to make the mass matrix all positive. Mathematica does not do it automatically. Notice there is a "-m" in the mass matrix above, so the code below just rewrite the above so that the mass matrix is all positive.

massMatrix[[2,2]]=-massMatrix[[2,2]];
stiffnessMatrix[[2,1]]=-stiffnessMatrix[[2,1]];
stiffnessMatrix[[2,2]]=-stiffnessMatrix[[2,2]];
Print[MatrixForm[massMatrix],MatrixForm[Transpose[List[vars]]],"+",MatrixForm[stiffnessMatrix],MatrixForm[Transpose[List[deps]]],"== 0"]

b) Calculate the normal mode frequencies for the vibrations of this system.

(charMatrix=stiffnessMatrix+λ massMatrix)//MatrixForm

eq=Det[charMatrix]==0

{(Subscript[λ, 1]),(Subscript[λ, 2])}=λ/.Solve[eq,λ]//Simplify

Finally find the harmonic frequencies. Take the positive square root of the above solutions

{(Subscript[ω, 1]),(Subscript[ω, 2])}=Map[Sqrt,%]

8. a) The unit sphere  and a cylinder defined by (((y-1))^2) +(z^2) =1 intersect. The region which is inside both the cylinder and the sphere is a peculiar shape. Use ContourPlot3D and RealTime3D to draw a picture of the region which can be rotated so that the region can be examined.

Solution outline:
Divide the volume that is common between the cylinder and the sphere into 2 parts. One is a standard volume shape (dome type shape) with a standard formula for its volume. The second volume, use volume integration to find.

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

cylinder=(((y-1))^2)+(z^2)-1;
sphere=(x^2)+(y^2)+(z^2)-1;

p1=ContourPlot3D[cylinder,{x,-1,1},{y,-1.999,1.999},{z,-1,1},PlotPoints->{4,5,5},ContourShading->False,DisplayFunction->Identity];

p2=ContourPlot3D[sphere,{x,-1,1},{y,-1,1},{z,-1,1},PlotPoints->{4,5,5},ContourShading->False,DisplayFunction->Identity];

<<RealTime3D`

Show[{p2,p1},DisplayFunction->$DisplayFunction];

<<Default3D`

b) What is the volume of this region?

IMPORTRANT:

Since I am not sure if my analytical solution is correct below, I just thought that I could have solved this problem much more easily numerically. 

I am not sure if we have to use analytical solution only. But I give this solution below, I get 40% as the covered volume which looks more closer to what I see in the plot. With the analytical solution I get 45%, and I think I over integrated there. So please accept this numerical solution below unless an analytical solution is needed in which case please accept my analytical solution below it.  Thank you.

Use MonteCarlo method.

v=NIntegrate[If[  (x^2)+(y^2)+(z^2)<=1 &&   (((y-1))^2)+(z^2)<1 ,1,0],{x,-1,1},{y,-1,1},{z,-1,1},Method->MonteCarlo,PrecisionGoal->3,MaxPoints->1000000]

vSphere=(4/3)π//N

Hence the percentage of the volume found to that of the sphere is

100(v/vSphere)

Analytical solution (In case the above is not acceptable)

First find the curve where the cylinder and the sphere meet

sol=Flatten[Solve[cylinder==sphere,{x,y,z}]]

Now find y for x=0

yp=y/.sol

y1=yp/.x->0

Now the above says y=.5 where the cylinder meets the sphere. Find the z coordinates where this happens

{z1,z2}=z/.Solve[sphere==0/.{y->y1,x->0},z]
N[%]

Hence y=0.5 is where the intersection started. That is, if I draw a flat surface at y=.5 that is perpendicular to the y-axis, it will meet the cylinder just as it entered the sphere, like this

Now imagine a flat surface y=.5 that cuts the sphere, then I have 2 volumes to calculate, one above this surface  and one below it as shown. The top surface has a known formula for it which I show below. The bottom surface I calculate using integration as will be done below

Now I need to find the volume of the segment below the cut-off plane. 
A section of the area is given by A=(∫_(y=0)^(y=.5))(Sqrt[1-(((y-1))^2)])dy  Hence the volume is (∫_(x=-1)^(x=1))Adx , next multiply this by 2 since it goes around all the way

v2=(∫_(-1)^1)(∫_0^(.5))(Sqrt[1-(((y-1))^2)])\[DifferentialD]y\[DifferentialD]x

v2=2*Chop[v2]

Now add to the above the volume of the top part as described above, which is given by (1/2)((πR)^2)h where R=((Sqrt[3])/2)as the solution for z from above and h=0.5 found above as well

v1=(1/2)π((Sqrt[3])/2)((1/2))

Hence total volume is

v=v1+v2//N

Note the volume of the unit sphere is

vSphere=(4/3)π//N

Hence the percentage of the volume found to that of the sphere is

100(v/vSphere)

The above is just a little less than half the sphere volume.

9. Express the electric field of a point charge at the origin using cartesian , cylindrical and spherical coordinates and basis vectors.  The JacobianMatrix, Scalefactors, and CoordinatesToCartesian functions from the VectorAnalysis package may be useful.

Solution outline:
The electric field is the electric field force per unit charge

                                       (Overscript[E, -])=((Overscript[F, -])/q)
                                       
where 
                                         (Overscript[F, -])=(Q/(4πϵ))(1/(r^2))(Overscript[u, -]) 
                                         
where (Overscript[u, -]) is a unit vector along the direction of the electric force vector (Overscript[F, -])

Since we are using a unit charge, and use a unit test charge as well, then

                                         (Overscript[E, -])=(1/(4πϵ))(1/(r^2))(Overscript[u, -]) 
                                         
But (Overscript[u, -])=((Overscript[F, -])/(|(Overscript[F, -])|)) and r is the length of (Overscript[F, -]) hence

                                         (Overscript[E, -])=(1/(4πϵ))(1/(|(Overscript[F, -])((|)^2)))((Overscript[F, -])/(|(Overscript[F, -])|)) 
                                         
So the problem is 2 parts: Express a vector (Overscript[F, -]) in different coordinates systems, and express the length of this vector in different coordinates systems. 

step1. Start by creating some symbols

Remove["Global`*"];
Needs["Calculus`VectorAnalysis`"]
Needs["Utilities`Notation`"];
DynamicSymbolize[f_,k_]:=Symbolize[NotationBoxTag[SubscriptBox[f,k]]]
DynamicSymbolize["e","1"];
DynamicSymbolize["e","2"];
DynamicSymbolize["e","1"];

DynamicSymbolize["e","x"];
DynamicSymbolize["e","y"];
DynamicSymbolize["e","z"];

DynamicSymbolize["e","r"];
DynamicSymbolize["e","θ"];
DynamicSymbolize["e","ϕ"];

Define a function which takes general vector of the form (Subscript[a, x]) (Subscript[e, x])+(Subscript[a, y]) (Subscript[e, y])+(Subscript[a, c]) (Subscript[e, c]) where the a's are coordinates and the e's are the base vectors associates with these coordinates. For example in Cartesians, the vector is (Subscript[a, x]) (Overscript[i, -])+(Subscript[a, y])(Overscript[j, -])+(Subscript[a, z])(Overscript[k, -]) then this function will return its length

length[{a_ ,e1_,b_ ,e2_,c_ ,e3_}]:=Module[{v},
v={a e1.e1,b e2.e2,c e3.e3};
Sqrt[v.v]
]

In Cartesian

ei={1,0,0};
ej={0,1,0};
ek={0,0,1};

electricVec=x  ei+y ej+z ek;
len=length[{x ,ei,y ,ej,z ,ek}];

unitVec=(electricVec/len);

electricFieldVec=(1/(4 π ϵ))((unitVec/(len^2)))//Simplify

ContourPlot[(Sqrt[(((electricFieldVec[[1]])^2)+((electricFieldVec[[2]])^2))/.{z->0.00001, ϵ->1}]),{x,-100,100},{y,-100,100}];

Now in  cylindrical

First find the unit vectors (Subscript[e, r]),(Subscript[e, θ]),(Subscript[e, z]) in cylindrical. Find define the relationship between the cartesian and the cylindrical coordinates

xCoord[r_,θ_,z_]:=r Cos[θ];
yCoord[r_,θ_,z_]:=r Sin[θ];
zCoord[r_,θ_,z_]:=z;

Find (Subscript[e, r])

(Subscript[e, r])=(∂_r){xCoord[r,θ,z],yCoord[r,θ,z],zCoord[r,θ,z]};
(Subscript[e, r])=((Subscript[e, r])/(Sqrt[(Subscript[e, r]).(Subscript[e, r])]))//Simplify

Find (Subscript[e, θ])

(Subscript[e, θ])=(∂_θ){xCoord[r,θ,z],yCoord[r,θ,z],zCoord[r,θ,z]};
(Subscript[e, θ])=((Subscript[e, θ])/(Sqrt[(Subscript[e, θ]).(Subscript[e, θ])]));
(Subscript[e, θ])=Simplify[%,Assumptions->{r≥0}]

Find (Subscript[e, z])

(Subscript[e, z])=(∂_z){xCoord[r,θ,z],yCoord[r,θ,z],zCoord[r,θ,z]};
(Subscript[e, z])=((Subscript[e, z])/(Sqrt[(Subscript[e, z]).(Subscript[e, z])]))//Simplify

Hence the electric field vector in cylindrical is

electricVec=xCoord[r,θ,z]  (Subscript[e, r])+yCoord[r,θ,z] (Subscript[e, θ])+zCoord[r,θ,z] (Subscript[e, z])

len=length[{xCoord[r,θ,z], (Subscript[e, r]),yCoord[r,θ,z] ,(Subscript[e, θ]),zCoord[r,θ,z] ,(Subscript[e, z])}]//Simplify

unitVec=(electricVec/len)

electricFieldVec=(1/(4 π ϵ))(unitVec/(len^2))//Simplify

Now do the same for the spherical coordinates. First define the relationship between the cartesian and the spherical coordinates

Now find the unit vectors in spherical coordinates

xCoord[r_,θ_,ϕ_]:=r Cos[ϕ]Sin[θ];
yCoord[r_,θ_,ϕ_]:=r Sin[θ] Sin[ϕ];
zCoord[r_,θ_,ϕ_]:=r Cos[θ];

Find (Subscript[e, r])

(Subscript[e, r])=(∂_r){xCoord[r,θ,ϕ],yCoord[r,θ,ϕ],zCoord[r,θ,ϕ]};
(Subscript[e, r])=((Subscript[e, r])/(Sqrt[(Subscript[e, r]).(Subscript[e, r])]))//Simplify

Find (Subscript[e, θ])

(Subscript[e, θ])=(∂_θ){xCoord[r,θ,ϕ],yCoord[r,θ,ϕ],zCoord[r,θ,ϕ]};
(Subscript[e, θ])=((Subscript[e, θ])/(Sqrt[(Subscript[e, θ]).(Subscript[e, θ])]))//Simplify;
(Subscript[e, θ])=Simplify[(Subscript[e, θ]),Assumptions->{r≥0}]

Find (Subscript[e, ϕ])

(Subscript[e, ϕ])=(∂_ϕ){xCoord[r,θ,ϕ],yCoord[r,θ,ϕ],zCoord[r,θ,ϕ]};
(Subscript[e, ϕ])=((Subscript[e, ϕ])/(Sqrt[(Subscript[e, ϕ]).(Subscript[e, ϕ])]))//Simplify;
(Subscript[e, ϕ])=Simplify[(Subscript[e, ϕ]),Assumptions->{r≥0}];
(Subscript[e, ϕ])=(Subscript[e, ϕ])/.Sqrt[any_^2]->any (*Strange mathematica*)

Hence the electric field vector in spherical is

electricVec=xCoord[r,θ,ϕ]  (Subscript[e, r])+yCoord[r,θ,ϕ] (Subscript[e, θ])+zCoord[r,θ,ϕ] (Subscript[e, ϕ])//Simplify

len=length[{xCoord[r,θ,ϕ], (Subscript[e, r]),yCoord[r,θ,ϕ] ,(Subscript[e, θ]),zCoord[r,θ,ϕ] ,(Subscript[e, ϕ])}]//Simplify

unitVec=(electricVec/len)//Simplify
(*unitVec=Simplify[unitVec,Assumptions->r≥0]*)

electricFieldVec=(1/(4 π ϵ))(unitVec/(len^2));
electricFieldVec=electricFieldVec/.Sqrt[any_^2]->any;
electricFieldVec=electricFieldVec/. (1/(Sqrt[any_^2]))->(1/any)

10. Data for the counts recorded in a radioactive decay experiment, in the form {time, counts/min} can be found  in a file called decay.dat, which is in the localdata directory defined by

Remove["Global`*"]
localdata=$TopDirectory<>"\\AddOns\\PTHBdata"

Import the data into a list.   The counter registers a background count rate even if there are no decays in the sample. The isotope has a decay constant of -0.3 , so a model for the count rate at time t is A +B (e^(-0.3 t)), where A is the background count rate.  Find A and B using Fit. Plot the data and the best fit function on the same plot using two different colors.

Needs["Graphics`Colors`"]

data=Import[localdata<>"\\decay.dat"];
{nRow,nCol}=Dimensions[data]

p1=ListPlot[data,PlotStyle->Yellow,DisplayFunction->Identity];

Fit the data

f=Fit[data,{1,Exp[-0.3  t]},t]

Find A and B

{A,B}=f/.a_+b_ Exp[any_]->{a,b}

p2=Plot[f,{t,data[[1,1]],data[[nRow,1]]},PlotStyle->Red,DisplayFunction->Identity];

Show[{p1,p2},AxesLabel->{"time","counts/min"},DisplayFunction->$DisplayFunction];

11. Consider the eigenvalue boundary value problem

(((d^2)y)/(d (x^2))) +λ (e^x) y=0         y[0]=0,  y[1]=0

Find the first 3 eigenvalues and plot the  corresponding solutions.

Solution outline: Since this is both an eigenvalue problem and a boundary value problem, DSolve nor NDSolve can solve it as is. Start by converting the problem to an initial value problem. Use y'[0]==1 since it does not matter which initial value for the "speed" we select, we might as well select y'[0]=1.

Next, try a range of eigenvalues, for each use NDSolve to find the solution and evaluate the solution at y[1].  Plot y[1] as a function of λ. Now use FindRoot to find where y[1] is zero on this plot. This gives us the values of λ which when used in the original ODE will result in the correct boundary condition met.

Remove["Global`*"];

ode=y''[x]+λ (E^x) y[x];
yp0=1;

yAtOne[yp0_?NumericQ,eigenValue_?NumericQ]:=y[1]/.NDSolve[{ode==0/.λ->eigenValue,y[0]==0,y'[0]==yp0},y,{x,0,2}]//First

Plot[yAtOne[yp0,λ],{λ,0,150},FrameLabel->{"λ","y[1]","Showing dependency of solution at 1 on eigenvalue",None},Frame->True];

Now by inspection and using FindRoot we find the 4 smallest eigenvalues

FindRoot[yAtOne[yp0,λ],{λ,0,10}]

FindRoot[yAtOne[yp0,λ],{λ,15,20}]

FindRoot[yAtOne[yp0,λ],{λ,45,55}]

FindRoot[yAtOne[yp0,λ],{λ,90,100}]

12. For what values of λ does the DE  x y''[x] +(1-x) y'[x] +λ y[x]=0  have polynomial solutions?

Remove["Global`*"];
ode=x y''[x]+(1-x)y'[x]+λ y[x]

Solution outline:
Assume the solution is a power series. Substitute this solution in the ODE. We find that we get one possible solution with one free parameter, say a[0].

Now look at the series in x which goes with this free parameter a[0].

This series will have λ in it as well.

We want to look at the conditions on λ which makes this series TERMINATE. 

This indicates a convergent solution. Using this condition, we can see for what values λ must have to terminate the a[0] series. Hence these will be the λ values which gives a polynomial solution in x.

step 1: Assume a power series solution

maxPower=4;
ys=(∑_(n=0)^maxPower)a[n] (x^n)

Plug in the ODE

seriesODE=ode/.y->Function[x,Evaluate[ys]]

Now since the ODE=0, then each power of x must be zero as well for this to happen. Use this to generate a set of equations

eqs=Table[Coefficient[seriesODE,x,n]==0,{n,0,maxPower}];
TableForm[eqs]

Drop the last equation since incomplete due to the side-effect of using series

eqs=Drop[eqs,-1];
TableForm[eqs]

We see from the above that if we solve for a[1] in terms of a[0], then we can solve for a[2] in terms of a[1], and then solve for a[3] in terms of a[2], etc... Hence we use a[0] as a free parameter.

sol1=Solve[eqs[[1]],a[1]]

sol2=Solve[eqs[[2]],a[2]]/.sol1//First//Simplify

sol3=Solve[eqs[[3]],a[3]]/.sol2//First//Simplify

sol4=Solve[eqs[[4]],a[4]]/.sol3//First//Simplify

replacments=Flatten[Union[sol1,sol2,sol3,sol4]]

Now plug in all the above into the power series trial  solution

ys=ys/.replacments

Collect[ys,a[0]]

To make it a little easier, let me just remove the a[0] term out of the way for now

main=%/.any_ a[0]->any

So now, the above is my main object to examine. I need the above series in x to TERMINATE. Here I am only looking at 4 terms because that is what I tried, but more terms can be added. But these are enough to see the pattern.The point is that we want to find the conditions which causes this series to terminate. Use collect to see it better

Collect[main,{x,λ}]

Factor[-(λ/4)+((λ^2)/4)]

Factor[-(λ/18)+((λ^2)/12)-((λ^3)/36)]

Factor[-(λ/96)+((11 (λ^2))/576)-((λ^3)/96)+((λ^4)/576)]

Ok, the pattern is clear now.

From the above we see the following. If λ=1, then that will remove the '(x^2)' term due to the (-1+λ) factor , and it will also remove all the terms in higher power of x after that, since we see that there is a factor (-1+λ) in each term after that. Hence λ=1 is a good eigenvalue in that it gives us a terminating series in x. 

For λ=2, see the same pattern results. It will make the (x^3) term vanish due it having a (-2+λ) in it. And all higher terms will all vanish as well since each of them has a (-2+λ) factor.

Hence we see that only for λ an integer will the series terminate, giving a polynomial solution. We continue this way.

So for λ=0, the solution is (Subscript[a, 0])(1)
So for λ=1, the solution is (Subscript[a, 0])(1-x)
For λ=2, the solution is (Subscript[a, 0])(1-2x+((x^2)/2))
For λ=3, the solution is

(Subscript[a, 0])(1-x λ+(x^2) (-(λ/4)+((λ^2)/4))+(x^3) (-(λ/18)+((λ^2)/12)-((λ^3)/36)))/.λ->3

Etc.....

Let (Subscript[a, 0])=1 (since free parameter), I get


So for λ=0, the solution is 1
So for λ=1, the solution is (1-x)
For λ=2, the solution is (1-2x+((x^2)/2))
For λ=3, the solution is (1-3 x+((3 (x^2))/2)-((x^3)/6))

This completes the solution. Hence only for 

                                         λ=0,1,2,.... 
                                         
will we get a polynomial solution in x. (for λ=0 it is a constant solution=1).

Note that we rejected non-integer λ values since these will not terminate the series. Also negative λ will not terminate the series.

Btw, it is possible to first see what type of power series solution for this ODE one can expect. The above ODE can be written in the form

                    y''(x)+((1-x)/x)y'(x)+(λ/x)y(x)=0

Now look at the functions that go with each coef. above. We see the coefficients functions are:   ((1-x)/x),(λ/x) and these are clearly not analytic at x=0, they blow up. Hence a series solution around x=0 can not be ordinary. It can be only regular singular or irregular singular.

If we multiply each coeff. function by (x-x0) or by x since we are expanding around x0=0, then we can remove all of the poles. 

Hence all of the coeff. functions become analytic. 

Hence the series solution is considered a regular singular. 

By Frobenius theorem, for a regular singular point, it guarantees the existence of a least one (or possible two) solutions of the form  (Subscript[y, 1])(x)=(x^(Subscript[b, 1]))(∑_(n=0)^∞)(Subscript[a, n])(x^n)where (Subscript[b, 1]) not necessarily integer. This is the solution found above.

In addition it says, if there is only one solution of this form, the second solution must be of the form (Subscript[y, 2])(x)=(Subscript[y, 1])Log[x]+(x^(Subscript[b, 2]))(∑_(n=0)^∞)(Subscript[c, n])(x^n)where (Subscript[b, 1]) and (Subscript[b, 2]) are roots to a quadratic equation called the indicial equation. I do not need to consider this second case, since the first solution (Subscript[y, 1]) gives us the conditions on λ we are asked to find.

13.  The force on a particle with charge q and mass m is given by the familiar formula F=q(E+v×B). Consider the case of a particle of mass m and charge q starting from rest at the origin in a region of space with a constant B field in the z direction and a constant E field in the x direction.
a) Use dimensional analysis to determine the characteristic length, time and energy in this problem.

Remove["Global`*"]
<<"MMHTools`DimTools`"
Needs["Miscellaneous`PhysicalConstants`"]
Needs["Miscellaneous`Units`"]

chargeUnits=Coulomb; (*q*)

electricFieldUnits=(Newton/Coulomb); (*E*)

velocityUnits=(Meter/Second); (*v*)

magneticFieldUnits=Tesla; (*B*)

massUnits = Kilogram; (*m*)

gravityUnits = (Meter/(Second^2)); (*g*)

Now start doing dimensional analysis to find the length scale of the problem (distance from the origin at some future time t). i.e. we are asking the question: 

What does the distance away from the origin that the particle will be  some time t depends on? 

This distance must depend on how strong the magnetic field B is, and on how strong E is, and on the mass of the particle, and on the amount of electric charge on the particle, and on the velocity of the particle. Hence we write

condp={q chargeUnits,ef electricFieldUnits,b magneticFieldUnits,m massUnits}//ToSymbolsUnits//ReduceUnits

What possible combination can give us a length variable out of the above parameters? use dimanal[]

dimanal[condp,{Kilogram,Meter,Second,Coulomb},Meter]

Verify the above answer

electricFieldUnits*massUnits/((magneticFieldUnits^2) chargeUnits)//ReduceUnits

Hence the length scale is given by 
                                           ((Subscript[L, r]))((E m)/((B^2) q))
where ((Subscript[L, r])) is a constant (non-dimensional factor)  i.e 

                                         length ~ ((E m)/((B^2) q))
                                         
To find what is the length is proportional to, I need to use some physical values. Assume the particle is an electron, we know the mass m and the charge q. For E and B, 

Now find a characteristic time

dimanal[condp,{Kilogram,Meter,Second,Coulomb},Second]

Find characteristic energy

dimanal[condp,{Kilogram,Meter,Second,Coulomb},Joule//ReduceUnits]

b)  Find the analytic solution to the equations of motion. What is the period of the motion? What is the maximum and minimum values of the kinetic energy of the particle as a function of time? Comment on how these quantities are related to the characteristic values found in part a.

(Subscript[e, k])={0,0,1};
(Subscript[e, i])={1,0,0};
(Subscript[e, j])={0,1,0};
eletricFieldVec={ex (Subscript[e, i])+ey (Subscript[e, j])+ez (Subscript[e, k])}//First
magneticFieldVec={bx (Subscript[e, i])+by (Subscript[e, j])+bz (Subscript[e, k])}//First
velocityVec={rx'[t] (Subscript[e, i])+ry'[t] (Subscript[e, j])+rz'[t] (Subscript[e, k])}//First
q(eletricFieldVec+Cross[velocityVec,magneticFieldVec])

Now since we are told that E is constant and along the x direction only, we set the electric field vector coordinates that way. Also since we are told that the magnetic field B is constant and along the z direction, we set its coordinates to point that way as well.

forceOnParticleVec=%/.{ey->0,ez->0,bx->0,by->0}

Declare a position vector

positionVec={rx[t] (Subscript[e, i])+ry[t] (Subscript[e, j])+rz[t] (Subscript[e, k])}//First

Declare the ODE using F=ma

ode=forceOnParticleVec==m Map[D[#,{t,2}]&,positionVec]

{lhs,rhs}=ode/.lhs_==rhs_->{lhs,rhs}

setOfODE=Thread[Equal[lhs,rhs]]

initialCoditions={rx[0]==0,rx'[0]==0,ry[0]==0,ry'[0]==0,rz[0]==0,rz'[0]==0};
sol=Flatten[DSolve[Join[setOfODE,initialCoditions],{rx[t],ry[t],rz[t]},t]]

Find the period.
From above we see that the x component of the position contains a frequency ω = bz q/m. Hence the period is ((2π m)/(B q))

Verify the units

massUnits/(magneticFieldUnits*chargeUnits)//ReduceUnits

For the motion along the y axis, the period is also given by the same value above

Find K.E.

vx=D[rx[t]/.sol,t]

vy=D[ry[t]/.sol,t]

v=Sqrt[(vx^2)+(vy^2)]

KE=(1/2)m (v^2)//Simplify

D[KE,t]//Simplify

The above is zero when Sin[((bz q t)/m)]=0, which means ((bz q t)/m)=nπ for n=0,1,2,3,...
n=0 means t=0, then K.E. is minimum because then KE=0 (initial condition, particle is still at rest)
look at the next value. When n=1 then t=((m π)/(B q))then KE is

KE/.t->((m π)/(bz q))

When n=2 then t=((2 m π)/(B q)) find what is KE now

KE/.t->((2 m π)/(bz q))

Hence we see that for n=1,3,5,7,....KE is max, and when n=0,2,4,6,... KE is min.

hence MAX KE is ((2 (E^2) m)/(B^2))which occurs at times t= n ((m π)/(bz q)) for ODD n

and MIN KE is 0 which occurs at times t=n ((m π)/(bz q)) for EVEN n

Now  I  Comment on how these quantities are related to the characteristic values found in part a. 
From part (a) we found that the energy is proportional to (((E^2)m)/(B^2)) and here we just found the maximum KE is ((2 (E^2) m)/(B^2)), Hence this agree with the finding for the KE energy part with a proportionality constant of 2.

In part(a) we found the time scale to be proportional to (m/(B q))and here we found the period to be ((2π m)/(B q))which agree with part (a), and the proportionality constant is 2 π

In part(a) we found the length scale to be ((E m)/((B^2) q)), and here we found the period to be ((2π m)/(B q)), and since the velocity is as shown above, I can find the distance travelled during one period, and this is taken as the length scale. Let me find this distance now:

period=((2 π m)/(bz q));
distanceTravelledInOnePeriod=v*period//Simplify

The stuff under the square root above can be written as K (Sqrt[(((ex^2) )/(bz^2))]) for some constant K, hence the distance above can be written as

((2 (Sqrt[2]) m π  K (ex/bz))/(bz q)) which is 

                              (Subscript[k, 1]) ((E m)/((B^2) q))

Where (Subscript[k, 1]) is a new constant which consumed (k 2 (Sqrt[2])π)

Comparing the above expression for the distance, we see it agrees which part (a), where the proportionality constant is (Subscript[k, 1])

c) For the special case of a proton in an E field of 1000 Volts/meter and a B field of 0.1 Tesla, draw a graph of the trajectory of the particle.

First check what the period is now that we have some numbers to work with

ProtonMass

ElectronCharge

Hence period is (in seconds):

period=((2 π m)/(b q))/.{m->1.67262158*^-27,q->1.602176462*^-19,b->0.1}

Now generate plot points to make a plot of the trajectory of the particle

x=rx[t]/.sol

y=ry[t]/.sol

pts=Table[{x,y}/.{ex->1000,bz->0.1,m->1.67262158*^-27,q->1.602176462*^-19},{t,1,1+100*period,period}];

ListPlot[pts,PlotRange->All,PlotJoined->True,AxesLabel->{"x","y"},PlotRange->All];

This plot is tricky, I make another plot with longer time intervals, it shows the periods better

pts=Table[{x,y}/.{ex->1000,bz->0.1,m->1.67262158*^-27,q->1.602176462*^-19},{t,0,10,0.1}];
ListPlot[pts,PlotRange->All,PlotJoined->True,AxesLabel->{"x","y"},PlotRange->All];

14. (grads) Let y[x] be a curve with endpoints symmetrically positioned about the y axis.  Consider the surface area of the surface generated by rotating y[x] about the x axis to make a surface of revolution.  Find the y[x] which yields the surface of minimum area. This is the shape adopted by a soap film supported by two circular hoops; the soap film  minimizes its surface area because there is an energetic cost per unit of area of film.

a) Write an expression for the surface area in terms of an integral.

Imagine a thin strip all the way around. This strip will have differential width of




ds = (Sqrt[(dx^2)+(dy^2)])

Hence the strip will have an differential area of 

dA= circumference * ds = 2π*y[x] * (Sqrt[(dx^2)+(dy^2)])

Now assume x goes from -1 to +1. (or we can write -L to +L, it does not matter). Now, due to symmetry, we integrate from x=0 to 1 and multiply by 2.

Hence the area A= 2(∫_0^(x=1))2πy[x]*(Sqrt[(dx^2)+(dy^2)])

But (Sqrt[(dx^2)+(dy^2)])=(Sqrt[((dx^2)/(dy^2))+1]) dy
OR
(Sqrt[(dx^2)+(dy^2)])=(Sqrt[((dy^2)/(dx^2))+1]) dx

I will pick the First form above.

Hence substitute the above into the area integral A we get

                                        A= 4π(∫_(Subscript[y, 0])^(Subscript[y, 1]))y[x](Sqrt[((dx^2)/(dy^2))+1])dy

b) Find the Euler equation for the extremal curves.

Remove["Global`*"];
A=4 π (∫_(Subscript[y, 0])^(Subscript[y, 1]))y(Sqrt[1+(((x'[y]))^2)])\[DifferentialD]y

Minimizing A leads to the Euler equation (which comes from writing the function we want to minimize as F(y,x,x'[y]), and it leads to the euler equation which is 

                                                (d/dy)(dF/(dx'))-(dF/dx)=0

F=y(Sqrt[1+(((x'[y]))^2)])

eulerEquation = Hold[D[D[F,x'[y]],y]-D[F,x[y]]==0]

Look at the last term above

(∂_(x[y]))F

Hence eulerEquation becomes now just the first term, which is

Hold[D[D[F,x'[y]],y]==0]

Now since the above derivative is ZERO, then this means (∂_(x'[y]))F must be a CONSTANT. Say k, so now we have the following

eulerEquation=(∂_(x'[y]))F==k

c)Solve the resulting DE and look for solutions which are symmetric about the y axis. Show that the solutions are of the form A Cosh[b x] .

Square both sides above, we get 

(((y x'))^2)=(k^2)(1+((x')^2))

(((y^2) (x'))^2)=(k^2)+(k^2)(( (x'))^2)

(((x'))^2)((y^2)-(k^2))=(k^2)

Hence 

x'[y]=(k/(Sqrt[((y^2)-(k^2))]))

sol=Flatten[DSolve[x'[y]==(k/(Sqrt[(y^2)-(k^2)])),x[y],y]//Simplify]

I need to convert the above to Cosh. I do not know how to make Mathematica do it, so searching the Tables reference Handbooks I found the following trig. Relationship

                       ArcCosh(Z) = Log( Z + (Sqrt[(Z^2)-1]))
                       
 Hence in the solution I have above, I can rewrite it as follows to convert things to Cosh as needed by this problem

x=x[y]/.sol

Rewrite as follows

x=C[1]+k Log[y+k(Sqrt[-1+((y^2)/(k^2))])]

or

x=C[1]+k Log[k((y/k)+(Sqrt[-1+((y^2)/(k^2))]))]

But Log[A*B] = some constant * Log[B] so above becomes

x=C[1]+k Log[((y/k)+(Sqrt[-1+((y^2)/(k^2))]))]

x=C[1]+k ArcCosh[(y/k)]

x=C[1]+k ArcCosh[(y/k)]

Hence, let C[1]=K1, we write the above as

x-K1=k ArcCosh[(y/k)]

Hence

((x-K1)/k)= ArcCosh[(y/k)]

or

(y/k)=Cosh[((x-K1)/k)]

hence

                                               y[x] = K Cosh[((x-K1)/K)]
                                               
 Now Plot it, assume some constants

Remove["Global`*"];
y=k Cosh[((x-k1)/k)];

by trial and error, for {k->1,k1->0} I get the symmeterial picture I needed. For different constants, the plots are shifted to left or right.

y=y/.{k->1,k1->0}

Plot[{y,-y},{x,-3,3},PlotRange->{{-2,2},{-2,2}},AxesOrigin->{0,0}];

15.(grads) Quantum tunneling is one of the characteristic non-intuitive  phenomena of quantum mechanics. Consider a beam of particles of unit density traveling to the right, represented by (e^(i k x)) , as shown in the figure below. The beam encounters a region of space where the potential energy is V0. For a classical particle, if the incident kinetic energy Ein=(((ℏ^2)(k^2))/(2 m))<V0 the particle would reflect off the barrier, while if (((ℏ^2)(k^2))/(2 m))>V0 , the particle would pass through. For a quantum particle, some fraction of the beam is reflected, as represented by the R (e^(-i k x))  term, and some fraction is transmitted, as represented by the T (e^(i k x)) term.

I wish I have time to attempt this problem as it looks interesting but I have not yet started on HW9 which is due in 2 days, I better start on that now.

The Schrodinger equation that governs this problem inside the barrier is

Schin=((-(ℏ^2))/(2 m))ψ''[x]+V0 ψ[x]==En ψ[x]

and outside:

Schout=((-(ℏ^2))/(2 m))ψ''[x]==En ψ[x]

T is the transmission amplitude. The object of primary interest in this problem is the behavior of the transmission probability given by (SuperStar[T])T=((Abs[T])^2)  as a function of the initial energy Ein=(((ℏ^2)(k^2))/(2 m)). Graphs of the transmission probability as a function of energy are shown in most quantum books, but the algebra required to construct the plot is not trivial.

a) Using k as the primary variable as suggested in the Figure above, write down solutions in each of the three regions. The central region will involve 2 unknown parameters. Invoke the requirement that the solutions must be continuous and have continuous first derivative at the boundaries of the regions. Use this to solve for all of the unknown parameters.

b) What are the natural units of length and energy for this problem? Also identify any dimensionless parameters.

c) Express T as found in part a) in terms of Ein rather than k.

d) Using the characteristic energy found in b) , express T in terms of a reduced dimensionless energy Er

e) T(Er)  still contains a and V0. Eliminate these using the dimensionless parameter found in b) . When this is done, you should have the transmission amplitude as a function of strictly dimensionless quantities. Note that this implies that the transmission amplitude is not a function of a and V0 separately, but only a special dimensionless combination of them.

d) Make a plot of Abs[T]^2 as a function of the reduced energy for values of the dimensionless parameter=4, 1, 1/4 .

e) extra credit. Find a compact expression for (SuperStar[T])T . Simply computing this using ComplexExpand and Conjugate yields a rather complicated expression containing  Arg or ArcTan which are necessary to specify the signs of various parameters. Judicious rules to deal with these terms will yield a pleasant expression. Compare plots of this expression with those you made d).