Mathematica Lesson III

Physics 80301  Fall 2005

Functions in Mathematica

Notice the use of the underscore on the left hand side of this function definition. The semicolen is
used to indicate "delayed" execution; that is to say, the function is evaluated only when it is called.

g[x_]  :=  (E^Sin)[x]

The function is treated in the same way as a built in function, Thus, for example

g[ArcSin[u]]

or, as a second exmple,

Simplify[D[g[y], {y, 2}]]

or, as a final example,

Plot[g[t], {t, -2, 2}, PlotStyle -> Thickness[0.01], 
      AxesLabel -> {"t", "g[t]"}]

Loops in Mathematica

Do loops

Let us look at a "do-loop" in Mathematica
Here is a loop that produces a set of plots suitable for animation. Select the bracket that contains
the compete set then choose "Animate Selected Graphics from the Cell drop-down menu.

Do[
        Plot[Sin[n x], {x, 0, 2  π}, 
          PlotStyle -> RGBColor[n/3, 0, 0]], {n, 1, 3, 
          0.25}];

Here is a summation loop 
Notice the C-type notation for summation (s += 1/2^i  is equivalent to s = s + 1/2^i).

Clear[s] ; s = 1 ; Do[
        s  +=  1/2^i, {i, 
          40}] ;

N[s, 20]

For loops

Here are a few examples of loops usimg For[start,test,incr,body]

For[i = 1, i < 5, i++, Print[1/i]]

Here is a more complicated For-loop (an you guess what it does?)

For[i = 1;  t =  1/x,  i^2 < 10, i++, t = t/x; 
      Print["Integral(1--infinity)[1/x^", i + 1, "]=", 
        Integrate[t, {x, 1, Infinity}]]]

While loops

While loops have the structure While[test, body]. As an example,

Clear[nn]; nn = 10; While[(nn = nn - 1) > 5, Print[nn]]

FixedPoint and NestWhile (see help browser)

Nonlinear Curve Fitting

In this section, we consider the problem of fitting experimental data to a nonlinear function. We have in mind the problem of fitting the intensity profile of a laser beam. Such a problem was assigned this semester in the Laser Physics course. The form of the intensity profile is assumed to be F(x) = a(e^(-(((x - c))^2)/(b^2))). We want to choose parameters (a,b,c) to fit experimental data to this functional form using the nonlinear fitting routine from the "Statistics" package. The package is loaded by the command:

<<  Statistics`NonlinearFit`

First, let's construct some psuedo data for the intensity:

intens  =  
        Table[(4 + 
                0.4 (Random[] - 0.5)) Exp[-(2  + 
                      0.2 (Random[] - 0.5)) (x - (2  + 
                          0.2 (Random[] - 0.5)))^2], {x, 0, 4, 
            0.05}];

Save the distance coordinates used to generate the intensity data

dist = Table[x, {x, 0, 4, 0.05}];

Make up a data table with two columns [dist, intens]

data  =  Transpose[{dist, intens}];

Plot the data table, but hold back the display:

plot1  = 
        ListPlot[data, 
          PlotStyle -> {PointSize[0.02], RGBColor[1, 0, 0]}, 
          AxesLabel -> {"ρ", "I(ρ)"}, 
          DisplayFunction  -> Identity];

Do the nonlinear fit:

f  = 
      NonlinearFit[data, a Exp[-(x - c)^2/b^2], 
        x, {a, b, c}]

Examine the fit parameters (a,b,c) :. [We expect approximately (a=4,b=1/(Sqrt[2]),c=2)]

NonlinearRegress[data, a Exp[-(x - c)^2/b^2], 
      x, {a, b, c}, 
RegressionReport -> 
        ParameterCITable]

Plot of fit function and hold it back the display.

plot2 = 
        Plot[f, {x, 0, 4}, PlotStyle -> Thickness[0.01], 
          DisplayFunction -> Identity];

Now display the two graphs together:

Modules in Mathematica

Here we give an example of a module. This is Mathematica's version of a Maple procedure or a
Fortran subroution. In this example, we solve a Laplace equation boundary-value problem in 2 dimensions and plot the results.

a)  In this example "imax" designates number of x intervals, "jmax" designates number of y intervals
amd "maxit" designates the maximum number of iterations permitted. These three quantities are input parameters. 

b) The terms {v,v0,del=1./jmax,err} are local variables, known only inside the module. 

c) The statements enclosed in (* .... *) are comments and disregarded in the execution block

d)  Notice the "do loop" used to initialize the arraya v and v0:
                  Do[v[i,j]=0;v0[i,j]=0,{i,0,imax},{j,0,jmax}];
                  
e)  Notice the "for loop" used in the iteration: the notation   m++  has the same meaning as in C               
  
f)  Notice how Show[] is used to display the table as a surface.

\(laplace[imax_, jmax_, maxit_]\  := \ 
      Module[{v, v0, del = 1. /jmax, err}, 
 (*\ 
          clear\ arrays\ v\ and\ v0\ to\ start\ *) 
\ \ \ \ \ \ Do[v[i, j] = 0; 
          v0[i, j] = 0, {i, 0, imax}, {j, 0, jmax}]; 
 (*\ 
          impose\ boundary\ condition\ V = \(1\ at\ x = \(0

                \ --\)\ all\ other\ boundary\ values\ are\ 0\)\ \ *) 
\ \ \ \ \ \ \ Do[v[0, j] = 1; 
          v0[0, j] = 1, {j, 0, jmax}]; 
 (*\ 
          Set\ up\ an\ iteration\ loop\ *) 
For[
          m = 1, 
\ 
          m < maxit\ , 
\(m++\), 
\ \ \ \ \ Do[Do[v[i, j]\  = \ 
                0.25\ \((v0[i - 1, j] + v0[i + 1, j] + v0[i, j - 1] + 
                      v0[i, j + 1])\), {i, 1, imax - 1}], {j, 1, 
              jmax - 1}]; 
\ \ \ \ \ err\  = \ 
            Max[Table[
                Abs[v[i, j] - v0[i, j]], {i, 1, imax - 1}, {j, 1, 
                  jmax - 1}]]; 
\ \ \ \ \ \ Do[
            Do[v0[i, j]\  = \ v[i, j], {i, 1, imax - 1}], {j, 1, 
              jmax - 1}];\ 
]; 
\ \ Print[" Number of iterations = ", maxit, "  Max Error = ", 
          err]; 
Show[
          SurfaceGraphics[
Table[
              v[i, j], {j, 0, jmax}, {i, 0, 
                imax}], 
MeshRange -> {{0, 
                  jmax + del}, {0, 
                  imax + del}}], 
BoxRatios\  -> \ {1, 
              1, 1}, \ PlotRange -> {0, 
              1}, 
Boxed -> False, \ 
          Axes -> 
            True, 
AxesLabel -> {"x", "y", "V(x,y)"}, 
TextStyle -> {FontFamily\  -> \ "Times-Roman", \ 
              FontSize\  -> \ 
                11}
]
]\)

Now, let's try it out:

laplace[20, 24, 80]