Lab 6: Physical Simulation II

At each location (x,y), the second derivative of the height of the water is equivalent to the Laplacian of the height field.

or, equivalently,

If you think of bit of water, it rises and falls with an acceleration equivalent to the Laplacian of the height field.

Last week we simulated the trajectory of a single particle under uniform acceleration, this week we will simulate water by calculating the acceleration everywhere and updating all the bits of water at once. In this lab we will continue the pattern of asking you to work out more and more of the details for implementation. The main tasks you need to accomplish are:

1. Write a function to compute the Laplacian of the height field. This can be approximated by the difference between the average of its 4 neighbors and the current value of the height-field at that location. You will have to make some decision about how to deal with the boundary (again!).
   1. Test this with an initial 5x5 height field array with zeros everyone and one location (in the middle) with value 1. After one iteration, you should get a negative value at the center and positive values immediately surrounding it.

      >> R 
      R =
           0     0     0     0     0
           0     0     0     0     0
           0     0     1     0     0
           0     0     0     0     0
           0     0     0     0     0
      >> laplacian(R)
      ans =
               0         0         0         0         0
               0         0    0.2500         0         0
               0    0.2500   -1.0000    0.2500         0
               0         0    0.2500         0         0
               0         0         0         0         0
      

   Find a TA if you are having trouble getting this.
2. Create a new matlab script (edit simulateWater) and start an initialization portion of the main script which:
   1. Defines a 50x50 array that represents the height of the water at each location.
   2. Defines a 50x50 array that represents the up/down velocity of the water at each location. (could start as all zeros)
   3. defines the constant "c". Later on, you should vary this to see how it affects the simulation, a value of 2 should get you started.
   4. defines the constant "dt", the temporal step size. Later, you should vary this, but 0.1 is a good value for this too.
   5. initializes the water height field. (all zeros with a 1 in the middle)
   6. displays the initial height field. (use imagesc)
3. Create an "main loop" which runs many times and:
   1. Calls an updateWater function (which you must create) which uses the velocity field to update the current height field, and the laplacian to update the velocity. As a hint, your updateWater function should start something like:

      function [water,waterV] = updateWater(water,waterV,c,dt)
      

      This should use the Laplacian of "water" and the constant "c" to compute the acceleration. Then use that acceleration to create the new velocity field.
   2. visualizes the current water height.
4. Now, you should be able to run your simulation, by running the simulateWater script, and see a circular wave radiating out from the center. Find a TA if you are having trouble getting this.
5. Now we are going to simulate "rain" by adding to your main loop. Randomly, about 1 time in 100, you should randomly choose some location (in your 50 x 50) grid, and set the velocity to be negative 1 (as would happen if a raindrop came from above and pushed down that bit of water).
   1. The matlab call rand(1) returns 1 number randomly chosen between zero and 1.
   2. You might use a statement like "if rand(1) < 0.1" to create a loop that would run, on average, 1 time in 10.
   3. You might use a statement like xx = 1+floor(rand(1) * 50), to create a random number between 1 and 50
6. let this simulation run for about 1000 steps and show the resulting image to a TA. This is the only required checkpoint of the day.
7. If the TA is occupied, or if you want to learn another interesting visualization tool in matlab, replace you "imagesc(water)" line with: h = surf(water); axis([0 50 0 50 -0.5 0.5])
8. You might notice that the water never seems to "calm down" after a bunch of drops have arrived. You might see how you could alter your udpateWater function to add a little bit of "friction". One way to do this would be just to "slow down" the water, with an additional term such as "waterV = waterV .* 0.99", which very slight slows down the water each step