马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。
您需要 登录 才可以下载或查看,没有账号?我要加入
x
<P>本是研究一个非线性常微分方程组的,因为带有控制变量,不知道怎么写程序,所以在写matlab程序时就把它当参数看待,也不知道是不是可以。这样就类似带参数r的洛伦茨方程了,但是针对r的某个具体值,在积分区间内变量x,y,z是有很多值的,且我发现变量长度不等,所以如果以r为横轴,某个状态变量为纵轴时就遇到了麻烦,不知道状态取哪个值,我在论坛看到一个程序,如下面,但我不是很明白函数getmax的意思,谁能给解释一下呢,不知道是不是当r为具体某个值时,取状态变量的最大值?但状态量有3个啊,怎么能画在一张图上,不明白,如果我只想看一个状态量随r的变化情况如何写程序。<br><br></P>
<P align=center>clear all<br>global r<br>t0=[0 150];%<FONT face=宋体>积分时间</FONT><br>y0=[1,0,0];
<br>
<p>
<P align=center>%bifurcation<br>for r=20:0.05:30 %r<FONT face=宋体>的变化精度</FONT><br>[t,y]=ode45('Lorenz',t0,y0);<br>[Xmax]=getmax(y);
<p>
<p>
<P align=center>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;
<p>
<p><br>
[此贴子已经被作者于2006-6-18 22:38:12编辑过]
|