|
<P>clear all<BR>global r<BR>t0=[0 150];%积分时间<BR>y0=[1,0,0];</P>
<P>%bifurcation<BR>for r=20:0.05:30 %r的变化精度<BR>[t,y]=ode45('Lorenz',t0,y0);<BR>[Xmax]=getmax(y);</P>
<P>plot(r,Xmax,'b','markersize',1)<BR>hold on<BR>clear Xmax<BR>end<BR>xlabel('r')<BR>ylabel('Xmax')<BR><BR><BR>function [Xmax] = getmax(y)<BR>a=length(y);<BR>j=1;<BR>for i=(a-1)/2:a<BR> <BR> b=(y(i,1)-y(i-2,1))/2;<BR> c=(y(i,1)+y(i-2,1))/2-y(i-1,1);<BR> <BR> if y(i-2,1)<=y(i-1,1)&y(i-1,1)>=y(i,1)&c==0<BR> Xmax(j)=y(i-1,1);<BR> j=j+1;<BR> elseif y(i-2,1)<=y(i-1,1)&y(i-1,1)>=y(i,1)<BR> Xmax(j)=y(i-1,1)-b^2/(4*c);<BR> j=j+1;<BR> end<BR>end<BR><BR>function dy = Lorenz(t,y)<BR>global r<BR>dy=zeros(3,1);<BR>dy(1)=-10*(y(1)-y(2));<BR>dy(2)=-y(1)*y(3)+r*y(1)-y(2);<BR>dy(3)=y(1)*y(2)-8*y(3)/3;<BR><BR><BR>仅供参考!要想图漂亮一点,可以把积分时间和r的变化精度增大,但是计算时间会变长。</P> |
|