Problem with a lab

I have a lab problem. It says to write a program to iteratively generate points in two-dimensional space using the following rules:

    (x_{n+1},y_{n+1}) = 

/ (0.5,0.27*y_n),                                                  with 2% probability;
|                
| (-0.139*x_n + 0.263*y_n + 0.57, 0.246*x_n + 0.224*y_n - 0.036),  with 15% probability;
<
| (0.17*x_n - 0.215*y_n + 0.408, 0.222*x_n + 0.176*y_n + 0.0893),  with 13% probability;
|
\ (0.781*x_n + 0.034*y_n + 0.1075, -0.032*x_n + 0.739*y_n + 0.27), with 70% probability.

[Click for pdf version of this problem, where you may find the equation easier to read.]

Start from an initial point (x_1,y_1)=(0.5,0.0). Carry out the iteration at least 30,000 times and plot all the data you obtain (as points) in an x-y plot.

I have written:

for iter = 1: 30000
    r = rand();
    x = 0.5;
    y = 0;
    if r < 0.02
       y = y * 0.27*r;
    elseif r < 0.02 + 0.15
        x = x -0.139 * r + 0.246 * r;
        y = y + 0.263 * r + 0.224 * r - 0.036;
    elseif r < 0.02 + 0.15 + 0.13
        x = x + 0.17 * r - 0.032 * r;
        y = y + 0.034 * r + 0.739 * r + 0.27;
    elseif r < 0.02 + 0.15 + 0.13 + 0.70
        x = x + 0.781 * r - 0.032 * r;
        y = y + 0.034 * r + 0.1075 + 0.739 * r + 0.27;
    end
end

...doesn't seem to work; only outputs one x, y, and r value without repeating 30000 times. What's the issue?

Answer 1

[James Tursa]

The intent of the assignment is to create vectors x and y, not just single values of x and y as your code is currently doing. So you need to use indexing on those x and y values in your loop. The starting point will be x(1) and y(1). The next point will be x(2) and y(2) which will depend on x(1) and y(1) per the random formula. In general, each successive x(n+1) and y(n+1) will depend on x(n) and y(n). E.g., an outline:

N = 30000;
x = zeros(N,1);  % <-- x is a vector, allocate N elements
y = zeros(N,1);  % <-- y is a vector, allocate N elements
x(1) = 0.5;  % <-- x starting point
y(1) = 0;    % <-- y starting point
for n = 1:N-1
    r = rand();
    if r < 0.02
        x(n+1) = 0.5;          % <-- the formula for x(n+1) in terms of x(n) and y(n)
        y(n+1) = 0.27 * y(n);  % <-- the formula for y(n+1) in terms of x(n) and y(n)
    elseif r < 0.02 + 0.15
        x(n+1) = x -0.139 * r + 0.246 * r;           % <-- FIX THIS TO BE IN TERMS OF x(n) AND y(n)
        y(n+1) = y + 0.263 * r + 0.224 * r - 0.036;  % <-- FIX THIS TO BE IN TERMS OF x(n) AND y(n)
    elseif r < 0.02 + 0.15 + 0.13
        x(n+1) = x + 0.17 * r - 0.032 * r;  % <-- FIX THIS TO BE IN TERMS OF x(n) AND y(n)
        y(n+1) = y + 0.034 * r + 0.739 * r + 0.27;  % <-- FIX THIS TO BE IN TERMS OF x(n) AND y(n)
    else
        x(n+1) = x + 0.781 * r - 0.032 * r;  % <-- FIX THIS TO BE IN TERMS OF x(n) AND y(n)
        y(n+1) = y + 0.034 * r + 0.1075 + 0.739 * r + 0.27;  % <-- FIX THIS TO BE IN TERMS OF x(n) AND y(n)
    end
end

I changed the loop counter to match the formula above so it would be easier for you to see the correlation between the code and the formula. I fixed up the 2% case. You need to fix up all of the other cases.

[Gaëtan Poirier]

Thank you very much. This helped a lot