3.三方演化博弈代码

论文标题：政府奖惩机制下药品安全质量监管三方演化博弈及仿真分析

视频可关注bilibili：谭小飞同学

yaopin.m

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

yaopinsanshitu.m

% the 1st line
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
% the 2nd line
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  
% the 3rd line
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,'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
% the small figure
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 figure
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]);