Physics 229& 100 Homework #1
Name: Nasser Abbasi

(converted to 6.0)

1. Practice in transcribing expressions into Mathematica syntax
a) Find the determinant of this matrix: ({{ }, {}})

a = {{1, 2, 3, 4}, {1, -1, 0, 1}, {-1, 2, 6, 7}, {0, 2, 4, 1}}; 
Print["Determinant of the matrix is ", Det[a]]

b) If y=(∫_( 0)^( (Sqrt[x])))sin((t^2))dt  , what is ((d y)/(d x)) ?

Clear[y,x]
y[x_]:=Integrate[Sin[t^2],{t,0, (Sqrt[x])}];
Print["(dy/dx)=",D[y[x],x]];

c) In your E&M book, you find the following equation which describes an LCR circuit : L(((d^2)I)/(d (t^2)))+R ((d I)/(d t))+(I/C)=(Subscript[V, 0])cos ω t   . Use DSolve to find the general solution to the equation.

Clear[r,L0 ,i,t,v0,w,c];
eq=L0 i''[t]+r i'[t]+((i[t])/c)==v0 Cos[w t];
sol=DSolve[eq,i[t],t];
Print["I(t)=",i[t]/.sol[[1]]/.{C[1]->k1,C[2]->k2}];

d) Find the first 3 terms in the Taylor's series expansion of     about x=0.

Clear[f,x,a]
f[x_]:=(Sqrt[a+Log[x+1]]);
Normal[Series[f[x],{x,0,2}]]

2

The electromagnetic energy density inside a cavity at temperature T has a distribution of frequencies f. The distribution is given by the famous Planck formula  u(f)=((8 π (f^2))/(c^3))((hf/((e^(hf/kT))-1))), where h is Planck's constant, k is Boltzmann's constant, and c is the velocity of light.

a) Make a plot of the energy density as a function of f for T=100K and 300K. Be sure to include the physically relevant range of f. Make the curves different colors. Make a substitution rule to provide numerical values for the constants.

Remove["Global`*"];
c=QuantityMagnitude@UnitConvert[Quantity["SpeedOfLight"]];
h=QuantityMagnitude@UnitConvert[Quantity["PlanckConstant"]];
k=QuantityMagnitude@UnitConvert[Quantity["BoltzmannConstant"]];
minFrequency=0.001;
maxFrequency=3*10^14;
u[f_,T0_]:=(((8 Pi f^2) (h f))/(c^3(Exp[h f/(k T0)]-1)));

Plot[{u[f,100],u[f,300]},{f,minFrequency,10^14},PlotRange->All,FrameLabel->{"frequency in Hz","Energy Density (J per unit volume)"},PlotStyle->{Red,Blue},Frame->True,PlotLabel->"Planck Formula for different tempreture"]

b) Find the total energy in the cavity, which is given by the integral of the energy density over frequency. Leave the result in symbolic form( i.e. in terms of h, k, c, etc.) . You might want to include an Assumptions->Re[ ]>0 to avoid If statements in the answer. See examples of this in Integrate.

Clear[h,k,c]
u[f_]:= 8 Pi f^2/c^3( h f/(Exp[h f /(k T0)]-1))
sol=Integrate[u[f],{f,0,Infinity},Assumptions->Re[h/(k T0)]>0];
Print["Total energy is ",sol , " Joules"]

c) The power flux in the radiation is the total energy times the rate at which it moves (which is the speed of light, c) times a geometrical factor (1/4). Evaluate this and find the numerical coefficient of (T^4) ; this is the Stefan-Boltzmann constant of radiative transfer.

Clear[f,T0];
sol=Integrate[u[f,T0],{f,0,Infinity},Assumptions->Re[(h/(k T0))]>0];
sol=(1/4)*sol*c;

Print["coefficient of (T^4)=",coeff=Coefficient[Collect[sol,T0],T0^4]];
c=QuantityMagnitude@UnitConvert[Quantity["SpeedOfLight"]];
h=QuantityMagnitude@UnitConvert[Quantity["PlanckConstant"]];
k=QuantityMagnitude@UnitConvert[Quantity["BoltzmannConstant"]];
Print["Its numerical value is =",coeff ," J (K^(-4)) (m^(-2)) (s^(-1))"];

3

Conic sections can be described by the equation (1/r)=1+e Cos[θ], where e is the eccentricity, and r and θ are polar coordinates, with x=r Cos[θ] and y=r Sin[θ]. Use ParametricPlot to draw the curves for e=0.8, e=1, and e=1.5

Remove["Global`*"];
eq=(1/r)==1+e Cos[θ];
r=(1/(1+e Cos[θ]));
x=r Cos[θ];
y=r Sin[θ];
plotit[v_]:=ParametricPlot[Evaluate[{x/.e->v,y/.e->v}],{θ,0,2 π},PlotLabel->{"PLot for e=",v},ImageSize->200]
Row[plotit/@{0.8,1,1.5}]

4

The temperature in degrees Fahrenheit F is related to the temperature in degrees centigrade C by the formula F=(9/5)C+32. The temperature in Kelvin K is related to the centigrade temperature by C-273=K. Use Table to  make a table which gives all three temperatures for K=0,300 in 10 degree increments. Use TableForm with TableHeadings to make a nice printout.

Remove["Global`*"];
convertToF[k_]:=(9/5)convertToC[k]+32.;
convertToC[k_]:=k-273.;
Table[{k,convertToC[k],convertToF[k]},{k,0,300,10}];
TableForm[%,TableHeadings->{None,{"Kelvin","Centigrade","Fahrenheit"}}]

5

Consider the system of linear equations

a) Solve for x,y,z,w using Solve.

Remove["Global`*"];
eq1=I x+y+z-I w==0;
eq2=x+I y-I z+w==I;
eq3=x-y+(E^(I δ)) z+(E^(-I δ)) w==1;
eq4=x-I y+(E^(I δ)) z-(E^(-I δ)) w==0;
sol1=Solve[{eq1,eq2,eq3,eq4},{x,y,z,w}];

Print["x = ",x1=First[x/.sol1]];
Print["y = ",y1=First[y/.sol1]];
Print["z = ",z1=First[z/.sol1]];
Print["w = ",w1=First[w/.sol1]];

b) Extract the matrix of coefficients.  One convenient way to do this is to use Input-> Create Table/Matrix/Palette->Matrix.  Solve for x,y,z,w using Inverse. Convince yourself you get the same answers as in part a.

eq1=I x+y+z-I w==0;
eq2=x+I y-I z+w==I;
eq3=x-y+(E^(I δ)) z+(E^(-I δ)) w==1;
eq4=x-I y+(E^(I δ)) z-(E^(-I δ)) w==0;
eqs={eq1,eq2,eq3,eq4};

{b,coeff}=Normal[CoefficientArrays[eqs,{x,y,z,w}]];

Print["Matrix of coefficients=",MatrixForm[coeff]]

sol2=LinearSolve[coeff,b];

Print["x = ",x2=sol2[[1]]];
Print["y = ",y2=sol2[[2]]];
Print["z = ",z2=sol2[[3]]];
Print["w = ",w2=sol2[[4]]];
Print["Verifying same as part a..."];
If[Simplify[x1-x2]==0,Print["x verfied OK"],Print["x NOT verified"]];
If[Simplify[y1-y2]==0,Print["y verfied OK"],Print["y NOT verified"]];
If[Simplify[z1-z2]==0,Print["z verfied OK"],Print["z NOT verified"]];
If[Simplify[w1-w2]==0,Print["w verfied OK"],Print["w NOT verified"]];

6

The Bernoulli numbers (Subscript[B, n]) appear in many formulas you will find in tables of integrals. They are defined as the coefficients of powers of z in the expansion

(z/((e^z)-1))=1+(Subscript[B, 1])z+((Subscript[B, 2])/(2!))(z^2)+((Subscript[B, 3])/(3!))(z^3)...

Write a function bern[n] that returns the nth Bernoulli number. Use Series and Coefficient[expr,form] to extract the coefficient. Use the function and the Table command to make a table of the first 20 Bernoulli numbers. Compare with the built-in function BernoulliB.

f[z_]:=(z/((E^z)-1));
bern[n_]:=Coefficient[Normal[Series[f[z],{z,0,n}]],(z^n)]*Factorial[n]
t=Table[{n,bern[n],BernoulliB[n]},{n,1,20}];
TableForm[t,TableHeadings->{None,{"n","bern[n]","BernoulliB[n]"}}]

7

Consider the quadratic equation (x^2)+ R1 x +R2=0 , where R1 and R2 are independent random numbers uniformly distributed on [0,1] .

a) For each choice of random coefficients R1 and R2, the equation has two solutions, x1 and x2. Write a function printsol[neq_] which prints out  x1 and x2  for neq randomly generated equations. Use a Do loop. Execute the function for neq->15.

Remove["Global`*"];

neq=15;
printsol[]:=Module[{eq,sol,x},
eq=x^2+RandomReal[]x+RandomReal[]==0;
sol=Solve[eq,x];
sol=x/.sol;
Print["(Subscript[x, 1]) = ",sol[[1]],"   (Subscript[x, 2]) = ",sol[[2]]];
];

Do[printsol[],{n,1,neq}];

printsol[npts_] := Module[{x,x1, x2, s}, 
   Do[s = Flatten[Solve[x^2 + RandomReal[]*x + RandomReal[] == 0, x]]; 
     x1 = x /. First[s]; x2 = x /. Last[s]; Print["x1=  ", x1, "    x2= ", x2], 
    {i, 1, npts}]]

printsol[15]

b) Write a function reimxy[neq_] which computes the solutions x1 and x2 discussed above, and plots the location of each point in the complex plane.  Initialize a list xy={}, and use Re, Im, AppendTo, ListPlot,  and Do. Execute the function for neq->5000. What fraction of the solutions are Real, i.e. have zero imaginary part? Use a counter and If.

Remove["Global`*"];
xy={};
re=0;
neq=5000;

reimxy[]:=Module[{eq,sol,x,x1,x2},
eq=(x^2)+RandomReal[]x+RandomReal[]==0;
sol=Solve[eq,x];
sol=x/.sol;
x1=sol[[1]]; x2=sol[[2]];
If[Im[x1]==0&&Im[x2]==0,re=re+1];
xy=Append[xy,{Re[x1],Im[x1]}];
xy=Append[xy,{Re[x2],Im[x2]}];
];

Do[reimxy[],{n,1,neq}];

ListPlot[xy,PlotLabel->"Solution in the complex plane",AxesLabel->{"Real axes","Imaginary axes"}]
Print["Out of ",neq," equations, there were ",re," real solutions. This is ",re*100./neq,"% of all solutions"]

8

Find the center of mass of a solid half sphere of radius R . Recall that the position of the  center of mass is (Subscript[(OverVector[x]), cm])=(1/M){∫x ρ dV,∫y ρ dV,∫z ρ dV}, where M is the total mass and ρ is the density. Include a diagram which specifies your coordinate system. Go here for a reminder about how to do this.

Remove["Global`*"]
(*By symmetry, the x and y coordinates of center of mass must be located at the origin. Hence we only need to determine the z coordinate*)

volumeOfSphere=(∫_0^R)(∫_0^π)(∫_0^(2π))(r^2)Sin[θ] \[DifferentialD]ϕ\[DifferentialD]θ\[DifferentialD]r ; (*using spherical coordinates*)
volumeOfHemiSphere=(volumeOfSphere/2);
z=r Cos[θ]; (*by expressing z in spherical coordinates, see diagram*)
zcoord=(1/volumeOfHemiSphere)(∫_0^R)(∫_0^(π/2))(∫_0^(2π))z (r^2)Sin[θ] \[DifferentialD]ϕ\[DifferentialD]θ\[DifferentialD]r;
Print["Coordinates of center of mass are X=0, Y=0, Z=", zcoord];

9 (grads only)

The wavefunctions for the quantum mechanical harmonic oscillator are proportional to (e^(-((x^2)/2)))(Subscript[H, n])(x) where (Subscript[H, n])(x) are the Hermite polynomials defined by (Subscript[H, n])(x)=(((-1))^n)(e^(x^2))((∂_x)^n)(e^(-(x^2))). Write a function hermite[n] which returns the nth Hermite polynomial. Compare your answer with the built-in function HermiteH[n,x].Write a function psi[n] that returns the nth harmonic oscillator wave function.  Make sure that your wave functions are normalized, ie that  (∫_(-∞)^∞)(SuperStar[Ψ])Ψdx=1. Plot psi[13] and psi[13]^2, the probability density of finding the particle at x. (See the examples section of Plot for the trick for getting Plot to work properly.)

10 (grads only)

van der Waals attraction to a wall
van der Waals interactions between molecules produce an attractive interaction described by a potential of the form -((Subscript[C, 6])/(r^6)), where r is the distance between 2 molecules. Show that the interaction energy between a molecule and a half space ( a solid wall) is ((-(Subscript[C, 6])π ρ)/(6 (h^3)))if the molecule is a distance h from a wall with density of molecules ρ.
The molecule is at vector position x1={0,0,h}

x2={x,y,z} is the position of a typical molecule in the wall, and ranges over the region z<0 . If the number density of molecules in the solid wall is ρ, the potential energy of interaction between the molecule at x1 and all of the molecules in the wall is ∫∫∫((-ρ (Subscript[C, 6]))/(|x1-x2(|^6)))dx dy dz, where the range of integration is over the solid half space.