下面的程序应该可以参考一下
- function feigenbaum
- % feigenbaum.m
- % Jeremie Korta
- % December 1, 2003
- % The following program will compute the Feigenbaum ratio using parameter
- % values r(n) at superstable fixed points of f^(2^n). As a check, the
- % program can be used to estimate the same ratio using period doubling
- % cascade from 3-cycle, which is superstable at r = 1.7549... taking the
- % unimodular map to be r-x^2.
- clear n r R c;
- global n q c
- max=10;
- q = 2; % if you want to find the Feigenbaum ratio for higher order maxima
- % change this parameter (the 3-cycle will no longer be stable at the same r!)
- % Known superstable r for 1 and 2-cycles for any q.
- c = 0; % signifies we're searching from 2-cycle
- r(1)=0; % for 1-cycle to be superstable, r = 0
- r(2)=1; % for 2-cycle, r = 1.
- % We can also calculate the Feigenbaum ratio at the 3-cycle
- % period-doubling cascade (for q=2). Hopefully, it's still 4.6692... !
- %c = 1; % signifies we're searching from 3-cycle
- %r(1) = 1.7549; % superstable 3-cycle occurs at r = 1.7549...
-
- for i=2:max
- n=i-1;
- del=10^(-n-c); % smaller for 3-cycle window
- next=r(i-1)+del;
- prev_sign=sign(func(next));
- while (sign(func(next))==prev_sign)
- next=next+del;
- end
- r(i)=fzero(@func, [r(i-1)+10^(-n-1-c), next]);
- end
- for j=1:(max-2)
- R(j)=(r(j)-r(j+1))/(r(j+1)-r(j+2));
- end
- fprintf('\n Best estimate of Feigenbaum ratio = %1.6f', R(max-2))
- % File: func.m
- %
- % func returns the f^(2^n) map of 0 where f = r-x^(2q) q=1,2,3... if
- % c = 0, or the f^(3*2^n) map of 0 if c = 1. The zero of this function
- % occurs at the superstable value r(n).
- function x=func(r)
-
- global n q c
- x=0;
- for i=1:((c > 0)*2+1)*2^n
- x = r - x^q;
- end
复制代码 |