Pretty Graph of the Day, Part II

solutionSnap

I like this series of graphs because it looks like someone is pulling on the solution like it's a string, until SNAP!  What's being shown here is the numerical solution to the initial value problem

y'=\frac{1}{(1-y)^{2}}, y(0)=0, t\in[0,1].

As you might suspect from the denominator, there's trouble when y=1.  This results in a problem that's not well-posed.  Using the same numerical scheme but only a different number of grid points produces wildly different behaviors.  What's basically happening is that if the numerical scheme hits the trouble region "just right" due to the grid point spacing, the gigantic derivative sends the solution flying off to some large and incorrect value.

You can reproduce this figure in Mathematica with the following code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
(* 4th Order Runge-Kutta Method *)
rungeKutta[y0_, npoints_, a_, b_, f_] := (
   h = (b - a)/(npoints - 1);
   y = Table[0, {npoints}];
   t = Table[0, {npoints}];
   y[[1]] = y0;
   t[[1]] = a;
   Do[
    yn = y[[n]];
    tn = t[[n]];
    m1 = f[tn, yn];
    m2 = f[tn + 0.5*h, yn + 0.5*h*m1];
    m3 = f[tn + 0.5*h, yn + 0.5*h*m2];
    m4 = f[tn + h, yn + h*m3];
    y[[n + 1]] = yn + (1/6)*h*(m1 + 2*m2 + 2*m3 + m4);
    t[[n + 1]] = tn + h;
    , {n, 1, npoints - 1}];
   solution = Table[{t[[i]], y[[i]]}, {i, 1, Length[t]}];
   Return[solution];
   );
 
 f[t_, y_] := 1/(y - 1)^2
y0 = 0;
a = 0; (* Starting t *)
b = 1; (* Ending t *)
 
rks = {};
Do[
  npoints = 10^5 + i; (* Number of grid points *) 
  sol = rungeKutta[y0, npoints, a, b, f];
  p = ListPlot[sol, PlotStyle -> {PointSize[0.001] , Orange}, 
    Frame -> True, AxesLabel -> {"t", "y"}, 
    PlotLabel -> 
     "10^5+ " <> ToString[i] <> 
      " grid points"];
  AppendTo[rks, p];
  , {i, 0, 4}];
 
GraphicsColumn[Table[rks[[j]], {j, 5}], Frame -> All, ImageSize -> 300]

It'll probably take several seconds to run because of the large number of grid points used in each solution.

Fun with Diffraction Gratings

A laser beam passing through a transmission diffraction grating straight on gives the standard diffraction pattern we all know and love.  It's a bit more interesting, though, if the beam hits the grating at an angle:

grating1

figure 1

Most intro optics books cover this situation, and the result is (equation 1):

a[\sin(\theta_{m}) - \sin(\theta_{i})] = m\lambda

where a is the grating spacing, \theta_{m} is the angle of the mth maxima, \theta_{i} is the incident angle, m is the maxima order, and \lambda is the wavelength. We can rewrite this (homework) in a more useful way as (equation 2):

\theta_{m}=\arcsin[\frac{m\lambda}{a} + \sin(\theta_{i})]

A similar, but slightly more complicated situation happens when you rotate the grating instead of the lazer:grating

figure 2

With a rotated grating (figure 2), the laser is still hitting the grating at an angle as it is in figure 1.  So, starting with equation 2 and using some geometry we get the angles that satisfy the maxima condition:

\theta_{m'}=\arcsin(\frac{m\lambda}{a}) + \sin(\theta_{g}) - \theta_{g}

This seemed liked a fun thing to model in Mathematica, especially since I had never played with any of the graphics features before.



You can view the source here.

You need Wolfram's CDF player, but it's totally worth because then you can look through all of the awesome demonstrations they have online.