yanghong
2009-03-05, 06:48
有哪位可以用MATLAB帮小弟做这个题:cool: 跪谢:lovely:
Program 1
Npre = 200; Nplot = 300;
y = zeros(Nplot,1);
for c = 0:0.005:2.0,
y(1) = 0.1;
for n = 1:Npre,
y(1) = y(1).^2 - c;
end,
for n = 1:Nplot-1,
y(n+1) = y(n).^2 - c;
end,
plot(c*ones(Nplot,1),y,'.','markersize',2);
hold on;
end,
title('Bifurcation diagram of the quadratic map');
xlabel('c'); ylabel('y_n');
xlim([0 2.0]);
hold off;
Comments to the program:
In the 1st line we specify: Npre, the number of orbit points we want to discard (number of pre-iterates), and Nplot, the number of orbit points we want to plot for each value of parameter c (number of iterates).
Line 2 allocates memory for an array y with Nplot rows and 1 column. The array is initially filled with zeros.
In the 3rd line we start a loop over the values of parameter c, starting from 0 and ending at 2.0 with the increment of 0.005.
Line 4 specifies the starting point for the orbit. The structure of the bifurcation diagram should not depend on a particular choice of the starting point, as long as the value is in the range from −2 to 2. Initial points outside this range result in orbits that diverge to infinity.
The next three lines are the loop where the pre-iterates are computed. Note that, since we do not want to store the values of the orbit points in the pre-iteration loop, we use y(1) on the left-hand-side, so that the old value in y(1) is replaced by the new value of the orbit point.
Next is the loop that computes the points that will be plotted. This loop is identical to the one in Program 1 of Lab 2.
Next line plots the orbit stored in y for a given value of c. Recall that in the Matlab function plot(x,y), the sizes of arrays x and y must be identical. Since c is a 1-by-1 array (a scalar) while y is an Nplot-by-1 array (a vector), we use expression c*ones(Nplot,1) to create a vector whose size is identical to that of y and whose elements are all equal to c. Alternatively, we can use the Matlab function repmat. That is, repmat(c,Nplot,1) produces the same result by replicating the scalar c into Nplot rows and 1 column, thus creating the Nplot-by-1 array.
The plot function also specifies the type of symbols to be used (dots), and contains a pair of arguments ('markersize',2) that change the size of the plotted dots. The default markersize is 6, which is too large and causes many dots to overlap.
We use hold on to retain all the plotted points while plotting orbits for subsequent values of parameter c.
Matlab function xlim is used to specify the lower and upper limits of the horizontal axis. Here it is used to ensure that the limits of the horizontal axis coinside with the range of c values. Similarly, the function ylim can be used to change the limits of the vertical axis. For example, the statement ylim([-2 0]); can be used to zoom in on the lower part of the bifurcation diagram. Try it.
Note that the larger the values of Npre and Nplot the longer the computation time. Sometimes you may need to wait a short while for the diagam to appear. If your program appears to be stuck, press Ctrl-C to interrupt the program execution.
The diagram produced by Program 1 combines information about the asymptotic behaviour of the quadratic map for different parameter values in the range from 0 to 2. Remember that each vertical line contains Nplot points. So, if we see only a small number of points appearing on a given vertical line, it means that these points are repeatedly visited by the orbit and, therefore, the orbit is periodic. It is thus clear that for 0 < c < 0.75 the orbits converge to a fixed point, while for 0.75 < c < 1.25 they converge to a period-2 orbit. At c = 1.25 there is another period-doubling bifurcation, where a period-2 orbit is replaced by a period-4 orbit. In fact, a closer look reveals that further increasing c results in a whole cascade of period-doubling bifurcations occuring closer and closer to each other and producing orbits of periods 8, 16, 32, 64, and so on.
Program 1 can be easily modified to create a bifurcation diagram for any range of parameter values with any increment and any number of points plotted for each c. However, specifying the increment directly is not very convenient, since it is difficult to estimate how many different values of c are necessary to obtain a good-looking diagram. Large increment values result in a very sparse diagram, while very small increment values increase computation time and may overload computer memory (if this ever happens, press Ctrl-C to interrupt program execution). It is more convenient to specify the number of different c values to be used in the given range and then tell the program to calculate the necessary increment. To this end,
Assignment 1. Transform Program 1 into a Matlab function bifdiag(Npre, Nplot, cmin, cmax, Nc), where cmin and cmax give the range of c, and Nc is the number of different c values in this range. Hint: Statement bifdiag(200, 300, 0, 2, 401) should produce a diagram identical to that produced by Program 1.
Assignment 2. Use function bifdiag or Program 1 (if you could not complete Assignment 1) to explore the structure of the bifurcation diagram.
a)By adjusting the values of cmin and cmax, produce a diagram of the region where the cascade of period-doubling bifurcations appears to end and new behaviour starts. Determine the value of parameter c at which this happens. Give the value with accuracy of four digits after the decimal point.
b) Find a place in the diagram where a stable period-3 orbit can be observed. Produce a diagram of this region. Determine the range of values of parameter c for which this orbit exists.
c) Find a place in the diagram which looks almost exactly like the whole diagram. For this you need to change the range of both horizontal and vertical axes. Produce a diagram of this region.
3. The Feigenbaum constant. Let us focus again on the cascade of period-doubling bifurcations producing orbits of periods 2, 4, 8, 16, 32, and so on. Let ck denote a bifurcation point where a period-2k orbit is created. For example, c1 = 0.75 and c2 = 1.25. We see from the bifurcation diagram that the intervals between successive period-doubling bifurcations become smaller and smaller. To quantify this decrease, we introduce ratios of the neighbouring intervals as follows:
δk = (ck+1 − ck) / (ck+2 − ck+1) .
It was shown in 1975 by Mitchell Feigenbaum that, in the limit of large k, the ratios δk converge to a single number δ called the Feigenbaum constant. The remarkable thing about this constant that it has the same value for all maps where a period-doubling bifurcation cascade can be observed. Because of this, it is often referred to as a universal constant and enjoys the same level of respect in the dynamical systems theory as π = 3.1415... in geometry or e = 2.718... in calculus.
Assignment 3. The subject of this exercise is to estimate the value of the Feigenbaum constant δ.
a) Use function bifdiag or Program 1 to construct a bifurcation diagram showing orbits of period 4, 8, 16 and 32 with enough detail to allow you to estimate the bifurcation points c3 , c4 and c5 with the accuracy of at least four digits after the decimal point. (To get better estimates of these values, it may be helpful to learn how to use the 'zoom in' feature available in Matlab figures and put a grid on the plot by typing grid on at the Matlab prompt. You may also need to increase the value of Npre in order to get better converging of orbits near the bifurcation points.) Use these values to calculate δ1, δ2 and δ3. Do these numbers appear to converge to a limiting value?
b) A better way to calculate the Feigenbaum constant is to look for the intersection points of the bifurcation diagram with the horizontal line corresponding to the extremum value (either a maximum or a minimum) of the map function. In case of the quadratic map, the minimum occurs at y = 0. One such intersection point occurs at c'1 = 0 and another at c'2 = 1. Determine the values of c'3, c'4, c'5, and c'6 with the accuracy of at least five digits after the decimal point. To do this, learn to use the Matlab function ginput and use format long to view results with higher accuracy. Calculate δ'1, δ'2, δ'3 and δ'4 using the same formula as in part (a). Do these numbers appear to converge to the same value as in (a)?
Program 1
Npre = 200; Nplot = 300;
y = zeros(Nplot,1);
for c = 0:0.005:2.0,
y(1) = 0.1;
for n = 1:Npre,
y(1) = y(1).^2 - c;
end,
for n = 1:Nplot-1,
y(n+1) = y(n).^2 - c;
end,
plot(c*ones(Nplot,1),y,'.','markersize',2);
hold on;
end,
title('Bifurcation diagram of the quadratic map');
xlabel('c'); ylabel('y_n');
xlim([0 2.0]);
hold off;
Comments to the program:
In the 1st line we specify: Npre, the number of orbit points we want to discard (number of pre-iterates), and Nplot, the number of orbit points we want to plot for each value of parameter c (number of iterates).
Line 2 allocates memory for an array y with Nplot rows and 1 column. The array is initially filled with zeros.
In the 3rd line we start a loop over the values of parameter c, starting from 0 and ending at 2.0 with the increment of 0.005.
Line 4 specifies the starting point for the orbit. The structure of the bifurcation diagram should not depend on a particular choice of the starting point, as long as the value is in the range from −2 to 2. Initial points outside this range result in orbits that diverge to infinity.
The next three lines are the loop where the pre-iterates are computed. Note that, since we do not want to store the values of the orbit points in the pre-iteration loop, we use y(1) on the left-hand-side, so that the old value in y(1) is replaced by the new value of the orbit point.
Next is the loop that computes the points that will be plotted. This loop is identical to the one in Program 1 of Lab 2.
Next line plots the orbit stored in y for a given value of c. Recall that in the Matlab function plot(x,y), the sizes of arrays x and y must be identical. Since c is a 1-by-1 array (a scalar) while y is an Nplot-by-1 array (a vector), we use expression c*ones(Nplot,1) to create a vector whose size is identical to that of y and whose elements are all equal to c. Alternatively, we can use the Matlab function repmat. That is, repmat(c,Nplot,1) produces the same result by replicating the scalar c into Nplot rows and 1 column, thus creating the Nplot-by-1 array.
The plot function also specifies the type of symbols to be used (dots), and contains a pair of arguments ('markersize',2) that change the size of the plotted dots. The default markersize is 6, which is too large and causes many dots to overlap.
We use hold on to retain all the plotted points while plotting orbits for subsequent values of parameter c.
Matlab function xlim is used to specify the lower and upper limits of the horizontal axis. Here it is used to ensure that the limits of the horizontal axis coinside with the range of c values. Similarly, the function ylim can be used to change the limits of the vertical axis. For example, the statement ylim([-2 0]); can be used to zoom in on the lower part of the bifurcation diagram. Try it.
Note that the larger the values of Npre and Nplot the longer the computation time. Sometimes you may need to wait a short while for the diagam to appear. If your program appears to be stuck, press Ctrl-C to interrupt the program execution.
The diagram produced by Program 1 combines information about the asymptotic behaviour of the quadratic map for different parameter values in the range from 0 to 2. Remember that each vertical line contains Nplot points. So, if we see only a small number of points appearing on a given vertical line, it means that these points are repeatedly visited by the orbit and, therefore, the orbit is periodic. It is thus clear that for 0 < c < 0.75 the orbits converge to a fixed point, while for 0.75 < c < 1.25 they converge to a period-2 orbit. At c = 1.25 there is another period-doubling bifurcation, where a period-2 orbit is replaced by a period-4 orbit. In fact, a closer look reveals that further increasing c results in a whole cascade of period-doubling bifurcations occuring closer and closer to each other and producing orbits of periods 8, 16, 32, 64, and so on.
Program 1 can be easily modified to create a bifurcation diagram for any range of parameter values with any increment and any number of points plotted for each c. However, specifying the increment directly is not very convenient, since it is difficult to estimate how many different values of c are necessary to obtain a good-looking diagram. Large increment values result in a very sparse diagram, while very small increment values increase computation time and may overload computer memory (if this ever happens, press Ctrl-C to interrupt program execution). It is more convenient to specify the number of different c values to be used in the given range and then tell the program to calculate the necessary increment. To this end,
Assignment 1. Transform Program 1 into a Matlab function bifdiag(Npre, Nplot, cmin, cmax, Nc), where cmin and cmax give the range of c, and Nc is the number of different c values in this range. Hint: Statement bifdiag(200, 300, 0, 2, 401) should produce a diagram identical to that produced by Program 1.
Assignment 2. Use function bifdiag or Program 1 (if you could not complete Assignment 1) to explore the structure of the bifurcation diagram.
a)By adjusting the values of cmin and cmax, produce a diagram of the region where the cascade of period-doubling bifurcations appears to end and new behaviour starts. Determine the value of parameter c at which this happens. Give the value with accuracy of four digits after the decimal point.
b) Find a place in the diagram where a stable period-3 orbit can be observed. Produce a diagram of this region. Determine the range of values of parameter c for which this orbit exists.
c) Find a place in the diagram which looks almost exactly like the whole diagram. For this you need to change the range of both horizontal and vertical axes. Produce a diagram of this region.
3. The Feigenbaum constant. Let us focus again on the cascade of period-doubling bifurcations producing orbits of periods 2, 4, 8, 16, 32, and so on. Let ck denote a bifurcation point where a period-2k orbit is created. For example, c1 = 0.75 and c2 = 1.25. We see from the bifurcation diagram that the intervals between successive period-doubling bifurcations become smaller and smaller. To quantify this decrease, we introduce ratios of the neighbouring intervals as follows:
δk = (ck+1 − ck) / (ck+2 − ck+1) .
It was shown in 1975 by Mitchell Feigenbaum that, in the limit of large k, the ratios δk converge to a single number δ called the Feigenbaum constant. The remarkable thing about this constant that it has the same value for all maps where a period-doubling bifurcation cascade can be observed. Because of this, it is often referred to as a universal constant and enjoys the same level of respect in the dynamical systems theory as π = 3.1415... in geometry or e = 2.718... in calculus.
Assignment 3. The subject of this exercise is to estimate the value of the Feigenbaum constant δ.
a) Use function bifdiag or Program 1 to construct a bifurcation diagram showing orbits of period 4, 8, 16 and 32 with enough detail to allow you to estimate the bifurcation points c3 , c4 and c5 with the accuracy of at least four digits after the decimal point. (To get better estimates of these values, it may be helpful to learn how to use the 'zoom in' feature available in Matlab figures and put a grid on the plot by typing grid on at the Matlab prompt. You may also need to increase the value of Npre in order to get better converging of orbits near the bifurcation points.) Use these values to calculate δ1, δ2 and δ3. Do these numbers appear to converge to a limiting value?
b) A better way to calculate the Feigenbaum constant is to look for the intersection points of the bifurcation diagram with the horizontal line corresponding to the extremum value (either a maximum or a minimum) of the map function. In case of the quadratic map, the minimum occurs at y = 0. One such intersection point occurs at c'1 = 0 and another at c'2 = 1. Determine the values of c'3, c'4, c'5, and c'6 with the accuracy of at least five digits after the decimal point. To do this, learn to use the Matlab function ginput and use format long to view results with higher accuracy. Calculate δ'1, δ'2, δ'3 and δ'4 using the same formula as in part (a). Do these numbers appear to converge to the same value as in (a)?