From CSE100B Wiki
More Fun with Simulations
Example from class ClassExample
A quick introduction to Conway's Game of Life (GoL) can be found on Wikipedia. The "game" is actually a set of update rules that govern the evolution of a 2D grid.
A cell in a 2D grid has 8 neighbors (including the 4 diagonal ones). In GoL, a cell can be in one of two states - either "alive" or "dead", and the rules determine the state in the next time step (also known as "generation"). At every generation, we look at all the cells, and for each one, we count how many of their neighbors (regardless of position) are alive. The rules are:
- If a cell has exactly 3 living neighbors, it will be alive in the next generation (regardless of its current state).
- If a cell is currently alive, and has exactly 2 living neighbors, it will remain alive in the next generation.
- All other configurations (4 neighbors or more, or one or zero neighbors) result in a dead cell.
lab9_1.m: Implement the Game of Life.
To get started create lab9_1.m with the following code:
grid_size = [10, 10]; simulation_steps = 100; grid = round(rand(grid_size)); %initialize the state grid for t=1:simulation_steps grid = game_of_life(grid); spy(grid); pause(0.25); end
Your job is to write the game_of_life function that creates a new generation from the current generation according to the rules above. Be sure to not modify the current state until you have computed the new state. The spy function shows which elements of a matrix are not zero.
If we have the matrix X = [1 2 3; 4 5 6; 7 8 9] then the call get_val (2, 2, X) returns 5. The call get_val (4, 4, X) returns 1. You can use this function to return the value of the 9 neighbors of each cell in the grid without worrying if the cell is on an edge or corner. The code for this function is:
function value = get_val(row,col,grid) [num_row num_col] = size(grid); x = mod(row-1,num_row)+1; y = mod(col-1,num_col)+1; value = grid(x,y); end
Next, we'll simulate a diffusion process in 2D. Diffusion makes nearby things more similar to each other. One common example is heat transfer where heat transfers from high temperature to low temprature. This means that an update step of a diffusion process should make the value of a cell more similar to the average of its neighbors. The following equation models a diffusion process:
grid_new(r,c) = a*grid(r,c) + (1-a) * (grid(r-1,c) + grid(r+1,c) + grid(r,c-1) + grid(r,c+1)) / 4;
The parameter a determines how fast grid_new(r,c) assimilates the value of its 4 immediate neighbors. The value of a should be between 0 and 1. Basically, this equation takes a weighted average of the cell's current value and the average of its 4 immediate neighbors to compute the new value of the cell.
lab9_2.m: Implement a diffusion process. You are free to choose the boundary conditions: 1) a wrap-around environment, 2) assume that the edge of our lattice has neighbors whose value is equal to their own, or 3) fixed to some constant like 0. Start with a random matrix and watch the values smooth out. You can plot the results using imagesc(grid, [0 1]); pause (.1);. It is tempting to write the update rule for the single lattice site and then just iterate over the entire grid. However, a MUCH faster version can be implemented using matrix notation!
Assume we have a matrix X = [1 2 3; 4 5 6; 7 8 9]. If we want to compute the bottom neighbor of every element (with wrap-around), we use bottom = [X(2:end,:); X(1,:)]. If we want to compute the left neighbor of every element (with 0s on the edge), we use left = [zeros(size(X,1),1) X(:,1:end-1)]. Use these ideas when implementing your diffusion process.
Wrap-around can also be realized using the function circshift or padarray (for example padarray(grid, [1 1], 'circular');.
So far we've dealt with vectors which take 1 index and matrices which take 2 indices. Multidimensional arrays are a simple extension of the regular matrix syntax and take 3 or more indices. A multidimensional array that takes 3 indices can be visualized as a cube. One practical example of such a structure is a color image. As you might know, the space of colors is 3-dimensional (often spanned by combinations of Red, Green and Blue, or RGB) Each pixel in a color image requires three numbers to be properly defined. We can store the three channels in the three "pages" of our data cube - the first page describes the matrix of the red values (so C(i,j,1) is the red value of the pixel at row i and column j), the second page - C(i,j,2) - holds the green channel, and the blue channel is stored in the third page of our cube. MATLAB will draw such an image with image(C).
lab9_3.m: Create a "movie" of three diffusion processes overlapping. Create a 100 by 100 by 3 random matrix and allow each channel (matrix page) to diffuse independently. By setting different rates for the different processes, we will get colorful patches.
Demonstrate lab9_1.m, lab9_2.m and lab9_3.m to the TA or instructor.
As you can see, in a diffusion system, the chemicals (simulated as colors) simply spread out and blend into each other, resulting in gray as the final color. In a reaction-diffusion system - a process in which two or more chemicals not only diffuse over a surface but also react with one another - the end result can be a stable patterns if the right formulas are used to simulate the reaction among the chemicals. Read more about simulating Reaction-Diffusion systems here if you're interested in using this idea for your course project.
Example of a simulation of Belosouv-Zhabotinsky reaction. Read more on Wiki.
Student Questions and Answers
Questions? Please post them on the Google Group http://groups.google.com/group/wustl-matlab.