画分岔图，出现Vectors must be the same lengths.

lxx244lxx

function tanqingweiqi
hold on;box on;%生成空的坐标轴
N=100;%设置迭代次数
Tn=[];%初始化
kg2a=5.83;
options=odeset('RelTol',1e-9);%设置微分方程求解选项
for kg1a=5:0.01:8;%计算不同的kg1a值
    x0=[0,0,500];%设置初值
    tspan=[0:0.001:3]
   y=ode4(@ODE317_2,tspan,x0,options,kg1a,kg2a)
    x0=y(end,:);
    for n=1:N;
        y=ode4(@ODE317_2,[0:0.001:0.005],x0,options,kg1a,kg2a);%求解微分方程
        x0=y(end,:);%更新初值
        xd=y(:,3);%提取求解结果
        [M,IK]=max(xd);%计算最大值
        if 1<IK & IK<length(xd);%判断是否为最大值内部点
            Tn=[Tn;M];%记录最大值
        end
    end
    plot(kg1a,Tn,'k.','markersize',1);%绘图
    Tn=[];%清空
   % pause(0.01);
end
end

function dy=ODE317_2(t,y,flag,kg1a,kg2a)
dy=zeros(3,1);
r=0.5;pb=1175;Fco=1400;Fco2=400;F=20000;
r1x=(kg1a*16.84*(1-y(1))*533529*exp(-76857/8.314/y(3))*(16.84*(1-y(1)))^0.674)/(kg1a*16.84*(1-y(1))+533529*exp(-76857/8.314/y(3))*(16.84*(1-y(1)))^0.674);
r2x=(kg2a*4.81*(1-y(2))*296452638*exp(-105858/8.314/y(3))*(4.81*(1-y(2)))^0.008)/(kg2a*4.81*(1-y(2))+296452638*exp(-105858/8.314/y(3))*(4.81*(1-y(2)))^0.008);
dy(1)=(pi*r*r*pb/Fco)*r1x;
dy(2)=(pi*r*r*pb/Fco2)*r2x;
dy(3)=(pi*r*r*pb/F)*(r1x*1000*(188.1276+0.0724*y(3)-4.2782*10^(-5)*y(3).^2+7.9724*10^(-9)*y(3).^3)+r2x*1000*(146.4633+0.06996*y(3)-2.7313*10^(-5)*y(3).^2+1.8928*10^(-9)*y(3).^3))/(1-0.14*y(1)-0.04*y(2))/(26.2932+0.01354*y(3)+6.1860*10^(-6)*y(3).^2-4.1355*10^(-9)*y(3).^3);

qibbxxt

1# lxx244lxx

1. 
   
2. function tanqingweiqi
   
3. hold on;box on;%生成空的坐标轴
   
4. N=100;%设置迭代次数
   
5. Tn=[];%初始化
   
6. kg2a=5.83;
   
7. options=odeset('RelTol',1e-9);%设置微分方程求解选项
   
8. for kg1a=5:0.01:8;%计算不同的kg1a值
   
9.     x0=[0,0,500];%设置初值
   
10.     tspan=0:0.001:3;
    
11.     [t,y]=ode45(@ODE317_2,tspan,x0,options,kg1a,kg2a);
    
12.     x0=y(end,:);
    
13.     for n=1:N;
    
14.         [t,y]=ode45(@ODE317_2,0:0.001:0.005,x0,options,kg1a,kg2a);%求解微分方程
    
15.         x0=y(end,:);%更新初值
    
16.         xd=y(:,3);%提取求解结果
    
17.         [M,IK]=max(xd);%计算最大值
    
18.         if 1<=IK & IK<=length(xd);%判断是否为最大值内部点
    
19.             Tn=[Tn;M];%记录最大值
    
20.         end
    
21.     end
    
22.     plot(kg1a,Tn,'k.','markersize',1);%绘图
    
23.     Tn=[];%清空
    
24.     % pause(0.01);
    
25. end
    
26. 
    
27. function dy=ODE317_2(t,y,kg1a,kg2a)
    
28. dy=zeros(3,1);
    
29. r=0.5;pb=1175;Fco=1400;Fco2=400;F=20000;
    
30. r1x=(kg1a*16.84*(1-y(1))*533529*exp(-76857/8.314/y(3))*(16.84*(1-y(1)))^0.674)/(kg1a*16.84*(1-y(1))+533529*exp(-76857/8.314/y(3))*(16.84*(1-y(1)))^0.674);
    
31. r2x=(kg2a*4.81*(1-y(2))*296452638*exp(-105858/8.314/y(3))*(4.81*(1-y(2)))^0.008)/(kg2a*4.81*(1-y(2))+296452638*exp(-105858/8.314/y(3))*(4.81*(1-y(2)))^0.008);
    
32. dy(1)=(pi*r*r*pb/Fco)*r1x;
    
33. dy(2)=(pi*r*r*pb/Fco2)*r2x;
    
34. dy(3)=(pi*r*r*pb/F)*(r1x*1000*(188.1276+0.0724*y(3)-4.2782*10^(-5)*y(3).^2+7.9724*10^(-9)*y(3).^3)+r2x*1000*(146.4633+0.06996*y(3)-2.7313*10^(-5)*y(3).^2+1.8928*10^(-9)*y(3).^3))/(1-0.14*y(1)-0.04*y(2))/(26.2932+0.01354*y(3)+6.1860*10^(-6)*y(3).^2-4.1355*10^(-9)*y(3).^3);
    

复制代码

lxx244lxx

2# qibbxxt
谢谢版主，想问下你的修改在什么地方，我看了半天没看出来。。还有就是这幅图应该怎么分析呢。。。因为刚开始学，有很多地方不懂。。