三方演化博弈复制动态方程matlab仿真（输入参数较多时）——matlab2016a版本

%作者：应该就我叫测差
%链接：https://zhuanlan.zhihu.com/p/363605531
%来源：知乎
%著作权归作者所有。商业转载请联系作者获得授权，非商业转载请注明出处。

1、function函数声明

function dydt=yaopin(t,y,rp,cph,cpl,cp,bt,fp,mp,ct,ft,mt,cg,tg)
dydt=zeros(3,1);
dydt(1)=y(1)*(y(1)-1)*(cph-cpl-cp-bt-y(2)*(rp-bt)-y(3)*(fp+mp));
dydt(2)=y(2)*(y(2)-1)*((1-y(1))*(bt-mt)-(ft+mt)*y(3)-ct);
dydt(3)=y(3)*(y(3)-1)*(cg-fp-ft-tg+(mp+fp+tg)*y(1)+(mt+ft+tg)*y(2)-tg*y(1)*y(2));
end

保存为yaopi.m文件

2.新建脚本文件

%作者：应该就我叫测差
%链接：https://zhuanlan.zhihu.com/p/363605531
%来源：知乎
%著作权归作者所有。商业转载请联系作者获得授权，非商业转载请注明出处。

%% 画图1中大图案“红+”的线
rp=100;
cph=85,cpl=0,cp=10,bt=40,fp=40,mp=20,ct=10,ft=20,mt=15,cg=15,tg=40;  %参数设置（输入）：function dydt=yaopin(t,y,rp,cph,cpl,cp,bt,fp,mp,ct,ft,mt,cg,tg)
for i=0.2
    figure(1)
    %set(0,'defaultfigurecolor','w') %将图1的背景色设置为白色,此处试验过，背景色没有改变。因此更改为whitebg('blac')
    whitebg('black') %自己不需要可以修改成white
    [t,y]=ode45(@(t,y) yaopin(t,y,rp,cph,cpl,cp,bt,fp,mp,ct,ft,mt,cg,tg),[0 50],[i 0.2 0.2]);
    %ode45（）：求解微分方程的函数；@(t,y) yaopin(t,y,rp,cph,cpl,cp,bt,fp,mp,ct,ft,mt,cg,tg)：函数句柄调用；[0 50]：表示t的范围；[i 0.2 0.2]：i指的是x，第一个0.2指的是y的值，第二个0.2指的是z的值
    plot3(y(:,1),y(:,2),y(:,3),'r+','linewidth',1);
    hold on
end
%% 画图1中“蓝-”的线
rp=150;
cph=85,cpl=0,cp=10,bt=40,fp=40,mp=20,ct=10,ft=20,mt=15,cg=15,tg=40;  
for j=0.2
    [t,y]=ode45(@(t,y) yaopin(t,y,rp,cph,cpl,cp,bt,fp,mp,ct,ft,mt,cg,tg),[0 50],[0.2 j 0.2]);
    plot3(y(:,1),y(:,2),y(:,3),'b-','linewidth',1);
    hold on     
end  
%% 画图1中“绿--”的线
rp=200;
cph=85,cpl=0,cp=10,bt=40,fp=40,mp=20,ct=10,ft=20,mt=15,cg=15,tg=40;  
for k=0.2
    [t,y]=ode45(@(t,y) yaopin(t,y,rp,cph,cpl,cp,bt,fp,mp,ct,ft,mt,cg,tg),[0 50],[0.2 0.2 k]);
    plot3(y(:,1),y(:,2),y(:,3),'g--','linewidth',1);
    axis([0 1 0 1 0 1])  %坐标轴x,y,z的控制范围
    set(gca,'XTick',[0:0.2:1],'YTick',[0:0.2:1],'ZTick',[0:0.2:1])
    xlabel('x');ylabel('y');zlabel('z','Rotation',360,'position',[-0.1 1 1.1]);
    grid on
    hold on  
    legend({'Rp=100','Rp=150','Rp=200'},'location','northeast');
    title('图  5 药品销售收入的影响','FontWeight','bold','position',[0 0 -0.2]);
end


%% 画图1中小图案的线
axes('position',[0.13 0.32 0.2 0.2]);
%axis off
grid on
hold on
rp=100,cph=85,cpl=0,cp=10,bt=40,fp=40,mp=20,ct=10,ft=20,mt=15,cg=15,tg=40;  
for i=0.2
    figure(1)
    set(0,'defaultfigurecolor','w')
    [t,y]=ode45(@(t,y) yaopin(t,y,rp,cph,cpl,cp,bt,fp,mp,ct,ft,mt,cg,tg),[0 50],[i 0.2 0.2]);
    plot3(y(:,1),y(:,2),y(:,3),'r+','linewidth',1);
    hold on 
end
rp=150,cph=85,cpl=0,cp=10,bt=40,fp=40,mp=20,ct=10,ft=20,mt=15,cg=15,tg=40;  
for j=0.2
    [t,y]=ode45(@(t,y) yaopin(t,y,rp,cph,cpl,cp,bt,fp,mp,ct,ft,mt,cg,tg),[0 50],[0.2 j 0.2]);
    plot3(y(:,1),y(:,2),y(:,3),'b-','linewidth',1);
    hold on     
end  
rp=200,cph=85,cpl=0,cp=10,bt=40,fp=40,mp=20,ct=10,ft=20,mt=15,cg=15,tg=40;  
for k=0.2
    [t,y]=ode45(@(t,y) yaopin(t,y,rp,cph,cpl,cp,bt,fp,mp,ct,ft,mt,cg,tg),[0 50],[0.2 0.2 k]);
    plot3(y(:,1),y(:,2),y(:,3),'g--','linewidth',1);
    axis([0 1 0 1 0 1])
    set(gca,'XTick',[0:0.2:1],'YTick',[0:0.2:1],'ZTick',[0:0.2:1])
    xlabel('x');ylabel('y');zlabel('z','Rotation',360);
    grid on
    hold on  
end
view([0 0]);
xlim([0 1]);ylim([0 1]);
set(gca,'XTickLabel','','YTickLabel','','ZTickLabel','')
xlabel('x','position',[0.8 0.4 0.3])
%ylabel('y','position',[0.8 0.4 0.3],'Rotation',360)
zlabel('z','position',[0.1 0.8 0.8],'Rotation',360)
hold on 

%% 3D图
clear; 
rp=150,cph=85,cpl=0,cp=10,bt=40,fp=40,mp=20,ct=10,ft=20,mt=15,cg=15,tg=40;  
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) yaopin(t,y,rp,cph,cpl,cp,bt,fp,mp,ct,ft,mt,cg,tg),[0 50],[i j k]);
        figure(2)
        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]);
title('图  11 数组1演化50次结果','FontWeight','bold','position',[1 0 -0.13]);

保存为yaopinsanshitu.m文件

3、run出图形