三方演化博弈复制动态方程matlab仿真——matlab2016a版本

三方演化博弈——matlab2016a版本

1.输入常微分方程组

%dx/dt=x(1-x)(0.2-yz)
%dy/dt=y(1-y)(0.2-xz)
%dz/dt=z(1-z)(0.2-xy),其中x换成x(1),y换成x(2),z换成x(3)
function dxdt=differential2(t,x)
dxdt=[x(1)*(1-x(1))*(0.2-x(2)*x(3));x(2)*(1-x(2))*(0.2-x(1)*x(3));x(3)*(1-x(3))*(0.2-x(1)*x(2))];
end
%% 尝试修改一
%function dxdt=differential2(t,x)
%dxdt=[x(1)*(1-x(1))*(-0.2-x(2)*x(3));x(2)*(1-x(2))*(-0.2-x(1)*x(3));x(3)*(1-x(3))*(0.2-x(1)*x(2))];
%end
%% 尝试修改二
%function dxdt=differential2(t,x)
%dxdt=zeros(3,1);
%dxdt(1)=x(1)*(1-x(1))*(0.2-x(2)*x(3));
%dxdt(1)=x(2)*(1-x(2))*(0.2-x(1)*x(3));
%dxdt(1)=x(3)*(1-x(3))*(0.2-x(1)*x(2));
%end
%% 尝试修改三
%function dxdt=differential2(x,t)
%dxdt=[x(1)*(1-x(1))*(0.2-x(2)*x(3)),x(2)*(1-x(2))*(0.2-x(1)*x(3)),x(3)*(1-x(3))*(0.2-x(1)*x(2))];
%end

保存为differential2.m文件

2.主函数

%% 三维图
for i=0.1:0.2:0.9    %循环步长为0.2,x的值
    for j=0.1:0.2:0.9     %循环步长为0.2,y的值
        for m=0.1:0.2:0.9    %循环步长为0.2,z的值
            [T,Y]=ode45('differential2',[0 20],[i j m]);    %常微分方程中常用ODE45（）函数，其调用格式为[T,Y]=ode45（'调用的函数名',[T自变量区间]，初始值向量，...）
            figure(1)    %图1
            grid on   %显示坐标轴网格线
            plot3(Y(:,1),Y(:,2),Y(:,3),'lineWidth',1);   %（逗号前是行，逗号前是列，冒号表示从第一行到最后一行；(Y(:,1)表示第一列的所有元素，相当于dx/dt,Y(:,2))表示第二列的所有元素相当于dy/dt，Y(:,2))表示第二列的所有元素相当于dz/dt,（'lineWidth',1）表示线宽为1
            % xlabel('x');    %x轴名称
            % ylabel('y')     %y轴名称
            % legend('ODE45 解','解析解')    曲线名称
            hold on
            view([0.5 0.5 0.5])    %调角度
        end
    end
end


%% 二维图
for i=0.1:0.2:0.9 
    for j=0.1:0.2:0.9  
        for m=0.1:0.2:0.9
            [T,Y]=ode45('differential2',[0 20],[i j m]);   
            figure(2)    %图2
            grid on   %显示坐标轴网格线
            plot(T,Y(:,1),'r-','lineWidth',1);    %x(1)和t的关系图，（'r-','lineWidth',1）表示曲线为“red”的“-”，线宽为1
            hold on
            plot(T,Y(:,2),'b--','lineWidth',1);    %x(2)和t的关系图,（'b--','lineWidth',1）表示曲线为“blue”的“--”，线宽为1
            hold on
            plot(T,Y(:,3),'g*','lineWidth',1);     %x(3)和t的关系图,（'g*','lineWidth',1）表示曲线为“green”的“*”，线宽为1
            % xlabel('x');  
            % ylabel('y')    
            % legend('ODE45 解','解析解')   
            hold on
        end
    end
end

保存为for2.m文件

3.run结果

run出来的图形跟原视频作者的图不一致，各位大神有什么发现吗，如何修改能得到原视频结果图

4.前面两个文件可修组合为一个

%% 三维图
for i=0.1:0.2:0.9    %循环步长为0.2,x的值
    for j=0.1:0.2:0.9     %循环步长为0.2,y的值
        for m=0.1:0.2:0.9    %循环步长为0.2,z的值
            f= @(t,x)[x(1)*(1-x(1))*(0.2-x(2)*x(3));x(2)*(1-x(2))*(0.2-x(1)*x(3));x(3)*(1-x(3))*(0.2-x(1)*x(2))];
            [T,Y]=ode45(f,[0 20],[i j m]);    %常微分方程中常用ODE45（）函数，其调用格式为[T,Y]=ode45（'调用的函数名',[T自变量区间]，初始值向量，...）
            figure(1)    %图1
            grid on   %显示坐标轴网格线
            plot3(Y(:,1),Y(:,2),Y(:,3),'lineWidth',1);   %（逗号前是行，逗号前是列，冒号表示从第一行到最后一行；(Y(:,1)表示第一列的所有元素，相当于dx/dt,Y(:,2))表示第二列的所有元素相当于dy/dt，Y(:,2))表示第二列的所有元素相当于dz/dt,（'lineWidth',1）表示线宽为1
            % xlabel('x');    %x轴名称
            % ylabel('y')     %y轴名称
            % legend('ODE45 解','解析解')    曲线名称
            hold on
            view([0.5 0.5 0.5])    %调角度
        end
    end
end


%% 二维图
for i=0.1:0.2:0.9 
    for j=0.1:0.2:0.9  
        for m=0.1:0.2:0.9
            f= @(t,x)[x(1)*(1-x(1))*(0.2-x(2)*x(3));x(2)*(1-x(2))*(0.2-x(1)*x(3));x(3)*(1-x(3))*(0.2-x(1)*x(2))];
            [T,Y]=ode45(f,[0 20],[i j m]);   
            figure(2)    %图2
            grid on   %显示坐标轴网格线
            plot(T,Y(:,1),'r-','lineWidth',1);    %x(1)和t的关系图，（'r-','lineWidth',1）表示曲线为“red”的“-”，线宽为1
            hold on
            plot(T,Y(:,2),'b--','lineWidth',1);    %x(2)和t的关系图,（'b--','lineWidth',1）表示曲线为“blue”的“--”，线宽为1
            hold on
            plot(T,Y(:,3),'g*','lineWidth',1);     %x(3)和t的关系图,（'g*','lineWidth',1）表示曲线为“green”的“*”，线宽为1
            % xlabel('x');  
            % ylabel('y')    
            % legend('ODE45 解','解析解')   
            hold on
        end
    end
end

保存为for2cs.m文件