Physics 229& 100 Homework #3
Name: Nasser Abbasi

1.a) Show that ^z , where z=x+i y, is analytic, i.e. it satisfies the Cauchy-Riemann conditions.

Remove["Global`*"] ; (* Define 2 functions to use for this problem *)

isAnalytic[u_, v_] := If[D[v, y] === D[u, x] && -D[u, y] === D[v, x], True, False] ;

showMsg[name_, result_] := If[result, Print[name<> " is Analytic function"], Print[name<> " not analytic function"]] ;

z = x +  y ;

{u, v} = {ComplexExpand[Re[^z]], ComplexExpand[Im[^z]]} ;

                x + \y showMsg[              , isAnalytic[u, v]] ;

        x +  y                is Analytic function

b)  Many familiar functions are not analytic. Show that the functions Conjugate[z] and Abs[z] do not satisfy the Cauchy Riemann conditions.

{u, v} = {ComplexExpand[Re[Conjugate[z]]], ComplexExpand[Im[Conjugate[z]]]} ;

showMsg["Conjugate[z]", isAnalytic[u, v]] ;

{u, v} = {ComplexExpand[Re[Abs[z]]], ComplexExpand[Im[Abs[z]]]} ;

showMsg["Abs[z]", isAnalytic[u, v]] ;

Conjugate[z] not analytic function

Abs[z] not analytic function

2

A 1D wave function is of the form ψ=A ^( k x) +B ^(- k x)

Remove["Global`*"] ;

ψ[x_, k_] := A ^( k x) + B ^(- k x)

a) find the probability density ψψ^*, assuming that A, B and k are real.

sol = Simplify[ComplexExpand[ψ[x, k] * Conjugate[ψ[x, k]]]]

A^2 + B^2 + 2 A B Cos[2 k x]

b) find the probability density ψψ^*, assuming that k is real but A and B are complex

sol = FullSimplify[ComplexExpand[ ψ[x, k] * Conjugate[ψ[x, k]], {A, B}]] ;

sol/.Conjugate[p_] ->p^*

Abs[A]^2 + Abs[B]^2 + B ^(-2  k x) A^* + A ^(2  k x) B^*

c) find the particle current density, which is given by (( ℏ)/(2 m))(ψ∇ψ^*-ψ^*∇ψ), assuming that A,B and k are real.( in 1D ∇ is simply d/dx)

sol = ( h)/(2 m) (ψ[x, k]   D[Conjugate[ψ[x, k]], x] - Conjugate[ψ[x, k]] D[ψ[x, k], x]) ;

sol = Simplify[ComplexExpand[sol]]

((A - B) h k Sin[k x] (- (A - B) Cos[k x] + (A + B) Sin[k x]))/m

{u, v} = {ComplexExpand[Re[sol]], ComplexExpand[Im[sol]]} ;

FullSimplify[Sqrt[u^2 + v^2]]

((A - B)^2 h^2 k^2 (A^2 + B^2 - 2 A B Cos[2 k x]) Sin[k x]^2)/m^2^(1/2)

3.  Complex numbers occur naturally in circuit theory. Recall that the impedance of an inductor is IωL and the impedance of a capacitor is 1/(I ω C) , and that V=I Z. Consider the following circuit, where Vin is a unit magnitude sinusoidal input e^iωt .

[Graphics:HTMLFiles/index_37.gif]

a) Find the symbolic expression for the ratio Vout/Vin.

Remove["Global`*"] ;

Vin = ^( ω t) ;

Vout = Vin - current ( ω L) ;

Vr = ComplexExpand[Vout/Vin]

- current L ω Cos[t ω] + Cos[t ω]^2 - current L ω Sin[t ω] + Sin[t ω]^2

Calculating the current:
Use the series RLC Z law, which says that total impedence = (R^2 + (Inductor_Z - Capacitor_Z)^2)^(1/2), then current=Vin/total_Impedence

inductorZ =  ω L ; capacitorZ = 1/( ω c) ; resistorZ = R ; totalZ = (R^2 + (inductorZ - capacitorZ )^2)^(1/2)

(R^2 + (/(c ω) +  L ω)^2)^(1/2)

vGround = 0 ;

voltageDrop = Vin - vGround ;

current = Simplify[voltageDrop/totalZ]

^( t ω)/(R^2 - (1 + c L ω^2)^2/(c^2 ω^2))^(1/2)

b) Vr=Vout/Vin is a complex quantity. Find an expression for the magnitude of Vr  in terms of L,R,C. The magnitude of Vr has a maximum at a certain frequency ω_max. What is ω_max in terms of R,L,C?

Vr

mag = Simplify[ComplexExpand[Abs[Vr], TargetFunctions→ {Re, Im}]]

Now to find where mag is maximum. The above is maximum when the denominator is zero.

expr = mag/.Sqrt[any_] →any

expr = Denominator[expr]

(R^2 - (1 + c L ω^2)^2/(c^2 ω^2))^2^(1/2)

expr = expr/.Sqrt[any_] →any

(R^2 - (1 + c L ω^2)^2/(c^2 ω^2))^2

expr = expr/.(a_ - b_)^2→ (a - b)

R^2 - (1 + c L ω^2)^2/(c^2 ω^2)

sol = FullSimplify[Flatten[Solve[expr == 0, ω]]]

sol = ReplaceList[ω, sol]

{-(R + (-4 L + c R^2)^(1/2)/c^(1/2))/(2 L), (R - (-4 L + c R^2)^(1/2)/c^(1/2))/(2 L), -(R - (-4 L + c R^2)^(1/2)/c^(1/2))/(2 L), (R + (-4 L + c R^2)^(1/2)/c^(1/2))/(2 L)}

Now find Maximum Frequency

params = {L->0.2, c->10^(-7), R->50} ;

Last[Position[nSol = ComplexExpand[Abs[sol/.params]], Max[nSol]]] ;

ωMax = sol[[%]]

{(R + (-4 L + c R^2)^(1/2)/c^(1/2))/(2 L)}

c) For the values L=0.2, C=10^(-7), R=50, make a plot of the magnitude of Vr as a function of ω over the physically relevant range.

[Graphics:HTMLFiles/index_75.gif]

d) Make a plot of the argument of Vr as a function of frequency for the same parameter values.

phase = Simplify[ComplexExpand[Arg[Vr], TargetFunctions→ {Re, Im}]]

phase = phase/.params ;

[Graphics:HTMLFiles/index_80.gif]

4. Consider  the function of a complex variable z given by  f(z)= 1/(z - (1 + )) .
Calculate Underscript[∫, C] f(z) dz , for the paths  C1 and C2 shown in the figure. Do you expect them to be the same?

[Graphics:HTMLFiles/index_83.gif]

a) first do the integrals along C1 and C2 numerically. See numerical contour integrals.

Remove["Global`*"] ;

maxX = 4 ; maxY = 3 ;

f[z_] := 1/(z - (1 + ))

c1 = Join[Table[i, {i, 0, maxX, .01}], Table[maxX + i, {i, 0, maxY I, 0.01 I}]] ;

sol = NIntegrate[f[z], Evaluate[Join[{z}, c1]]] ;

Print["Integral along C1 = ", sol, "\nWhose magnitude is ", Abs[sol]] ;

c2 = Join[Table[i, {i, 0, maxY I, 0.01 I}], Table[maxY * I + i, {i, 0, maxX, .01}]] ;

sol = NIntegrate[f[z], Evaluate[Join[{z}, c2]]] ;

Print["Integral along C2 = ", sol, "\nWhose magnitude is ", Abs[sol]] ;

General :: spell1 : Possible spelling error: new symbol name \"maxY\" is similar to existing symbol \"maxX\".  More…

Integral along C1 = 0.935901 + 2.9442 <br />Whose magnitude is 3.08937

Integral along C2 = 0.935901 - 3.33899 <br />Whose magnitude is 3.46767

The answer is not the same, since the function f(z) is not analytical over the closed region, as there is a pole inside the rectangular area.

b) Now calculate each integral analytically. Do this by setting z=x+i y , dz=dx+i dy, and expanding f(z) dz into terms which depend only on the real variables x and y. The path will determine the values of x,y,dx,dy along it.

z = x +  y ;

dz = D[z, x] dx + D[z, y] dy ;

g = ComplexExpand[f[z] dz] ;

{re, im} = {ComplexExpand[Re[g]], ComplexExpand[Im[g]]} ;

{re, im} = {Collect[re, {dx, dy}], Collect[im, {dx, dy}]} ;

{realdXpart, realdYpart} = re/.dx (pat1_) + dy (pat2_) → {pat1, pat2} ;

{imdXpart, imdYpart} = im/.dx (pat1_) + dy (pat2_) → {pat1, pat2} ; (*now do C1 integral*)

Print["Analytical Integral over C1 is \n", c1] ;

Print["Numerical value of the above is ", N[c1, 16]] ; (*now do C2 integral *)

Print["Analytical Integral over C2 is \n", c2] ;

Print["Numerical value of the above is ", N[c2, 16]] ;

General :: spell1 : Possible spelling error: new symbol name \"realdYpart\" is similar to existing symbol \"realdXpart\".  More…

General :: spell1 : Possible spelling error: new symbol name \"imdYpart\" is similar to existing symbol \"imdXpart\".  More…

Analytical Integral over C1 is <br /> (π/4 + ArcTan[1/3] + ArcTan[2/3] + ArcTan[3]) + 1/2 Log[13/10] + Log[5]/2

Numerical value of the above is 0.935901088450796 + 2.944197093739912 

Analytical Integral over C2 is <br /> (-π/4 - ArcTan[1/2] - ArcTan[3/2] - ArcTan[2]) + 1/2 Log[5/2] + 1/2 Log[13/5]

Numerical value of the above is 0.935901088450796 - 3.338988213439674 

5. Consider the integral              ∫_ (-∞)^∞1/(x^4 + 1)dx

a) Find the value using the Integrate command.

Remove["Global`*"] ;

f[x_] := 1/(x^4 + 1) ;

Print["Result = ", Integrate[f[x], {x, -∞, ∞}]] ;

Result = π/2^(1/2)

b) Locate the poles (the places where the function blows up).

poles = x/.Solve[Denominator[f[x]] == 0] ;

Print["Poles are ", N[poles]//MatrixForm] ;

Poles are  ( {{-0.707107 - 0.707107 }, {0.707107 + 0.707107 }, {0.707107 - 0.707107 }, {-0.707107 + 0.707107 }} )

{xcoord, ycoord} = {Map[Re, poles], Map[Im, poles]} ;

pts = Transpose[{xcoord, ycoord}] ;

ListPlot[pts, PlotStyle->PointSize[.02], PlotLabel->"Poles locations"] ;

General :: spell1 : Possible spelling error: new symbol name \"ycoord\" is similar to existing symbol \"xcoord\".  More…

[Graphics:HTMLFiles/index_128.gif]

c) Find the residues of the function at the poles in the upper half plane using Residue.

upperPoles = Select[poles, Im[#] >0&] ;

residue = Map[Residue [f[x], {x, #}] &, upperPoles] ;

Print["Residues for poles in upper plane are ", ComplexExpand[residue]]

General :: spell1 : Possible spelling error: new symbol name \"residue\" is similar to existing symbol \"Residue\".  More…

Residues for poles in upper plane are  {-(1/4 + /4)/2^(1/2), (1/4 - /4)/2^(1/2)}

d) Find the residues of the function at the poles in the upper half plane by doing a series expansion around the poles using Series. The residues are the coefficients of  1/(x-x_pole)

s = Map[Series[f[x], {x, #, 5}] &, upperPoles] ;

r = ComplexExpand[MapThread[Coefficient[#1, 1/(x - #2)] &, {s, upperPoles}]] ;

Print["Residues in upper plane found using Series method are ", r] ;

(*Notice that answer matches part c *)

Residues in upper plane found using Series method are  {-(1/4 + /4)/2^(1/2), (1/4 - /4)/2^(1/2)}

e) Calculate the integral using the fact that it equals 2πi ( sum of residues). Conclude that a) was much easier than b+d+e

Print["Integral value is ", ComplexExpand[2 π I Total[residue]]] ;

Print["Answer matches part (a)"] ;

Integral value is π/2^(1/2)

Answer matches part (a)

Answer using part(a) method was simpler since we did not have to determine the poles, followed by the residues around the poles which will be inside the region of integration (we also had which poles to pick).

6. Here is a list of complex numbers.

cnum = Table[Random[] + I ( Random[] - 0.5), {i, 1, 20}]

a) Produce a list with real and imaginary parts of each element exchanged. Use a pattern.

flip = Transpose[{Re[cnum], Im[cnum]}]/.{p1_, p2_} -> {p2, p1} ;

newList = Map[First, flip] + I * Map[Last, flip]

b) In one line of code, produce a list of the elements of cnum which have negative imaginary parts. Use Select and a pure function.

Select[cnum, Im[#] <0&]

c) Here is a longer list of complex numbers. Note the ; at the end, so it is not printed out. In one  line of code, count the number of elements whose magnitude is greater than 1/2.

cnum2 = Table[Random[] + I Random[], {i, 1, 250}] ;

Length[Select[cnum2, Abs[#] >1/2&]]

200

7  Here is a list which represents data in {x,y,z} triplets. In a single line of code, produce a plot of 1/x versus Log[z].

Remove["Global`*"] ;

data = Table[{Random[], 3 * Random[], 1/Random[]}, {i, 1, 20}] ;

[Graphics:HTMLFiles/index_157.gif]

8 a) The list rth contains the {r,θ} values for a number of points in the plane. In one line of code, convert the list into a list of the {x,y} coordinates of the points. Use a pattern.

Remove["Global`*"] ;

rth = Table[{Random[] * 5, Random[] * 2 Pi}, {15}]

xy = rth/.{r_, θ_} → {r Cos[θ], r Sin[θ]}

(*ps . Do I need to worry about which quadrant it is ? *)

b) rvec is a list of vectors. In one line of code, convert it into a list of normalized vectors. Use Map and a pure function.

rvec = Table[{Random[], 2 * Random[], 1/Random[]}, {15}]

Map[#/Norm[#] &, rvec]

9. The following table is a list of the positions of a ball that has been thrown in the air.

Remove["Global`*"] ;

traj = Table[{4, 6} t - 9.8 {0, 1} t^2/2, {t, 0, 1.2, .05}]

Construct a picture which shows a circle of radius of  0.1 at each of these locations. Use Circle, Graphics and Show and Map and a pure function.  

Show[Graphics[Map[Circle[#, .1] &, traj]]] ;

[Graphics:HTMLFiles/index_172.gif]

10.(grads)  Here is a list of harmonic oscillator wave functions:

Please note: I am not signed for the graduate course, but making an attempt to solve this.

Remove["Global`*"] ;

Show that they are properly normalized. Use Map and a pure function.

{1, 1, 1, 1, 1, 1/4 (-11 + 15 Cos[2 Arg[m ω]]), 1}

Map[Integrate[ ComplexExpand[ Conjugate[#] #], {x, -∞, ∞}, Assumptions→ {ω>0, ℏ>0, m>0}] &, howf]

{1, 1, 1, 1, 1, 1, 1}

11 (grads)

a) Using Residues to compute an integral that Mathematica gets wrong.  Consider the integral ∫_0^(2π) ^Cos[x]Cos[ x-Sin[x]]dx . Use Integrate to compute the integral.

Remove["Global`*"] ;

f[x_] := ^Cos[x] Cos[x - Sin[x]] ;

Integrate[f[x], {x, 0, 2π}]

π

b) Now do the integral numerically and convince yourself that  the analytic answer (in Version≤ 5.2) is wrong.

NIntegrate[f[x], {x, 0, 2π}]

(*we see the answer in (a) is therefor wrong*)

6.28319

c) Show that the integral of ^z/z^2 around the unit circle is related to the integral in part a). Use the residue theorem applied to this relation to compute the value of the integral.

g[z_] := ^z/z^2 ;

poles = Solve[Denominator[g[z]] == 0] ;

Series[g[z], {z, 0, 5}] ;

residue = Coefficient[%, 1/z] ;

sol = Abs[2.π I residue]

6.28319

Note answer matches part (a)

Well done !


Created by Mathematica  (January 27, 2007) Valid XHTML 1.1!