Homework 1 — Solution

3.016 Mathematical Methods for Materials Scientists and Engineers
S. M. Allen
September 30, 2005

Individual Exercise I-1:

First, recognize that you will need to use ListPlot to make the final plot of your data.  ListPlot produces a plot from (x,y) pairs that are stored in a Table.  The program you need to write has to produce output that will be written to a Table with the proper format.
The function Random[Integer] is used here; it returns either 0 or 1 with probability 1/2.
Quoting Charles Cantrell: "All of the problem is contained in one line of input.  First I create a table of data members that contains the trial number and the number of required steps to reach the endpoint.  Then I plot the data using the ListPlot command.  The While loop runs until the position is either 100 or -100, and the Table command is run 100 times--giving us 100 data points."

data = 
        Table[{i, location = 0; timestep = 0; 
            While[Abs[location] < 100, 
              If[Random[Integer] == 0, location = location + 1, 
                location = location - 1]; timestep = timestep + 1]; 
            timestep}, {i, 100}];
ListPlot[data, 
      AxesLabel -> {"Trial number", "TimeSteps"}]

Individual Exercise I-2:

Note: Some features of Omar Fabian's work were incorporated in the solution to this exercise.

1.

Input the LJ potential, differentiate and set derivative == 0 to find minima.

LJ = (a/r)^12 - (b/r)^6

MinimaOfLJ  =  Solve[D[LJ, r] ==  0, r]

Note that there are six roots but only one of them is "physical" — i.e., real and positive (the second root). Use this to set rmin. Note how Replace is used to extract the value of the second root.

Minr = r /.  MinimaOfLJ[[2]]

2.

Replace r with rmin in LJ to get EMin = LJ(rmin).

MinE =   LJ /. r -> Minr

3.

The Solve function can be used to eliminate specific variables from a set of equations.  Here we have the equation for rMin in terms of (a,b) and the equation for EMin in terms of (a,b).  Solve is used to find expressions for a and b in terms of rMin and EMin.

Clear[rMin, EMin]

NewVars = 
      Solve[{rMin == Minr, EMin == MinE}, {a, b}]

All that remains is to Replace a and b with these new variables.

LJModified = LJ /. NewVars[[1]]

4.

The force is obtained by differentiation

Force  :=  -(∂_r) LJModified
Force

5.

Define new "normalized" variables rNorm = r / rMin and LJNorm = LJ / EMin and use them to define a normalized LJ function LJNorm(rNorm). Recongnizing that I will eventually want to plot these functions, I define them with delayed assignment :=

LJNorm :=  
      Simplify[-(LJModified/EMin) /. rMin -> r/rNorm, 
        Assumptions -> {a > 0  && b > 0}]
LJNorm

ForceNorm := 
      Simplify[-(Force rMin/EMin) /. rMin -> r/rNorm, 
        Assumptions -> {a > 0  && b > 0}]
ForceNorm

To demonstrate that LJNorm and ForceNorm are dimensionless:

EMinUnits = Mass Length^2/Time^2;
ForceUnits = 
        Mass Length^2/(Time^2 Length);
rUnits  = Length;
rMinUnits  =  Length;

rBarUnits  =  rUnits/rMinUnits;
rBarUnits  == 1

FBarUnits = 
        ForceUnits rMinUnits/EMinUnits;
FBarUnits  == 1

So, the variables \(TraditionalForm\`r\&_\) and \(TraditionalForm\`F\&_\) are dimensionless.

6.

Now for the plot. Note in the use of the DisplayFunction option to control whether or not the graphics are rendered to the screen.

NewPlot  =  
      Plot[{LJNorm, ForceNorm},  {rNorm,  0,  2},  
        PlotRange  ->  {-4, 
            4},    PlotStyle  ->  {{Thickness[ .01],  
              Hue[0.3]}, {Thickness[ .01],  Hue[1]}}, 
        AxesLabel -> {"r/rMin", "Normalized LJ and F"}, 
        DisplayFunction -> Identity]

Show[NewPlot, Graphics[Text["Force", {1.5 , -2}]], 
      Graphics[Text["Potential", { .6, 2}]], 
      DisplayFunction -> $DisplayFunction]

7.

I will use the same approach as above, using as my constants the equilibrium position, xMin, and minimum potential energy,

PotEnergy  =  A  + B x  +  k x^2/2

SpringForceEq =  (D[PotEnergy, x] /. x ->  xMin ) ==  
        0

MinEnergyEq =   ( 
          PotEnergy /. x -> xMin) ==  -PotMin

SolutionPot =  
      Solve[{SpringForceEq, MinEnergyEq}, {A, B}]

ModifiedPotEnergy  =  
      PotEnergy  /.  SolutionPot[[1]]

8.

Use function Series to get the Taylor series representation, include Normal to eliminate the O[r - rMin]^3 term.

HarmonicApprox = 
      Normal[Series[LJModified, {r, rMin, 2}]]

Note that this expression is a parabolic fit to the bottom of the potential well for the LJ function.

9.

The coefficient of the second-order term of HarmonicApprox is equal to k/2, where k is the spring constant. We can use the relation between vibrational frequency and spring constant to get the desired expression.

k  =  
      2 SeriesCoefficient[Series[LJModified, {r, rMin, 2}], 2]

VibFrequency  =  (1/(2  π)) (Sqrt[(k/Mass)])