matlab ode45函数演化博弈仿真，报错 索引超出组元素的数目，请帮我看下

Question

matlab ode45函数演化博弈仿真，报错 索引超出组元素的数目

以下是定义函数

function dydt=czc(y,t,c2,c3,t2,e1,e2,e3,d,j,w,f,h,r1,a,b,v,o,q1,q2)
dydt=zeros(3,1);
dydt(1)=y(1)*(1-y(1))*(y(2)*(-q1^(o)*v*f^(b)/(q1^(o)+(1-q1)^o)^(1/o)-q1^(o)*v*h^(b)/(q1^(o)+(1-q1)^o)^(1/o)+v*h^(b))+y(3)*(-q2^(o)*v*f^(b)/(q2^(o)+(1-q2)^o)^(1/o)-q2^(o)*v*h^(b)/(q2^(o)+(1-q2)^o)^(1/o)+v*h^(b))+y(2)*y(3)*(-v*f^(b)+q1^(o)*v*f^(b)/(q1^(o)+(1-q1)^o)^(1/o)+q2^(o)*v*f^(b)/(q2^(o)+(1-q2)^o)^(1/o)+q1^(o)*v*h^(b)/(q1^(o)+(1-q1)^o)^(1/o)+q2^(o)*v*h^(b)/(q2^(o)+(1-q2)^o)^(1/o)+v*h^(b))+r1-v*h^(b));
dydt(2)=y(2)*(1-y(2))*(y(1)*q1^(o)*w^(a)/(q1^(o)+(1-q1)^o)^(1/o)-y(1)*y(3)*q1^(o)*w^(a)/(q1^(o)+(1-q1)^o)^(1/o)-c2+t2-e1+e2+j);
dydt(3)=y(3)*(1-y(3))*(y(1)*q2^(o)*f^(a)/(q2^(o)+(1-q2)^o)^(1/o)+y(2)*e3+y(1)*y(2)*(f^(a)-q1^(o)*f^(a)/(q1^(o)+(1-q1)^o)^(1/o)-q2^(o)*f^(a)/(q2^(o)+(1-q2)^o)^(1/o)+q1^(o)*v*w^(b)/(q1^(o)+(1-q1)^o)^(1/o))-c3-e3+d);
end

以下是绘图过程

a=0.88,b=0.88,o=0.69,v=1.1,c2=5,c3=15,t2=10,e1=10,e2=10,e3=10,j=1,d=1,w=5,q1=0.5,q2=0.8,f=5,h=5,r1=4;
for i=0.1:0.2:1
    for j=0.1:0.2:1
        for k=0.1:0.2:1
            [t,y]=ode45('czc',[0 50],[i j k]);
            figure(1)
            grid on
            plot3(y(:,1),y(:,2),y(:,3),'linewidth',1);
            set(gca,'XTick',[0:0.2:1],'YTick',[0:0.2:1],'ZTick',[0:0.2:1])
            hold on
            axis([0 1 0 1 0 1])
            view([45 10])
        end
    end
end
xlabel('x','Rotation',0);
ylabel('y','Rotation',0);
zlabel('z','Rotation',360,'position',[0 0 1.05]);

以下是报错内容

索引超出数组元素的数目(1)。

出错 czc (第 3 行)
dydt(1)=y(1)*(1-y(1))*(y(2)*(-q1^(o)*v*f^(b)/(q1^(o)+(1-q1)^o)^(1/o)-q1^(o)*v*h^(b)/(q1^(o)+(1-q1)^o)^(1/o)+v*h^(b))+y(3)*(-q2^(o)*v*f^(b)/(q2^(o)+(1-q2)^o)^(1/o)-q2^(o)*v*h^(b)/(q2^(o)+(1-q2)^o)^(1/o)+v*h^(b))+y(2)*y(3)*(-v*f^(b)+q1^(o)*v*f^(b)/(q1^(o)+(1-q1)^o)^(1/o)+q2^(o)*v*f^(b)/(q2^(o)+(1-q2)^o)^(1/o)+q1^(o)*v*h^(b)/(q1^(o)+(1-q1)^o)^(1/o)+q2^(o)*v*h^(b)/(q2^(o)+(1-q2)^o)^(1/o)+v*h^(b))+r1-v*h^(b));

出错 odearguments (第 90 行)
f0 = feval(ode,t0,y0,args{:});   % ODE15I sets args{1} to yp0.

出错 ode45 (第 115 行)
  odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

出错 czctu (第 8 行)
            [t,y]=ode45('czc',[0 50],[i j k]);

对于定义函数中的符号，实际上是想完成以下运算

辛苦了，请帮我解答一下

Answer 1

你好，帮你更正如下，主要是改了ode45调用czc函数的方式，我修改之后的格式是最经典，最能反应matlab求ode本质的了

a=0.88; b=0.88; o=0.69; v=1.1; c2=5; c3=15;
t2=10; e1=10; e2=10; e3=10; j=1; d=1; w=5; 
q1=0.5; q2=0.8; f=5; h=5; r1=4;
for i=0.1:0.2:1
    for j=0.1:0.2:1
        for k=0.1:0.2:1
            [t,y]=ode45(@(t,y)czc(y,t,c2,c3,t2,e1,e2,e3,d,j,w,f,h,r1,a,b,v,o,q1,q2),[0 50],[i j k]); % 这里修改了一下，一般是先t后y，然后其它参数放后面
            figure(1)
            grid on
            plot3(y(:,1),y(:,2),y(:,3),'linewidth',1);
            set(gca,'XTick',[0:0.2:1],'YTick',[0:0.2:1],'ZTick',[0:0.2:1])
            hold on
            axis([0 1 0 1 0 1])
            view([45 10])
        end
    end
end
xlabel('x','Rotation',0);
ylabel('y','Rotation',0);
zlabel('z','Rotation',360,'position',[0 0 1.05]);

function dydt=czc(y,t,c2,c3,t2,e1,e2,e3,d,j,w,f,h,r1,a,b,v,o,q1,q2)
dydt=zeros(3,1);
dydt(1)=y(1)*(1-y(1))*(y(2)*(-q1^(o)*v*f^(b)/(q1^(o)+(1-q1)^o)^(1/o)-q1^(o)*v*h^(b)/(q1^(o)+(1-q1)^o)^(1/o)+v*h^(b))+y(3)*(-q2^(o)*v*f^(b)/(q2^(o)+(1-q2)^o)^(1/o)-q2^(o)*v*h^(b)/(q2^(o)+(1-q2)^o)^(1/o)+v*h^(b))+y(2)*y(3)*(-v*f^(b)+q1^(o)*v*f^(b)/(q1^(o)+(1-q1)^o)^(1/o)+q2^(o)*v*f^(b)/(q2^(o)+(1-q2)^o)^(1/o)+q1^(o)*v*h^(b)/(q1^(o)+(1-q1)^o)^(1/o)+q2^(o)*v*h^(b)/(q2^(o)+(1-q2)^o)^(1/o)+v*h^(b))+r1-v*h^(b));
dydt(2)=y(2)*(1-y(2))*(y(1)*q1^(o)*w^(a)/(q1^(o)+(1-q1)^o)^(1/o)-y(1)*y(3)*q1^(o)*w^(a)/(q1^(o)+(1-q1)^o)^(1/o)-c2+t2-e1+e2+j);
dydt(3)=y(3)*(1-y(3))*(y(1)*q2^(o)*f^(a)/(q2^(o)+(1-q2)^o)^(1/o)+y(2)*e3+y(1)*y(2)*(f^(a)-q1^(o)*f^(a)/(q1^(o)+(1-q1)^o)^(1/o)-q2^(o)*f^(a)/(q2^(o)+(1-q2)^o)^(1/o)+q1^(o)*v*w^(b)/(q1^(o)+(1-q1)^o)^(1/o))-c3-e3+d);
end

效果：