|
马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。
您需要 登录 才可以下载或查看,没有账号?我要加入
x
小弟目前正在画一分岔图,关注论坛也有很长一段时间,根据各位大侠的指点,我也进行了程序的编写,无奈就是没有理想的结果出现,还望各位大侠指点一二,先谢过了
微分方程:d2y/dt2+1/RC*dy/dt+1/LC*y=s/LC*Vin
其中R,L,C是已知量
当y>11.3时,s=0;当y<=11.3时,s=1
以下是我编写的程序:
function buck_bifur_Vin_getmax
clear all
t0 = [0 10];
for Vin=linspace(20,40,10)
[t,y] = ode113('buck',t0,[0;0;22;0.02;0.000047;Vin;1;11.3]);
ymax = getmax(y(:,1));
plot(Vin,ymax,'b','markersize',1)
hold on
clear ymax
end
function [ymax] = getmax(y)
a = length(y);
j = 1;
for i=(a-1)/2:a
b=(y(i,1)-y(i-2,1))/2;
c=(y(i,1)+y(i-2,1))/2-y(i-1,1);
if (y(i-2,1)<=y(i-1,1) && y(i-1,1)>=y(i,1) && c==0)
ymax(j) = y(i-1,1);
j=j+1;
if(y(i-2,1)<=y(i-1,1) && y(i-1,1)>=y(i,1))
ymax(j) = y(i-1,1)-b^2/(4*c);
j=j+1;
end
end
end
function dy = buck(t,y)
if y(1)>y(8) y(7)=0;
else y(7)=1;
end
dy = zeros(8,1);
dy(1) = y(2);
dy(2) = -y(2)/(y(3)*y(5)) - y(1)/(y(4)*y(5)) + y(6)*y(7)/(y(4)*y(5));
dy(3) = 0;
dy(4) = 0;
dy(5) = 0;
dy(6) = 0;
dy(7) = 0;
dy(8) = 0;
[ 本帖最后由 咕噜噜 于 2009-9-27 16:27 编辑 ] |
|