DAY24
Thursday, 05 July 2012

Today is a LAB day.

We have a LAB #7 to be done during class time, for credit. When you have completed the lab, please let Detelina know so that she can approve and record your work.

LAB#7:

I have shown you some simplified random number generators, but it is time to become familiar with the built-in function, called rand().

Task 1: Write a program that calls rand() 10 times, printing out the value each time. The rand() function returns an integer, so use the proper format. The heart of your program might look like this:

        loop on j
          i = random int between 0 and RAND_MAX
          print i
      

Task 2: Inside the same program, and just before the loop you've already written, print out the value of RAND_MAX. RAND_MAX is a value defined by C in the include file <stdlib.h>, so your program must now add the appropriate include statement. RAND_MAX is an integer, and its value is the largest integer that will be returned by rand().

        print RAND_MAX
        loop on j
          i = random int between 0 and RAND_MAX
          print i
      

Task 3: Ordinarily, it's fine if we get different random numbers each time we run the program. But for debugging, or comparing results, it is useful to start the random sequence at the same place. This is done with the srand() function. We input a "seed" value, and that specifies how the random numbers will be determined. Let's all use a seed of 12345.

        srand ( 12345 );
        print RAND_MAX
        loop on j
          i = random int between 0 and RAND_MAX
          print i
      

Task 4: If we want random real numbers between 0 and 1, then one way to get them is to call rand(), divide the result by RAND_MAX, and store that in a real variable. Add a new set of statements to your program that calls rand() 10 times, and using this idea to compute a real number r that you have declared as a float. Each time, print the value of r.

        srand ( 12345 );
        print RAND_MAX
        loop on j
          r = random real number (float) between 0 and 1
          print r
      

Task 5: If we want a random real number s between a and b, we can compute a random real number r between 0 and 1, and then compute

        s = a + r * ( b - a )
      

        set a
        set b
        srand ( 12345 );
        print RAND_MAX
        loop on j
          r = random real number (float) between 0 and 1
          s = random real number between a and b
          print s
      

Task 6: We are finally ready to do something interesting. We are going to simulate a random walk in one dimension. The idea is that someone starts at the position x=0, and takes a number of steps, each of a random size. Our question is, if you take n random steps, how far away from your starting point are you likely to be?

Now we are going to try this experiment 1,000 times. On each try, we are going to take 100 random steps. Each step will be a random number between -1 and +1. Our position is the sum of these 100 steps. Theoretically, we could end up 100 units away.

      set a
      set b
      srand ( 12345 );
      print RAND_MAX

      loop on i (number of experiments):

        Initialize x = 0;

        loop on j
          r = random number between 0 and 1
          s = random number between -1 and +1
          x = x + s.

        print x     <-- This is after the J loop is completed!
      

Run your program for 1,000 experiments (I), with each experiment taking 100 steps (J). Your program should print out a list of 1,000 numbers. You will probably believe that the average value is 0, but is there a pattern to the values? Save the output of your program to a file, so that we can feed it to gnuplot:

        ./a.out > data.txt
      

Task 7: The data is not sorted or binned in any way. We could have our program do this for us. However, we can also ask gnuplot to automatically bin our data with the following pair of somewhat crazy commands:

      gnuplot
        bin(x) = 0.5 * floor(x/0.5) + 0.25
        plot "data.txt" using (bin($1)):(1.0) smooth freq with boxes
      

HOMEWORK #7 (must be turned in by next Thursday):

Homework 7.1:

If y is a function of x, it is sometime possible to describe the behavior of y, as x varies, by giving an initial value y(x0)=y0, and a formula for the derivative of y with respect to x:

        y(x0) = y0
        dy/dx = f(x,y)
      

As we have discussed in class, we can think about approximating the exact solution y(x) by computing estimates for y at equally spaced points x0, x0+dx, x0+2*dx, ... We start by approximating the derivative by a difference:

        ( ynew - y ) / dx = f (x,y)
      

        ynew = y + dx * f(x,y)
      

        x[0] = x0     (Initialization)
        y[0] = y0

        x[1] = x[0] + dx     (Step 1)
        y[1] = y[0] + dx * f(x[0],y[0])

        x[2] = x[1] + dx     (Step 2)
        y[2] = y[1] + dx * f(x[1],y[1])
      

Write a program that can solve differential equations using this idea. Your program should assume that the right hand side is computed in a C function whose definition is

        float f ( float x, float y )
      

        x0 = 0
        y0 = 3
        dy/dx = 1 - 3 * sin(x)
      

Print the values of x and y, and then use gnuplot to make a plot of your solution.

Homework 7.2:

The secant method can be used to try to estimate a value x for which f(x) is zero. Such a value is called a root of the equation f(x)=0. The method works by starting with two estimates for the root, which we will call a and b; it then computes fa=f(a) and fb=f(b). An improved estimate of the root should be c=b-f(b)*(b-a)/(fb-fa).

Once we have the estimate c, we want to compute fc = f(c). If this number is small (in absolute value!) then we will say it is close enough to the root. Otherwise, we want to try again. To try again, we need to update the old information:

        a = b
        fa = fb
        b = c
        fb = fc
      

The outline of your code should be something like this:

        initialize a, b, fa, fb;
        loop:
          compute c
          evaluate fc
          print c and fc
          if |fc| is small enough, then leave the loop
          update a, fa, b, fb, and prepare to loop again
      

Write a program that carries out the secant method. Your function f(x) should be defined using a C function. The information to use is f(x) = sin(x)-x/2, a = -4, b = 2. Leave the loop if you find that |f(c)| <= 0.000001.

(This loop has the potential to run forever, so remember that you can always kill your program using Control-C (hold down the Ctrl key and hit C). You could also protect yourself by stopping if the loop takes more than 20 steps.)

Homework 7.3: (Graduate students only!) Consider the mathematical function f(x) = 1/(1+x^2). We want to approximate the integral of this function over the interval [0,1]. We will use the Monte Carlo method to do this.

The Monte Carlo method needs a sequence of random numbers in the interval [0,1]. We can get these by calling the rand() function (which returns random nonnegative integers) and dividing by RAND_MAX, a symbolic name for the largest integer returned by rand(). In other words, to compute 100 random numbers, you would do something like this:

        for ( i = 0; i < 100; i++ )
        {
          r = ( double ) rand ( ) / ( double ) RAND_MAX;
        }
      

To use the Monte Carlo method to approximate the integral, we want to start with a loop like the one above, but add statements that evaluate the function f() at each random value r, sum all these values, and then, (after the loop!) take the average. For this particular simple problem, that's all we need to do to estimate the integral. The exact value of the integral is Pi/4, or roughly 0.785398.

Write a program that uses the Monte Carlo approach to estimate the integral of f(x)=1/(1+x^2), using 10,000 random numbers. Print your estimate of the integral, and print the error, that is, the difference between your estimate and the exact result.

For each of your homework programs, turn in a copy of the program, and a copy of the output.

Programs we might discuss: