Homework 3I — Solution to Individual Exercises

3.016 Mathematical Methods for Materials Scientists and Engineers
S. M. Allen
October 31, 2005

Clear["Global`*"]

Individual Exercise I3-1: (Kreyszig MCG 6.4, p. 77)

Using the matrices provided

A = {{1, 3, 2}, {3, 5, 0}, {2, 0, 4}}
B = {{0, 2, 1}, {-2, 0, -3}, {-1, 3, 
          0}}
c = {{1}, {0}, {-2}}
d = {{3}, {1}, {2}}

Going to it:

B . c

% // MatrixForm

Transpose[c] . A . c

A . B . c

% // MatrixForm

Flatten[A . c] . Flatten[Transpose[B . d]]

Transpose[d] . B . d

B . d

As to why the result is zero, look at the eigensystem of B

Eigensystem[B]

Note that the third eigenvector is equal to d and its corresponding eigenvalue is 0. Multiplication of a matrix by its eigenvector returns a parallel vector scaled in length by the corresponding eigenvalue.

Individual Exercise I3-2: (Kreyszig MCG 7.10, p. 87)

Clear[A, B]

A  =  {{3, 4}, {4, -3}}
P  = {{-4, 2}, {3, -1}}

B = Inverse[P] . A . P

Note that by definition, B is similar to A (B was obtained from A by a similarity transformation).

EigA = Eigensystem[A]
EigB = Eigensystem[B]

Obviously, the eigenvalues are the same.

Inverse[P] . EigA[[2, 1]]
Inverse[P] . EigA[[2, 2]]

We see also that multiplication of Inverse[P] by an eigenvector of A yields an eigenvector of B (see Kreyszig AEM Theorem 1, p. 392).

Individual Exercise I3-3: (Kreyszig MCG 7.12, p. 87)

Clear[A]

Define the matrix A by

A  =  {{5, 4}, {1, 2}}

EiVecA = Eigenvectors[A]

Make a square matrix of eigenvectors, each eigenvector a column of the matrix.

X = Transpose[EiVecA]

DMatrix = Inverse[X] . A . X

Eigenvalues[A]

Homework 3G — Solution to Group Exercises

3.016 Mathematical Methods for Materials Scientists and Engineers
S. M. Allen
November 14, 2005

Clear["Global`*"]

Homework 3G — Solution to Group Exercises

3.016 Mathematical Methods for Materials Scientists and Engineers
S. M. Allen
November 14, 2005

Clear["Global`*"]

Group Exercise G3-1: Edge Dislocation

1.

Enter the stress matrix in cartesian coordinates; use Replace to convert to Polar form.

(Subscript[σ, xx]) = ((G b)/(2  π (1 - ν))) ((-y (3 (x^2) + (y^2)))/((((x^2) + (y^2)))^2));
(Subscript[σ, xy]) = ((G b)/(2  π (1 - ν))) ((x ( (x^2) - (y^2)))/((((x^2) + (y^2)))^2));
(Subscript[σ, yx]) = ((G b)/(2  π (1 - ν))) ((x ( (x^2) - (y^2)))/((((x^2) + (y^2)))^2));
(Subscript[σ, yy]) = ((G b)/(2  π (1 - ν))) ((y ( (x^2) - (y^2)))/((((x^2) + (y^2)))^2));

SigMatCart = {{(Subscript[σ, xx]), (Subscript[σ, xy])}, {(Subscript[σ, yx]), (Subscript[σ, yy])}};
MatrixForm[%]

(Subscript[σ, rr]) = (Subscript[σ, xx]) /. {x -> r Sin[θ], 
            y -> r Cos[θ]};
(Subscript[σ, r])θ = (Subscript[σ, xy]) /. {x -> r Sin[θ], 
            y -> r Cos[θ]};
(Subscript[σ, θ])r = (Subscript[σ, yx]) /. {x -> r Sin[θ], 
            y -> r Cos[θ]};
(Subscript[σ, θ])θ = (Subscript[σ, yy]) /. {x -> 
              r Sin[θ], y -> r Cos[θ]};

SigMatPolar = {{(Subscript[σ, rr]), (Subscript[σ, r])θ}, {(Subscript[σ, θ])r, (Subscript[σ, θ])θ}};
MatrixForm[%]

2.

HydroPress[x_, 
        y_] := (2/3) ((Subscript[σ, xx]) + (Subscript[σ, yy]))

ContourPlot[
      HydroPress[x, y] /. {G -> 1, b -> 1, ν -> 1/3}, {x, -4, 
        4}, {y, -4, 4}, Contours -> 50]

3.

Principal axes are given by eigenvectors of the stress matrix

evs = Eigenvectors[SigMatCart]

The orientation of principal axes can be illustrated by plotting the first eigenvector(normalized) as a function of position using PlotVectorField in the Graphics package:

NormEV1[x_, y_] := 
      FullSimplify[
        evs[[1]]/Sqrt[evs[[1, 1]]^2 + evs[[1, 2]]^2], 
        Assumptions -> {x ∈ Reals, 
            y ∈ Reals}]

<< Graphics`PlotField`

PlotVectorField[
      NormEV1[x, y], {x, -4, 4}, {y, -4, 4}]

Note that there are singularities along the line x = -y, but Mathematica will continue and make the remainder of the plot anyway...

Note how the vectors point toward each other along x = -y in the second quadrant, and away from each other along x = -y in the fourth quadrant.

4.

From page 61 of the class notes (Lecture 10), it is seen that the maximum shear stress is equal to the radius of Mohr's circle and thus is 1/2 the difference between the eigenvalues of the stress matrix.

PrincipalStresses[x_, y_] = Eigenvalues[SigMatCart]

MaxShear[x_, 
        y_] = (1/2) (PrincipalStresses[x, 
                y][[1]] - PrincipalStresses[x, 
                y][[2]])

ContourPlot[
      MaxShear[x, y] /. {G -> 1, b -> 1, ν -> 1/3}, {x, -4, 
        4}, {y, -4, 4}, Contours -> 50]

Note that the maximum shear stress is concentrated in the dislocation's glide plane. (See The Structure of Materials to learn more about dislocations and glide planes.)

Group Exercise G3-2: Entropy of Small System

Credits go to the Woolf group: Jill A. Rowehl, Jin Suk Kim, Katherine Hartman, Talia Gershon ; and to Prof. Carter, who assisted them.

Each atom can have its electron in either the n = 1 or the n = 2 states. The energy of an electron in one of these states depends on n as

Energy[n_] := (Subscript[ϵ, 0]) (1 - (1/(n^2)))

1 and 2.

The distinct quantum states for the three-atom system can be represented

States = Tuples[{1, 2}, 3]

Make a list of energies of each of these eight states with this double summation

EnergyDist = {};
EnergyList = 
      Do[AppendTo[EnergyDist, 
          Sum[1 - (1/(States[[i, j]]^2)), {j, 
              Length[States[[i]]]}]], {i, Length[States]}]
EnergyDist

UniqueEnergies = Union[EnergyDist]

Here is a function that will calculate the degeneracy of energy states in a distribution:

EntropyVec[Elist_] := 
      Module[{uniques = Union[Elist], counts = {}}, 
        For[i = 1, i ≤ Length[uniques], i++, 
          AppendTo[
            counts, {Count[Elist, 
                uniques[[i]]]}]]; 
Return[
          counts]]

Degeneracy = Flatten[EntropyVec[EnergyDist]]

Here is the entropy of each state

SOverBoltzList =  Log[Degeneracy]

Make a function that will combine two lists and use it to make ordered pairs of {energy,entropy} which can then be plotted with ListPlot:

CombList[list1_, list2_] := 
      Module[{big = {}}, 
        Do[AppendTo[big, {list1[[i]], list2[[i]]}], {i, 
            Length[list1]}]; Return[big]]

SvsE = CombList[UniqueEnergies, SOverBoltzList]

ListPlot[SvsE, PlotJoined -> True, 
      AxesLabel -> {"Energy", "S/k"}]

3.

Proceed as above but for atoms that can have three energy levels n = 1, 2, or 3:

States = Tuples[{1, 2, 3}, 3];

EnergyDist = {};
EnergyList = 
      Do[AppendTo[EnergyDist, 
          Sum[1 - (1/(States[[i, j]]^2)), {j, 
              Length[States[[i]]]}]], {i, Length[States]}]

UniqueEnergies = Union[EnergyDist];

Degeneracy = Flatten[EntropyVec[EnergyDist]];

SOverBoltzList =  Log[Degeneracy];

SvsE = CombList[UniqueEnergies, SOverBoltzList];

ListPlot[SvsE, PlotJoined -> True, 
      AxesLabel -> {"Energy", "S/k"}]

3.

Proceed as above but for 5 atoms that can have three energy levels n = 1, 2, 3:

States = Tuples[{1, 2, 3}, 5];

EnergyDist = {};
EnergyList = 
      Do[AppendTo[EnergyDist, 
          Sum[1 - (1/(States[[i, j]]^2)), {j, 
              Length[States[[i]]]}]], {i, Length[States]}]

UniqueEnergies = Union[EnergyDist];

Degeneracy = Flatten[EntropyVec[EnergyDist]];

SOverBoltzList =  Log[Degeneracy];

SvsE = CombList[UniqueEnergies, SOverBoltzList];

ListPlot[SvsE, PlotJoined -> True, 
      AxesLabel -> {"Energy", "S/k"}]

4 (Extra Credit).

Because dU = T dS - P dV, for fixed volume (Subscript[((dS/dU)), V]) = 1/T and the slope of the plots above of S/k vs. Energy should be positive. Lack of accounting for high-energy states underestimates the entropy and at high temperatures S/k necessarily goes through a maximum and then falls to zero. The maximum in S/k would correspond to infinite temperature, and beyond the maximum would correspond to negative temperatures. This is unphysical.

Group Exercise G3-1: Hamiltonian

solution by WCC

Background Material

Here is the free electron Hamiltonian

Ham[
        phi_]  := (( -(ℏ^2))/(2 m)) D[
          phi, {x, 2}]

Normalize x with the characteristic distance L, and define the normalized Hamiltonian operator as
HamNorm = \(TraditionalForm\`H\/\(\(\(\ \)\(-ℏ\^2\)\)\/\(2\ m\ L\^2\)\)\) and z = x/L

HamNorm[phi_]  :=  -D[phi, {z, 2}]

Define the approximating basis functions; the Hamiltonian matrix elements, and the normalizing elements

phi[i_]  :=  z^i (1 - z)^i

H[i_,  j_]  := 
      Evaluate[ 
        Integrate[phi[i] HamNorm[phi[j]], {z, 0, 1}, 
          Assumptions -> 
            i  ∈  Integers  &&  j  ∈  Integers  &&  
              i  >  1  &&  j  >  1]]

F[i_,  j_]  := 
      Evaluate[ 
        Integrate[phi[i] phi[j], {z, 0, 1}, 
          Assumptions -> 
            i  ∈  Integers  &&  j  ∈  Integers  &&  
              i  >  1  &&  j  >  1]]

Part 1

Simplify[H[i, j], 
      Assumptions -> 
        i  ∈  Integers  &&  j  ∈  Integers  &&  
          i  >  1  &&  j  >  1]

We see that it is not at all obvious that H[i,j] == H[j,i]
but it is plausible given the explicit values below

H[2, 3]

H[3, 2]

To see if this is symmetric, try and integrate by parts:

dv  =  HamNorm[phi[j]]

v  =  
      Integrate[dv, z, 
        Assumptions ->  
          j  ∈  Integers  &&  j  >  1]

u  =  phi[i]

du  =  D[phi[i], z]

HNormNew =  
      Simplify[(u v) /. z -> 1, 
          Assumptions  -> 
            i  ∈  Integers  &&  j  ∈  Integers  &&  
              i  >  1  &&  
              j  >  1] - 
Simplify[(u v) /. z -> 0, 
          Assumptions  -> 
            i  ∈  Integers  &&  j  ∈  Integers  &&  
              i  >  1  &&  j  >  1] - 
Integrate[
          v du, {z, 0, 1}, 
          Assumptions  -> 
            i  ∈  Integers  &&  j  ∈  Integers  &&  
              i  >  1  &&  j  >  1]

Which is obviously symmetric with respect to interchange of i and j.
Test for consistency with above result

HNormNew /. {i -> 3, j -> 2}
HNormNew /. {j -> 2,  i ->  3}

Part 2

Write a function to calculate E(c1,c2,... cn)

Eform[
        approxlevel_Integer]  :=  
Sum[c[index1] c[index2] HNormNew /. {i -> index1, j -> index2}, {index1, approxlevel}, {index2, approxlevel}]/Sum[c[index1] c[index2] F[index1, index2], {index1, approxlevel}, {index2, approxlevel}]

Eform[2]

Find the minimum

The hard way:

eqs[M_]  :=  
      Table[D[Eform[M], c[k]] ==  0, {k, M}]
vars[M_]  :=  Table[c[k], {k, M}]

eqs[2]

Clear[sols]
sols[M_]  :=  Solve[eqs[M], vars[M]]

sols[2]

Simplify[Eform[2] /. sols[2]]

The medium way:

solsmed[M_]  :=  Minimize[Eform[M], vars[M]]

solsmed[2]

Hmm, are these two results the same?

FullSimplify[ ((28 (95 + 8 (Sqrt[133])))/(266 + 23 (Sqrt[133])))]

Yes...  And now the easy way

solseasy[M_]  :=  NMinimize[Eform[M], vars[M]]

solseasy[2]

\(The\ exact\ ground\ \ state\ energy\ is\ E\  = \ \(\(-\((2\ π)\)\^2\)\ ℏ\^2\)\/\(8\ m\ L\^2\)\ which\ corresponds\ to\ \ π\^\(\(2\)\(\ \)\)\ in\ our\ normalized\ energy\ units, \ 
    therefore\ we\ should\ compare\ π\^\(\(2\)\(\ \)\)\ with\ \(28\ \((95 + 8\ \@133)\)\)\/\(266 + 23\ \@133\)\)

N[((28 (95 + 8 (Sqrt[133])))/(266 + 23 (Sqrt[133])))]
N[(π^2)]

Part 3

The equation between E((Overscript[c,  ->]) ), H((Overscript[c,  ->]) ), and    can be rewritten as
E((Overscript[c,  ->]) ) F((Overscript[c,  ->]) ) = H((Overscript[c,  ->]) )

in our approximating scheme, this becomes a relationship between two quadratic forms:
E((Overscript[c,  ->]) )    (Overscript[c,  ->]) •(Underscript[F, ─])\(TraditionalForm\`\(\(• 
             c\&→\)\&\ \ \ \ \ \ \ \)\) =   (Overscript[c,  ->]) •(Underscript[H, ─])\(TraditionalForm\`\(\(• 
             c\&→\)\&\ \ \ \ \ \ \ \)\) 

Take the gradient with respect to (Overscript[c, →]) and use the extremum condition ∇E = (Overscript[0, →])
E((Overscript[c,  ->]) )  (Underscript[F, ─])\(TraditionalForm\`\(\(• 
             c\&→\)\&\ \ \ \ \ \ \ \)\) = (Underscript[H, ─])\(TraditionalForm\`\(\(• 
             c\&→\)\&\ \ \ \ \ \ \ \)\) 
or as a Homogeneous equation,
 ((Underscript[H, ─]) - E((Overscript[c,  ->]) )  (Underscript[F, ─])) \(TraditionalForm\`\(\(• 
             c\&→\)\&\ \ \ \ \ \ \ \)\)  = 0
 which has solutions if and only if det( ((Underscript[H, ─]) - E((Overscript[c,  ->]) )  (Underscript[F, ─])) = 0

Part 4

The larger eigenvalue is the same as the the maximum above:

NMaximize[Eform[2], vars[2]]

it doesn't correlate to the other energy levels at all.

Part 5

\(To\ make\ comparisons, \ 
    let' s\ divide\ the\ difference\ between\ the\ solutions\ by\ exact\ the\ ground\ state\ \((normalized)\)\ energy\ π\^2\ and\ write\ a\ function\ that\ returns\ this\ \(\(ratio\)\(:\)\)\)

solsNormeasy[
        M_]  := (( NMinimize[Eform[M], vars[M]][[1]] - (π^2) )/((π^2) ))

compareData  =  Table[{m, solsNormeasy[m]}, {m, 2, 20}]

<< Graphics`Graphics`

DisplayLater  =  
        DisplayFunction  ->  Identity;
DisplayNow  =  
        DisplayFunction  ->  $DisplayFunction;

p1  =  
      LogListPlot[compareData, PlotStyle -> {PointSize[0.05], Hue[1]},  
        DisplayLater]

p2  =  
      LogListPlot[compareData, PlotJoined -> True, DisplayLater]

Show[p1, p2, DisplayNow]

We hit machine precision difficulties at around M=5...  There is no need to go further unless we take measures to improve computation accuracy.