Physics 229& 100 Homework #3
Name: Nasser Abbasi

1.a) Show that (e^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 + I y;
{u, v} = {ComplexExpand[Re[(E^z)]], 
          ComplexExpand[
            Im[(E^z)]]};
\(\(showMsg[\*""\!\(E\^\(x + I\\ y\)\)"", isAnalytic[u, v]];\)\)

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]];

2

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

Remove["Global`*"];
ψ[x_, k_] := 
      A (E^(I k x)) + 
        B (E^(-I k x))

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

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

b) find the probability density ψ\(TraditionalForm\`\(ψ\^*\)\), 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\^*\)\)

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

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

{u, v} = {ComplexExpand[Re[sol]], 
          ComplexExpand[Im[sol]]};
FullSimplify[Sqrt[(u^2) + (v^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 .

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

Remove["Global`*"];
Vin = (E^(I ω t));
Vout = 
        Vin - current (I ω L);
Vr = ComplexExpand[(Vout/Vin)]

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

inductorZ = I ω L;
capacitorZ = 
        (1/(I ω c));
resistorZ  = R;
totalZ=(Sqrt[(R^2) + (((inductorZ - capacitorZ ))^2)])

vGround = 0;
voltageDrop  =  Vin  -  vGround;
current  =  Simplify[(voltageDrop/totalZ)]

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 (Subscript[ω, max]). What is (Subscript[ω, 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]

expr = expr /. Sqrt[any_] -> any

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

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

sol = ReplaceList[ω, sol]

Now find Maximum Frequency

params = {L -> 0.2, c -> (10^(-7)), 
          R -> 50};
Last[
        Position[nSol = ComplexExpand[Abs[sol /. params]], 
          Max[nSol]]];
ωMax = sol[[%]]

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.

\(\(Plot[mag /. params, {ω, 0, 90000}, 
        PlotLabel -> \*""(Subscript[V, r])(ω)"", 
        AxesLabel -> {"ω rad/sec", \*""|(Subscript[V, r])|""}];\)\)

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;
\(\(Plot[180\/π\ phase, {ω, 0, 90000}, Frame -> True, 
        FrameLabel -> {"ω rad/sec", \*""Phase (Subscript[V, r]) in degrees"", \*""Phase plot of (Subscript[V, r])(ω)"", 
            None}];\)\)

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

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 + I)))

    
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]];

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 + I 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};

  

c1 = 
        Integrate[realdXpart /. y -> 0, {x, 0, maxX}] + 
          Integrate[realdYpart /. x -> maxX, {y, 0, maxY}] + 
          I*(Integrate[imdXpart /. y -> 0, {x, 0, maxX}] + 
                Integrate[
                  imdYpart /. x -> maxX, {y, 0, 
                    maxY}]);

    
Print["Analytical Integral over C1 is \n", 
        c1];
Print["Numerical value of the above is ", 
        N[c1, 16]];

  

c2 = 
        Integrate[realdXpart /. y -> maxY, {x, 0, maxX}] + 
          Integrate[realdYpart /. x -> 0, {y, 0, maxY}] + 
          I*(Integrate[imdXpart /. y -> maxY, {x, 0, maxX}] + 
                Integrate[
                  imdYpart /. x -> 0, {y, 0, 
                    maxY}]);

    
Print["Analytical Integral over C2 is \n", 
        c2];
Print["Numerical value of the above is ", 
        N[c2, 16]];

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, -∞, ∞}]];

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

poles = 
        x /. Solve[Denominator[f[x]] == 0];
Print["Poles are ", 
        N[poles] // MatrixForm];

{xcoord, ycoord} = {Map[Re, poles], 
          Map[Im, poles]};
pts = Transpose[{xcoord, ycoord}];
ListPlot[pts, PlotStyle -> PointSize[ .02], 
        PlotLabel -> "Poles locations"];

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]]

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-(Subscript[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];

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)"];

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) &]]

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}];

\(\(ListPlot[data /. {x_, y_, z_} -> {1\/x, Log[z]}, 
        PlotStyle -> PointSize[ .015], 
        AxesLabel -> {\*""\!\(1\/x\)"", "Log[x]"}];\)\)

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]]];

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`*"];
howf = {(((E^(-((m (x^2) ω)/(2 ℏ)))) (((((m ω)/ℏ)))^(1/4)))/(π^(1/4))), (((Sqrt[2]) (E^(-((m (x^2) ω)/(2 ℏ)))) x (((((m ω)/ℏ)))^(3/4)))/(π^(1/4))), (((E^(-((m (x^2) ω)/(2 ℏ)))) (-2 + ((4 m (x^2) ω)/ℏ)) (((((m ω)/ℏ)))^(1/4)))/(2 (Sqrt[2]) (π^(1/4)))), (((E^(-((m (x^2) ω)/(2 ℏ)))) (-12 x (Sqrt[((m ω)/ℏ)]) + 8 (x^3) (((((m ω)/ℏ)))^(3/2))) (((((m ω)/ℏ)))^(1/4)))/(4 (Sqrt[3]) (π^(1/4)))), (((E^(-((m (x^2) ω)/(2 ℏ)))) (12 + ((16 (m^2) (x^4) (ω^2))/(ℏ^2)) - ((48 m (x^2) ω)/ℏ)) (((((m ω)/ℏ)))^(1/4)))/(8 (Sqrt[6]) (π^(1/4)))), (((E^(-((m (x^2) ω)/(2 ℏ)))) (120 x (Sqrt[((m ω)/ℏ)]) - 160 (x^3) (((((m ω)/ℏ)))^(3/2)) + 32 (x^5) (((((m ω)/ℏ)))^(5/2))) (((((m ω)/ℏ)))^(1/4)))/(16 (Sqrt[15]) (π^(1/4)))), (((E^(-((m (x^2) ω)/(2 ℏ)))) (-120 + ((64 (m^3) (x^6) (ω^3))/(ℏ^3)) - ((480 (m^2) (x^4) (ω^2))/(ℏ^2)) + ((720 m (x^2) ω)/ℏ)) (((((m ω)/ℏ)))^(1/4)))/(96 (Sqrt[5]) (π^(1/4))))};

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

Map[
      Integrate[ 
          ComplexExpand[ 
            Conjugate[#] #], {x, -∞, ∞}, 
          Assumptions -> {Re[((m ω)/ℏ)] > 0 && 
                m ≥ 0 && ω ≥ 
                  0 && ℏ ≥ 0}] &, howf]

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

11 (grads)

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

Remove["Global`*"];
f[x_] := (E^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  π}]

c) Show that the integral of ((e^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_] := ((E^z)/(z^2));
poles = Solve[Denominator[g[z]] == 0];
Series[g[z], {z, 0, 5}];
residue = Coefficient[%, (1/z)];
sol = Abs[2.  π I residue]

Note answer matches part (a)

Well done!