美文网首页
4. 四方演化博弈代码

4. 四方演化博弈代码

作者: 谭小飞同学 | 来源:发表于2021-05-07 08:27 被阅读0次

    论文标题：重大疫情期医药研究报道质量监管四方演化博弈分析

    关注b站：谭小飞同学

    %sifang.m

    % 特殊的希腊字符换成其读音字母表示，文中的r,m,g,p依次表示为y(1),y(2),y(3),y(4)
% 论文标题为《重大疫情期医药研究报道质量监管四方演化博弈分析》
% 由于公式较长，本人虽然几番认真仔细核对，难免也可能出错，建议读者自行录入公式，看看是否有出入
function dydt=sifang(t,y,a,phi,b,Rr,Fr,Crh,Crl,Rm,Fm,psi,Cml,Cmh,Cgl,Cgh,Ng,Tg,Np,Tp,Cp)
dydt=zeros(4,1);
dydt(1)=a*y(1)*(1-y(1))*((1-(1-y(2))*(1-y(3))*(1-y(4)))*phi+(1-y(2))*(1-y(3))*(y(4)+(1-y(4))*b)*Rr+(1-y(2))*y(3)*Fr-Crh+Crl);
dydt(2)=a*y(2)*(1-y(2))*((1-y(1))*(1-y(3))*(y(4)+(1-y(4))*b)*Rm+(1-y(1))*y(3)*Fm-(1-y(1))*(1-y(3))*(1-y(4))*psi+Cml-Cmh);
dydt(3)=a*y(3)*(1-y(3))*((y(1)+(1-y(1))*(1-y(2)))*(Cgl-Cgh)+(1-y(1))*(1-y(2))*(Fr+Fm+(1-y(4))*(Ng-b*Tg)));
dydt(4)=a*y(4)*(1-y(4))*((1-y(1))*(1-y(2))*(1-y(3))*(Np-b*Tp-Cp)-y(1)*Cp);
end

    sifang2.m

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%政府大图
figure(1)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%子图1
clc;clear;
a=0.7,phi=22,b=0.3,Rr=25,Fr=10,Crh=10,Crl=2,Rm=15,Fm=6,psi=12,Cml=1,Cmh=5,Cgl=1,Cgh=12,Ng=20,Tg=10,Np=12,Tp=5,Cp=10;
%subplot(3,1,1)
set(0,'defaultfigurecolor','w')
[t,y]=ode45(@(t,y) sifang(t,y,a,phi,b,Rr,Fr,Crh,Crl,Rm,Fm,psi,Cml,Cmh,Cgl,Cgh,Ng,Tg,Np,Tp,Cp),[0,200],[0.4,0.3,0.2,0.3]);
points=1:1:length(t);
plot(t,y(:,1),'r^-','linewidth',1,'markersize',3,'markerfacecolor','r','markerindices',points);
% hold on放在函数ode45前面的话，生成的图像不是封闭的方框，而是坐标系的第一象限
hold on
plot(t,y(:,2),'b-','linewidth',1);
hold on
plot(t,y(:,3),'y-.','linewidth',1);
hold on
plot(t,y(:,4),'g--','linewidth',1);
hold on
set(gca,'XTick',[0:50:200],'YTick',[0.0:0.2:1.0])
set(gca,'YTickLabel',num2str(get(gca,'YTick')','%.1f'));
axis([0 200 -0.05 1.05])
xlabel('$t$','interpreter','latex');
ylabel('概率');
zhuti=title('$C_{gh}=12$');
set(zhuti,'interpreter','latex')
legend('医药研究机构({\it\fontname{Bodoni MT}r})','媒体平台({\it\fontname{Bodoni MT}m})','政府监管部门({\it\fontname{Bodoni MT}g})','群众({\it\fontname{Bodoni MT}p})');
%%%%%%%%%%%%%%%子图2
clc;clear;
a=0.7,phi=22,b=0.3,Rr=25,Fr=10,Crh=10,Crl=2,Rm=15,Fm=6,psi=12,Cml=1,Cmh=5,Cgl=1,Cgh=6,Ng=20,Tg=10,Np=12,Tp=5,Cp=10;
figure(2)
%subplot(3,1,2)
set(0,'defaultfigurecolor','w')
[t,y]=ode45(@(t,y) sifang(t,y,a,phi,b,Rr,Fr,Crh,Crl,Rm,Fm,psi,Cml,Cmh,Cgl,Cgh,Ng,Tg,Np,Tp,Cp),[0,200],[0.5,0.3,0.2,0.3]);
points=1:1:length(t);
plot(t,y(:,1),'r^-','linewidth',1,'markersize',3,'markerfacecolor','r','markerindices',points);
hold on
plot(t,y(:,2),'b-','linewidth',1);
hold on
plot(t,y(:,3),'y-.','linewidth',1);
hold on
plot(t,y(:,4),'g--','linewidth',1);
hold on
set(gca,'XTick',[0:50:200],'YTick',[0.0:0.2:1.0])
set(gca,'YTickLabel',num2str(get(gca,'YTick')','%.1f'));
axis([0 200 -0.05 1.05])
xlabel('$t$','interpreter','latex');
ylabel('概率');
zhuti=title('$C_{gh}=6$');
set(zhuti,'interpreter','latex')
legend('医药研究机构({\it\fontname{Bodoni MT}r})','媒体平台({\it\fontname{Bodoni MT}m})','政府监管部门({\it\fontname{Bodoni MT}g})','群众({\it\fontname{Bodoni MT}p})');
%%%%%%%%%%%%%%%子图3
clc;clear;
a=0.7,phi=22,b=0.3,Rr=25,Fr=10,Crh=10,Crl=2,Rm=15,Fm=6,psi=12,Cml=1,Cmh=5,Cgl=1,Cgh=2,Ng=20,Tg=10,Np=12,Tp=5,Cp=10;
figure(3)
%subplot(3,1,3)
set(0,'defaultfigurecolor','w')
[t,y]=ode45(@(t,y) sifang(t,y,a,phi,b,Rr,Fr,Crh,Crl,Rm,Fm,psi,Cml,Cmh,Cgl,Cgh,Ng,Tg,Np,Tp,Cp),[0,200],[0.5,0.3,0.2,0.3]);
points=1:1:length(t);
plot(t,y(:,1),'r^-','linewidth',1,'markersize',3,'markerfacecolor','r','markerindices',points);
hold on
plot(t,y(:,2),'b-','linewidth',1);
hold on
plot(t,y(:,3),'y-.','linewidth',1);
hold on
plot(t,y(:,4),'g--','linewidth',1);
hold on
set(gca,'XTick',[0:50:200],'YTick',[0.0:0.2:1.0])
set(gca,'YTickLabel',num2str(get(gca,'YTick')','%.1f'));
axis([0 200 -0.05 1.05])
xlabel('$t$','interpreter','latex');
ylabel('概率');
zhuti=title('$C_{gh}=2$');
set(zhuti,'interpreter','latex')
legend('医药研究机构({\it\fontname{Bodoni MT}r})','媒体平台({\it\fontname{Bodoni MT}m})','政府监管部门({\it\fontname{Bodoni MT}g})','群众({\it\fontname{Bodoni MT}p})');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%群众大图
figure(4)
%%%%%%%%%%%%%%%子图1
clc;clear;
a=0.7,phi=22,b=0.3,Rr=12,Fr=10,Crh=10,Crl=2,Rm=6,Fm=6,psi=12,Cml=1,Cmh=5,Cgl=1,Cgh=6,Ng=20,Tg=10,Np=12,Tp=5,Cp=10;
%subplot(3,1,1)
set(0,'defaultfigurecolor','w')
[t,y]=ode45(@(t,y) sifang(t,y,a,phi,b,Rr,Fr,Crh,Crl,Rm,Fm,psi,Cml,Cmh,Cgl,Cgh,Ng,Tg,Np,Tp,Cp),[0,100],[0.4,0.3,0.2,0.3]);
points=1:1:length(t);
plot(t,y(:,1),'r^-','linewidth',1,'markersize',3,'markerfacecolor','r','markerindices',points);
hold on
plot(t,y(:,2),'b-','linewidth',1);
hold on
plot(t,y(:,3),'y-.','linewidth',1);
hold on
plot(t,y(:,4),'g--','linewidth',1);
hold on
set(gca,'XTick',[0:20:100],'YTick',[0.0:0.2:1.0])
set(gca,'YTickLabel',num2str(get(gca,'YTick')','%.1f'));
axis([0 100 -0.05 1.05])
xlabel('$t$','interpreter','latex');
ylabel('概率');
zhuti=title('$R_{r}=12,R_{m}=6$');
set(zhuti,'interpreter','latex')
legend('医药研究机构({\it\fontname{Bodoni MT}r})','媒体平台({\it\fontname{Bodoni MT}m})','政府监管部门({\it\fontname{Bodoni MT}g})','群众({\it\fontname{Bodoni MT}p})');
%%%%%%%%%%%%%%%子图2
clc;clear;
a=0.7,phi=22,b=0.3,Rr=25,Fr=10,Crh=10,Crl=2,Rm=15,Fm=6,psi=12,Cml=1,Cmh=5,Cgl=1,Cgh=6,Ng=20,Tg=10,Np=12,Tp=5,Cp=10;
figure(5)
%subplot(3,1,2)
set(0,'defaultfigurecolor','w')
[t,y]=ode45(@(t,y) sifang(t,y,a,phi,b,Rr,Fr,Crh,Crl,Rm,Fm,psi,Cml,Cmh,Cgl,Cgh,Ng,Tg,Np,Tp,Cp),[0,100],[0.5,0.3,0.2,0.3]);
points=1:1:length(t);
plot(t,y(:,1),'r^-','linewidth',1,'markersize',3,'markerfacecolor','r','markerindices',points);
hold on
plot(t,y(:,2),'b-','linewidth',1);
hold on
plot(t,y(:,3),'y-.','linewidth',1);
hold on
plot(t,y(:,4),'g--','linewidth',1);
hold on
set(gca,'XTick',[0:20:100],'YTick',[0.0:0.2:1.0])
set(gca,'YTickLabel',num2str(get(gca,'YTick')','%.1f'));
axis([0 100 -0.05 1.05])
xlabel('$t$','interpreter','latex');
ylabel('概率');
zhuti=title('$R_{r}=25,R_{m}=15$');
set(zhuti,'interpreter','latex')
legend('医药研究机构({\it\fontname{Bodoni MT}r})','媒体平台({\it\fontname{Bodoni MT}m})','政府监管部门({\it\fontname{Bodoni MT}g})','群众({\it\fontname{Bodoni MT}p})');
%%%%%%%%%%%%%%%子图3
clc;clear;
a=0.7,phi=22,b=0.3,Rr=50,Fr=10,Crh=10,Crl=2,Rm=35,Fm=6,psi=12,Cml=1,Cmh=5,Cgl=1,Cgh=6,Ng=20,Tg=10,Np=12,Tp=5,Cp=10;
figure(6)
%subplot(3,1,3)
set(0,'defaultfigurecolor','w')
[t,y]=ode45(@(t,y) sifang(t,y,a,phi,b,Rr,Fr,Crh,Crl,Rm,Fm,psi,Cml,Cmh,Cgl,Cgh,Ng,Tg,Np,Tp,Cp),[0,100],[0.5,0.3,0.2,0.3]);
points=1:1:length(t);
plot(t,y(:,1),'r^-','linewidth',1,'markersize',3,'markerfacecolor','r','markerindices',points);
hold on
plot(t,y(:,2),'b-','linewidth',1);
hold on
plot(t,y(:,3),'y-.','linewidth',1);
hold on
plot(t,y(:,4),'g--','linewidth',1);
hold on
%下面一行是加方框的，不过也可以在成图界面去手动制图
%rectangle('position',[0 0 2 1],'edgecolor','g','linewidth',1)
hold on
set(gca,'XTick',[0:20:100],'YTick',[0.0:0.2:1.0])
set(gca,'YTickLabel',num2str(get(gca,'YTick')','%.1f'));
axis([0 100 -0.05 1.05])
xlabel('$t$','interpreter','latex');
ylabel('概率');
zhuti=title('$R_{r}=50,R_{m}=35$');
set(zhuti,'interpreter','latex')
legend('医药研究机构({\it\fontname{Bodoni MT}r})','媒体平台({\it\fontname{Bodoni MT}m})','政府监管部门({\it\fontname{Bodoni MT}g})','群众({\it\fontname{Bodoni MT}p})');
%%%%%%%%%%%%%%%%%% 群众第三张图的小图
%下面一行决定小图的位置及其大小
axes('position',[0.35 0.2 0.2 0.2]);
clc;clear;
a=0.7,phi=22,b=0.3,Rr=50,Fr=10,Crh=10,Crl=2,Rm=35,Fm=6,psi=12,Cml=1,Cmh=5,Cgl=1,Cgh=6,Ng=20,Tg=10,Np=12,Tp=5,Cp=10;
set(0,'defaultfigurecolor','w')
[t,y]=ode45(@(t,y) sifang(t,y,a,phi,b,Rr,Fr,Crh,Crl,Rm,Fm,psi,Cml,Cmh,Cgl,Cgh,Ng,Tg,Np,Tp,Cp),[0,100],[0.5,0.3,0.2,0.3]);
points=1:1:length(t);
plot(t,y(:,1),'r^-','linewidth',1,'markersize',3,'markerfacecolor','r','markerindices',points);
hold on
plot(t,y(:,2),'b-','linewidth',1);
hold on
plot(t,y(:,3),'y-.','linewidth',1);
hold on
plot(t,y(:,4),'g--','linewidth',1);
hold on
set(gca,'XTick',[0:1:3],'YTick',[0:0.2:1])
hold on
set(gca,'YTickLabel',num2str(get(gca,'YTick')','%.1f'));
axis([0 3 -0.05 1.05])
hold on 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%民间组织大图
figure(7)
%%%%%%%%%%%%%%%子图1
clc;clear;
a=0.7,phi=22,b=0.3,Rr=25,Fr=10,Crh=10,Crl=2,Rm=15,Fm=6,psi=12,Cml=1,Cmh=5,Cgl=1,Cgh=6,Ng=20,Tg=10,Np=12,Tp=5,Cp=10;
%subplot(3,1,1)
set(0,'defaultfigurecolor','w')
[t,y]=ode45(@(t,y) sifang(t,y,a,phi,b,Rr,Fr,Crh,Crl,Rm,Fm,psi,Cml,Cmh,Cgl,Cgh,Ng,Tg,Np,Tp,Cp),[0,100],[0.5,0.3,0.2,0.3]);
points=1:1:length(t);
plot(t,y(:,1),'r^-','linewidth',1,'markersize',3,'markerfacecolor','r','markerindices',points);
hold on
plot(t,y(:,2),'b-','linewidth',1);
hold on
plot(t,y(:,3),'y-.','linewidth',1);
hold on
plot(t,y(:,4),'g--','linewidth',1);
hold on
set(gca,'XTick',[0:20:100],'YTick',[0.0:0.2:1.0])
set(gca,'YTickLabel',num2str(get(gca,'YTick')','%.1f'));
axis([0 100 -0.05 1.05])
xlabel('$t$','interpreter','latex');
ylabel('概率');
zhuti=title('$\beta=0.3$');
set(zhuti,'interpreter','latex')
legend('医药研究机构({\it\fontname{Bodoni MT}r})','媒体平台({\it\fontname{Bodoni MT}m})','政府监管部门({\it\fontname{Bodoni MT}g})','群众({\it\fontname{Bodoni MT}p})');
%%%%%%%%%%%%%%%子图2
clc;clear;
a=0.7,phi=22,b=0.31,Rr=25,Fr=10,Crh=10,Crl=2,Rm=15,Fm=6,psi=12,Cml=1,Cmh=5,Cgl=1,Cgh=6,Ng=20,Tg=10,Np=12,Tp=5,Cp=10;
figure(8)
%subplot(3,1,2)
set(0,'defaultfigurecolor','w')
[t,y]=ode45(@(t,y) sifang(t,y,a,phi,b,Rr,Fr,Crh,Crl,Rm,Fm,psi,Cml,Cmh,Cgl,Cgh,Ng,Tg,Np,Tp,Cp),[0,100],[0.5,0.3,0.2,0.3]);
points=1:1:length(t);
plot(t,y(:,1),'r^-','linewidth',1,'markersize',3,'markerfacecolor','r','markerindices',points);
hold on
plot(t,y(:,2),'b-','linewidth',1);
hold on
plot(t,y(:,3),'y-.','linewidth',1);
hold on
plot(t,y(:,4),'g--','linewidth',1);
hold on
set(gca,'XTick',[0:20:100],'YTick',[0.0:0.2:1.0])
set(gca,'YTickLabel',num2str(get(gca,'YTick')','%.1f'));
axis([0 100 -0.05 1.05])
xlabel('$t$','interpreter','latex');
ylabel('概率');
zhuti=title('$\beta=0.31$');
set(zhuti,'interpreter','latex')
legend('医药研究机构({\it\fontname{Bodoni MT}r})','媒体平台({\it\fontname{Bodoni MT}m})','政府监管部门({\it\fontname{Bodoni MT}g})','群众({\it\fontname{Bodoni MT}p})');
%%%%%%%%%%%%%%%子图3
clc;clear;
a=0.7,phi=22,b=0.5,Rr=50,Fr=10,Crh=10,Crl=2,Rm=35,Fm=6,psi=12,Cml=1,Cmh=5,Cgl=1,Cgh=6,Ng=20,Tg=10,Np=12,Tp=5,Cp=10;
figure(9)
%subplot(3,1,3)
set(0,'defaultfigurecolor','w')
[t,y]=ode45(@(t,y) sifang(t,y,a,phi,b,Rr,Fr,Crh,Crl,Rm,Fm,psi,Cml,Cmh,Cgl,Cgh,Ng,Tg,Np,Tp,Cp),[0,100],[0.5,0.3,0.2,0.3]);
points=1:1:length(t);
plot(t,y(:,1),'r^-','linewidth',1,'markersize',3,'markerfacecolor','r','markerindices',points);
hold on
plot(t,y(:,2),'b-','linewidth',1);
hold on
plot(t,y(:,3),'y-.','linewidth',1);
hold on
plot(t,y(:,4),'g--','linewidth',1);
hold on
%又到了加方框的时候了
%rectangle('position',[0 0 2 1],'edgecolor','g','linewidth',1)
hold on
set(gca,'XTick',[0:20:100],'YTick',[0.0:0.2:1.0])
set(gca,'YTickLabel',num2str(get(gca,'YTick')','%.1f'));
axis([0 100 -0.05 1.05])
xlabel('$t$','interpreter','latex');
ylabel('概率');
zhuti=title('$\beta=0.5$');
set(zhuti,'interpreter','latex')
legend('医药研究机构({\it\fontname{Bodoni MT}r})','媒体平台({\it\fontname{Bodoni MT}m})','政府监管部门({\it\fontname{Bodoni MT}g})','群众({\it\fontname{Bodoni MT}p})');
%%%%%%%%%%%%%%%%%% 民间组织第三张图的小图
axes('position',[0.35 0.2 0.2 0.2]);
clc;clear;
a=0.7,phi=22,b=0.5,Rr=50,Fr=10,Crh=10,Crl=2,Rm=35,Fm=6,psi=12,Cml=1,Cmh=5,Cgl=1,Cgh=6,Ng=20,Tg=10,Np=12,Tp=5,Cp=10;
set(0,'defaultfigurecolor','w')
[t,y]=ode45(@(t,y) sifang(t,y,a,phi,b,Rr,Fr,Crh,Crl,Rm,Fm,psi,Cml,Cmh,Cgl,Cgh,Ng,Tg,Np,Tp,Cp),[0,100],[0.5,0.3,0.2,0.3]);
points=1:1:length(t);
plot(t,y(:,1),'r^-','linewidth',1,'markersize',3,'markerfacecolor','r','markerindices',points);
hold on
plot(t,y(:,2),'b-','linewidth',1);
hold on
plot(t,y(:,3),'y-.','linewidth',1);
hold on
plot(t,y(:,4),'g--','linewidth',1);
hold on
set(gca,'XTick',[0:1:3],'YTick',[0.0:0.2:1.0])
set(gca,'YTickLabel',num2str(get(gca,'YTick')','%.1f'));
axis([0 3 -0.05 1.05])
hold on 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%3d figure
%%%%%%%%%%%%%%%%%%%%%%%%%g=0的子图
clc;clear; 
a=0.7,phi=22,b=0.3,Rr=25,Fr=10,Crh=10,Crl=2,Rm=15,Fm=6,psi=12,Cml=1,Cmh=5,Cgl=1,Cgh=6,Ng=20,Tg=10,Np=12,Tp=5,Cp=10;  
figure(10)
%subplot(2,1,1)
for i=0.1:0.2:1
    for j=0.1:0.2:1
        for k=0
            for l=0.1:0.2:1
            [t,y]=ode45(@(t,y) sifang(t,y,a,phi,b,Rr,Fr,Crh,Crl,Rm,Fm,psi,Cml,Cmh,Cgl,Cgh,Ng,Tg,Np,Tp,Cp),[0 50],[i j k l]);
            grid on
            %搞清楚你画的是哪三维的关系
            plot3(y(:,1),y(:,2),y(:,4),'linewidth',1);
            set(gca,'XTick',[0:0.2:1],'YTick',[0:0.2:1],'ZTick',[0:0.2:1])
            set(gca,'XTickLabel',num2str(get(gca,'XTick')','%.1f'));
            set(gca,'YTickLabel',num2str(get(gca,'YTick')','%.1f'));
            set(gca,'ZTickLabel',num2str(get(gca,'ZTick')','%.1f'));
            hold on
            axis([0 1 0 1 0 1])
            end
        end
    end
end
xlabel('$r$','interpreter','latex');
ylabel('$m$','interpreter','latex');
zlabel('$p$','interpreter','latex','Rotation',360);
title('g=0','interpreter','latex');
%%%%%%%%%%%%%%%%%%%%%%%%%g=0.8的子图
clc;clear; 
a=0.7,phi=22,b=0.3,Rr=25,Fr=10,Crh=10,Crl=2,Rm=15,Fm=6,psi=12,Cml=1,Cmh=5,Cgl=1,Cgh=6,Ng=20,Tg=10,Np=12,Tp=5,Cp=10;  
figure(11)
%subplot(2,1,2)
for i=0.1:0.2:1
    for j=0.1:0.2:1
        for k=0.8
            for l=0.1:0.2:1
            [t,y]=ode45(@(t,y) sifang(t,y,a,phi,b,Rr,Fr,Crh,Crl,Rm,Fm,psi,Cml,Cmh,Cgl,Cgh,Ng,Tg,Np,Tp,Cp),[0 50],[i j k l]);
            grid on
            %搞清楚你画的是哪三维的关系
            plot3(y(:,1),y(:,2),y(:,4),'linewidth',1);
            set(gca,'XTick',[0:0.2:1],'YTick',[0:0.2:1],'ZTick',[0:0.2:1])
            set(gca,'XTickLabel',num2str(get(gca,'XTick')','%.1f'));
            set(gca,'YTickLabel',num2str(get(gca,'YTick')','%.1f'));
            set(gca,'ZTickLabel',num2str(get(gca,'ZTick')','%.1f'));
            hold on
            axis([0 1 0 1 0 1])
            end
        end
    end
end
xlabel('$r$','interpreter','latex');
ylabel('$m$','interpreter','latex');
zlabel('$p$','interpreter','latex','Rotation',360);
title('{\it\fontname{Bodoni MT}g}=0.8');
    image.png
    image.png
    image.png
    image.png
    image.png
    image.png
    image.png
    image.png image.png
    image.png
    image.png

    请注意，上面图片与原论文中的样图还是有很大的差别，下面我们使用subplot()函数及MATLAB的图像编辑器来逼近原样图（具体操作看上方链接视频）。

    稍微改变下代码如下：

    % 特殊的希腊字符换成其读音字母表示，文中的r,m,g,p依次表示为y(1),y(2),y(3),y(4)
% 论文标题为《重大疫情期医药研究报道质量监管四方演化博弈分析》
% 由于公式较长，本人虽然几番认真仔细核对，难免也可能出错，建议读者自行录入公式，看看是否有出入
function dydt=sifang(t,y,a,phi,b,Rr,Fr,Crh,Crl,Rm,Fm,psi,Cml,Cmh,Cgl,Cgh,Ng,Tg,Np,Tp,Cp)
dydt=zeros(4,1);
dydt(1)=a*y(1)*(1-y(1))*((1-(1-y(2))*(1-y(3))*(1-y(4)))*phi+(1-y(2))*(1-y(3))*(y(4)+(1-y(4))*b)*Rr+(1-y(2))*y(3)*Fr-Crh+Crl);
dydt(2)=a*y(2)*(1-y(2))*((1-y(1))*(1-y(3))*(y(4)+(1-y(4))*b)*Rm+(1-y(1))*y(3)*Fm-(1-y(1))*(1-y(3))*(1-y(4))*psi+Cml-Cmh);
dydt(3)=a*y(3)*(1-y(3))*((y(1)+(1-y(1))*(1-y(2)))*(Cgl-Cgh)+(1-y(1))*(1-y(2))*(Fr+Fm+(1-y(4))*(Ng-b*Tg)));
dydt(4)=a*y(4)*(1-y(4))*((1-y(1))*(1-y(2))*(1-y(3))*(Np-b*Tp-Cp)-y(1)*Cp);
end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%政府大图
figure(1)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%子图1
clc;clear;
a=0.7,phi=22,b=0.3,Rr=25,Fr=10,Crh=10,Crl=2,Rm=15,Fm=6,psi=12,Cml=1,Cmh=5,Cgl=1,Cgh=12,Ng=20,Tg=10,Np=12,Tp=5,Cp=10;
subplot(3,1,1)
set(0,'defaultfigurecolor','w')
[t,y]=ode45(@(t,y) sifang(t,y,a,phi,b,Rr,Fr,Crh,Crl,Rm,Fm,psi,Cml,Cmh,Cgl,Cgh,Ng,Tg,Np,Tp,Cp),[0,200],[0.4,0.3,0.2,0.3]);
points=1:1:length(t);
plot(t,y(:,1),'r^-','linewidth',1,'markersize',3,'markerfacecolor','r','markerindices',points);
% hold on放在函数ode45前面的话，生成的图像不是封闭的方框，而是坐标系的第一象限
hold on
plot(t,y(:,2),'b-','linewidth',1);
hold on
plot(t,y(:,3),'y-.','linewidth',1);
hold on
plot(t,y(:,4),'g--','linewidth',1);
hold on
set(gca,'XTick',[0:50:200],'YTick',[0.0:0.2:1.0])
set(gca,'YTickLabel',num2str(get(gca,'YTick')','%.1f'));
axis([0 200 -0.05 1.05])
xlabel('$t$','interpreter','latex');
ylabel('概率');
zhuti=title('$C_{gh}=12$');
set(zhuti,'interpreter','latex')
legend('医药研究机构({\it\fontname{Bodoni MT}r})','媒体平台({\it\fontname{Bodoni MT}m})','政府监管部门({\it\fontname{Bodoni MT}g})','群众({\it\fontname{Bodoni MT}p})');
%%%%%%%%%%%%%%%子图2
clc;clear;
a=0.7,phi=22,b=0.3,Rr=25,Fr=10,Crh=10,Crl=2,Rm=15,Fm=6,psi=12,Cml=1,Cmh=5,Cgl=1,Cgh=6,Ng=20,Tg=10,Np=12,Tp=5,Cp=10;
%figure(2)
subplot(3,1,2)
set(0,'defaultfigurecolor','w')
[t,y]=ode45(@(t,y) sifang(t,y,a,phi,b,Rr,Fr,Crh,Crl,Rm,Fm,psi,Cml,Cmh,Cgl,Cgh,Ng,Tg,Np,Tp,Cp),[0,200],[0.5,0.3,0.2,0.3]);
points=1:1:length(t);
plot(t,y(:,1),'r^-','linewidth',1,'markersize',3,'markerfacecolor','r','markerindices',points);
hold on
plot(t,y(:,2),'b-','linewidth',1);
hold on
plot(t,y(:,3),'y-.','linewidth',1);
hold on
plot(t,y(:,4),'g--','linewidth',1);
hold on
set(gca,'XTick',[0:50:200],'YTick',[0.0:0.2:1.0])
set(gca,'YTickLabel',num2str(get(gca,'YTick')','%.1f'));
axis([0 200 -0.05 1.05])
xlabel('$t$','interpreter','latex');
ylabel('概率');
zhuti=title('$C_{gh}=6$');
set(zhuti,'interpreter','latex')
legend('医药研究机构({\it\fontname{Bodoni MT}r})','媒体平台({\it\fontname{Bodoni MT}m})','政府监管部门({\it\fontname{Bodoni MT}g})','群众({\it\fontname{Bodoni MT}p})');
%%%%%%%%%%%%%%%子图3
clc;clear;
a=0.7,phi=22,b=0.3,Rr=25,Fr=10,Crh=10,Crl=2,Rm=15,Fm=6,psi=12,Cml=1,Cmh=5,Cgl=1,Cgh=2,Ng=20,Tg=10,Np=12,Tp=5,Cp=10;
%figure(3)
subplot(3,1,3)
set(0,'defaultfigurecolor','w')
[t,y]=ode45(@(t,y) sifang(t,y,a,phi,b,Rr,Fr,Crh,Crl,Rm,Fm,psi,Cml,Cmh,Cgl,Cgh,Ng,Tg,Np,Tp,Cp),[0,200],[0.5,0.3,0.2,0.3]);
points=1:1:length(t);
plot(t,y(:,1),'r^-','linewidth',1,'markersize',3,'markerfacecolor','r','markerindices',points);
hold on
plot(t,y(:,2),'b-','linewidth',1);
hold on
plot(t,y(:,3),'y-.','linewidth',1);
hold on
plot(t,y(:,4),'g--','linewidth',1);
hold on
set(gca,'XTick',[0:50:200],'YTick',[0.0:0.2:1.0])
set(gca,'YTickLabel',num2str(get(gca,'YTick')','%.1f'));
axis([0 200 -0.05 1.05])
xlabel('$t$','interpreter','latex');
ylabel('概率');
zhuti=title('$C_{gh}=2$');
set(zhuti,'interpreter','latex')
legend('医药研究机构({\it\fontname{Bodoni MT}r})','媒体平台({\it\fontname{Bodoni MT}m})','政府监管部门({\it\fontname{Bodoni MT}g})','群众({\it\fontname{Bodoni MT}p})');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%群众大图
figure(4)
%%%%%%%%%%%%%%%子图1
clc;clear;
a=0.7,phi=22,b=0.3,Rr=12,Fr=10,Crh=10,Crl=2,Rm=6,Fm=6,psi=12,Cml=1,Cmh=5,Cgl=1,Cgh=6,Ng=20,Tg=10,Np=12,Tp=5,Cp=10;
subplot(3,1,1)
set(0,'defaultfigurecolor','w')
[t,y]=ode45(@(t,y) sifang(t,y,a,phi,b,Rr,Fr,Crh,Crl,Rm,Fm,psi,Cml,Cmh,Cgl,Cgh,Ng,Tg,Np,Tp,Cp),[0,100],[0.4,0.3,0.2,0.3]);
points=1:1:length(t);
plot(t,y(:,1),'r^-','linewidth',1,'markersize',3,'markerfacecolor','r','markerindices',points);
hold on
plot(t,y(:,2),'b-','linewidth',1);
hold on
plot(t,y(:,3),'y-.','linewidth',1);
hold on
plot(t,y(:,4),'g--','linewidth',1);
hold on
set(gca,'XTick',[0:20:100],'YTick',[0.0:0.2:1.0])
set(gca,'YTickLabel',num2str(get(gca,'YTick')','%.1f'));
axis([0 100 -0.05 1.05])
xlabel('$t$','interpreter','latex');
ylabel('概率');
zhuti=title('$R_{r}=12,R_{m}=6$');
set(zhuti,'interpreter','latex')
legend('医药研究机构({\it\fontname{Bodoni MT}r})','媒体平台({\it\fontname{Bodoni MT}m})','政府监管部门({\it\fontname{Bodoni MT}g})','群众({\it\fontname{Bodoni MT}p})');
%%%%%%%%%%%%%%%子图2
clc;clear;
a=0.7,phi=22,b=0.3,Rr=25,Fr=10,Crh=10,Crl=2,Rm=15,Fm=6,psi=12,Cml=1,Cmh=5,Cgl=1,Cgh=6,Ng=20,Tg=10,Np=12,Tp=5,Cp=10;
%figure(5)
subplot(3,1,2)
set(0,'defaultfigurecolor','w')
[t,y]=ode45(@(t,y) sifang(t,y,a,phi,b,Rr,Fr,Crh,Crl,Rm,Fm,psi,Cml,Cmh,Cgl,Cgh,Ng,Tg,Np,Tp,Cp),[0,100],[0.5,0.3,0.2,0.3]);
points=1:1:length(t);
plot(t,y(:,1),'r^-','linewidth',1,'markersize',3,'markerfacecolor','r','markerindices',points);
hold on
plot(t,y(:,2),'b-','linewidth',1);
hold on
plot(t,y(:,3),'y-.','linewidth',1);
hold on
plot(t,y(:,4),'g--','linewidth',1);
hold on
set(gca,'XTick',[0:20:100],'YTick',[0.0:0.2:1.0])
set(gca,'YTickLabel',num2str(get(gca,'YTick')','%.1f'));
axis([0 100 -0.05 1.05])
xlabel('$t$','interpreter','latex');
ylabel('概率');
zhuti=title('$R_{r}=25,R_{m}=15$');
set(zhuti,'interpreter','latex')
legend('医药研究机构({\it\fontname{Bodoni MT}r})','媒体平台({\it\fontname{Bodoni MT}m})','政府监管部门({\it\fontname{Bodoni MT}g})','群众({\it\fontname{Bodoni MT}p})');
%%%%%%%%%%%%%%%子图3
clc;clear;
a=0.7,phi=22,b=0.3,Rr=50,Fr=10,Crh=10,Crl=2,Rm=35,Fm=6,psi=12,Cml=1,Cmh=5,Cgl=1,Cgh=6,Ng=20,Tg=10,Np=12,Tp=5,Cp=10;
%figure(6)
subplot(3,1,3)
set(0,'defaultfigurecolor','w')
[t,y]=ode45(@(t,y) sifang(t,y,a,phi,b,Rr,Fr,Crh,Crl,Rm,Fm,psi,Cml,Cmh,Cgl,Cgh,Ng,Tg,Np,Tp,Cp),[0,100],[0.5,0.3,0.2,0.3]);
points=1:1:length(t);
plot(t,y(:,1),'r^-','linewidth',1,'markersize',3,'markerfacecolor','r','markerindices',points);
hold on
plot(t,y(:,2),'b-','linewidth',1);
hold on
plot(t,y(:,3),'y-.','linewidth',1);
hold on
plot(t,y(:,4),'g--','linewidth',1);
hold on
%下面一行是加方框的，不过也可以在成图界面去手动制图
%rectangle('position',[0 0 2 1],'edgecolor','g','linewidth',1)
hold on
set(gca,'XTick',[0:20:100],'YTick',[0.0:0.2:1.0])
set(gca,'YTickLabel',num2str(get(gca,'YTick')','%.1f'));
axis([0 100 -0.05 1.05])
xlabel('$t$','interpreter','latex');
ylabel('概率');
zhuti=title('$R_{r}=50,R_{m}=35$');
set(zhuti,'interpreter','latex')
legend('医药研究机构({\it\fontname{Bodoni MT}r})','媒体平台({\it\fontname{Bodoni MT}m})','政府监管部门({\it\fontname{Bodoni MT}g})','群众({\it\fontname{Bodoni MT}p})');
%%%%%%%%%%%%%%%%%% 群众第三张图的小图
%下面一行决定小图的位置及其大小
axes('position',[0.35 0.2 0.1 0.1]);
clc;clear;
a=0.7,phi=22,b=0.3,Rr=50,Fr=10,Crh=10,Crl=2,Rm=35,Fm=6,psi=12,Cml=1,Cmh=5,Cgl=1,Cgh=6,Ng=20,Tg=10,Np=12,Tp=5,Cp=10;
set(0,'defaultfigurecolor','w')
[t,y]=ode45(@(t,y) sifang(t,y,a,phi,b,Rr,Fr,Crh,Crl,Rm,Fm,psi,Cml,Cmh,Cgl,Cgh,Ng,Tg,Np,Tp,Cp),[0,100],[0.5,0.3,0.2,0.3]);
points=1:1:length(t);
plot(t,y(:,1),'r^-','linewidth',1,'markersize',3,'markerfacecolor','r','markerindices',points);
hold on
plot(t,y(:,2),'b-','linewidth',1);
hold on
plot(t,y(:,3),'y-.','linewidth',1);
hold on
plot(t,y(:,4),'g--','linewidth',1);
hold on
set(gca,'XTick',[0:1:3],'YTick',[0:0.2:1])
hold on
set(gca,'YTickLabel',num2str(get(gca,'YTick')','%.1f'));
axis([0 3 -0.05 1.05])
hold on 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%民间组织大图
figure(7)
%%%%%%%%%%%%%%%子图1
clc;clear;
a=0.7,phi=22,b=0.3,Rr=25,Fr=10,Crh=10,Crl=2,Rm=15,Fm=6,psi=12,Cml=1,Cmh=5,Cgl=1,Cgh=6,Ng=20,Tg=10,Np=12,Tp=5,Cp=10;
subplot(3,1,1)
set(0,'defaultfigurecolor','w')
[t,y]=ode45(@(t,y) sifang(t,y,a,phi,b,Rr,Fr,Crh,Crl,Rm,Fm,psi,Cml,Cmh,Cgl,Cgh,Ng,Tg,Np,Tp,Cp),[0,100],[0.5,0.3,0.2,0.3]);
points=1:1:length(t);
plot(t,y(:,1),'r^-','linewidth',1,'markersize',3,'markerfacecolor','r','markerindices',points);
hold on
plot(t,y(:,2),'b-','linewidth',1);
hold on
plot(t,y(:,3),'y-.','linewidth',1);
hold on
plot(t,y(:,4),'g--','linewidth',1);
hold on
set(gca,'XTick',[0:20:100],'YTick',[0.0:0.2:1.0])
set(gca,'YTickLabel',num2str(get(gca,'YTick')','%.1f'));
axis([0 100 -0.05 1.05])
xlabel('$t$','interpreter','latex');
ylabel('概率');
zhuti=title('$\beta=0.3$');
set(zhuti,'interpreter','latex')
legend('医药研究机构({\it\fontname{Bodoni MT}r})','媒体平台({\it\fontname{Bodoni MT}m})','政府监管部门({\it\fontname{Bodoni MT}g})','群众({\it\fontname{Bodoni MT}p})');
%%%%%%%%%%%%%%%子图2
clc;clear;
a=0.7,phi=22,b=0.31,Rr=25,Fr=10,Crh=10,Crl=2,Rm=15,Fm=6,psi=12,Cml=1,Cmh=5,Cgl=1,Cgh=6,Ng=20,Tg=10,Np=12,Tp=5,Cp=10;
%figure(8)
subplot(3,1,2)
set(0,'defaultfigurecolor','w')
[t,y]=ode45(@(t,y) sifang(t,y,a,phi,b,Rr,Fr,Crh,Crl,Rm,Fm,psi,Cml,Cmh,Cgl,Cgh,Ng,Tg,Np,Tp,Cp),[0,100],[0.5,0.3,0.2,0.3]);
points=1:1:length(t);
plot(t,y(:,1),'r^-','linewidth',1,'markersize',3,'markerfacecolor','r','markerindices',points);
hold on
plot(t,y(:,2),'b-','linewidth',1);
hold on
plot(t,y(:,3),'y-.','linewidth',1);
hold on
plot(t,y(:,4),'g--','linewidth',1);
hold on
set(gca,'XTick',[0:20:100],'YTick',[0.0:0.2:1.0])
set(gca,'YTickLabel',num2str(get(gca,'YTick')','%.1f'));
axis([0 100 -0.05 1.05])
xlabel('$t$','interpreter','latex');
ylabel('概率');
zhuti=title('$\beta=0.31$');
set(zhuti,'interpreter','latex')
legend('医药研究机构({\it\fontname{Bodoni MT}r})','媒体平台({\it\fontname{Bodoni MT}m})','政府监管部门({\it\fontname{Bodoni MT}g})','群众({\it\fontname{Bodoni MT}p})');
%%%%%%%%%%%%%%%子图3
clc;clear;
a=0.7,phi=22,b=0.5,Rr=50,Fr=10,Crh=10,Crl=2,Rm=35,Fm=6,psi=12,Cml=1,Cmh=5,Cgl=1,Cgh=6,Ng=20,Tg=10,Np=12,Tp=5,Cp=10;
%figure(9)
subplot(3,1,3)
set(0,'defaultfigurecolor','w')
[t,y]=ode45(@(t,y) sifang(t,y,a,phi,b,Rr,Fr,Crh,Crl,Rm,Fm,psi,Cml,Cmh,Cgl,Cgh,Ng,Tg,Np,Tp,Cp),[0,100],[0.5,0.3,0.2,0.3]);
points=1:1:length(t);
plot(t,y(:,1),'r^-','linewidth',1,'markersize',3,'markerfacecolor','r','markerindices',points);
hold on
plot(t,y(:,2),'b-','linewidth',1);
hold on
plot(t,y(:,3),'y-.','linewidth',1);
hold on
plot(t,y(:,4),'g--','linewidth',1);
hold on
%又到了加方框的时候了
%rectangle('position',[0 0 2 1],'edgecolor','g','linewidth',1)
hold on
set(gca,'XTick',[0:20:100],'YTick',[0.0:0.2:1.0])
set(gca,'YTickLabel',num2str(get(gca,'YTick')','%.1f'));
axis([0 100 -0.05 1.05])
xlabel('$t$','interpreter','latex');
ylabel('概率');
zhuti=title('$\beta=0.5$');
set(zhuti,'interpreter','latex')
legend('医药研究机构({\it\fontname{Bodoni MT}r})','媒体平台({\it\fontname{Bodoni MT}m})','政府监管部门({\it\fontname{Bodoni MT}g})','群众({\it\fontname{Bodoni MT}p})');
%%%%%%%%%%%%%%%%%% 民间组织第三张图的小图
axes('position',[0.35 0.2 0.1 0.1]);
clc;clear;
a=0.7,phi=22,b=0.5,Rr=50,Fr=10,Crh=10,Crl=2,Rm=35,Fm=6,psi=12,Cml=1,Cmh=5,Cgl=1,Cgh=6,Ng=20,Tg=10,Np=12,Tp=5,Cp=10;
set(0,'defaultfigurecolor','w')
[t,y]=ode45(@(t,y) sifang(t,y,a,phi,b,Rr,Fr,Crh,Crl,Rm,Fm,psi,Cml,Cmh,Cgl,Cgh,Ng,Tg,Np,Tp,Cp),[0,100],[0.5,0.3,0.2,0.3]);
points=1:1:length(t);
plot(t,y(:,1),'r^-','linewidth',1,'markersize',3,'markerfacecolor','r','markerindices',points);
hold on
plot(t,y(:,2),'b-','linewidth',1);
hold on
plot(t,y(:,3),'y-.','linewidth',1);
hold on
plot(t,y(:,4),'g--','linewidth',1);
hold on
set(gca,'XTick',[0:1:3],'YTick',[0.0:0.2:1.0])
set(gca,'YTickLabel',num2str(get(gca,'YTick')','%.1f'));
axis([0 3 -0.05 1.05])
hold on 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%3d figure
%%%%%%%%%%%%%%%%%%%%%%%%%g=0的子图
clc;clear; 
a=0.7,phi=22,b=0.3,Rr=25,Fr=10,Crh=10,Crl=2,Rm=15,Fm=6,psi=12,Cml=1,Cmh=5,Cgl=1,Cgh=6,Ng=20,Tg=10,Np=12,Tp=5,Cp=10;  
figure(10)
subplot(2,1,1)
for i=0.1:0.2:1
    for j=0.1:0.2:1
        for k=0
            for l=0.1:0.2:1
            [t,y]=ode45(@(t,y) sifang(t,y,a,phi,b,Rr,Fr,Crh,Crl,Rm,Fm,psi,Cml,Cmh,Cgl,Cgh,Ng,Tg,Np,Tp,Cp),[0 50],[i j k l]);
            grid on
            %搞清楚你画的是哪三维的关系
            plot3(y(:,1),y(:,2),y(:,4),'linewidth',1);
            set(gca,'XTick',[0:0.2:1],'YTick',[0:0.2:1],'ZTick',[0:0.2:1])
            set(gca,'XTickLabel',num2str(get(gca,'XTick')','%.1f'));
            set(gca,'YTickLabel',num2str(get(gca,'YTick')','%.1f'));
            set(gca,'ZTickLabel',num2str(get(gca,'ZTick')','%.1f'));
            hold on
            axis([0 1 0 1 0 1])
            end
        end
    end
end
xlabel('$r$','interpreter','latex');
ylabel('$m$','interpreter','latex');
zlabel('$p$','interpreter','latex','Rotation',360);
title('g=0','interpreter','latex');
%%%%%%%%%%%%%%%%%%%%%%%%%g=0.8的子图
clc;clear; 
a=0.7,phi=22,b=0.3,Rr=25,Fr=10,Crh=10,Crl=2,Rm=15,Fm=6,psi=12,Cml=1,Cmh=5,Cgl=1,Cgh=6,Ng=20,Tg=10,Np=12,Tp=5,Cp=10;  
%figure(11)
subplot(2,1,2)
for i=0.1:0.2:1
    for j=0.1:0.2:1
        for k=0.8
            for l=0.1:0.2:1
            [t,y]=ode45(@(t,y) sifang(t,y,a,phi,b,Rr,Fr,Crh,Crl,Rm,Fm,psi,Cml,Cmh,Cgl,Cgh,Ng,Tg,Np,Tp,Cp),[0 50],[i j k l]);
            grid on
            %搞清楚你画的是哪三维的关系
            plot3(y(:,1),y(:,2),y(:,4),'linewidth',1);
            set(gca,'XTick',[0:0.2:1],'YTick',[0:0.2:1],'ZTick',[0:0.2:1])
            set(gca,'XTickLabel',num2str(get(gca,'XTick')','%.1f'));
            set(gca,'YTickLabel',num2str(get(gca,'YTick')','%.1f'));
            set(gca,'ZTickLabel',num2str(get(gca,'ZTick')','%.1f'));
            hold on
            axis([0 1 0 1 0 1])
            end
        end
    end
end
xlabel('$r$','interpreter','latex');
ylabel('$m$','interpreter','latex');
zlabel('$p$','interpreter','latex','Rotation',360);
title('{\it\fontname{Bodoni MT}g}=0.8');
    image.png
    image.png image.png image.png

    出图效果还是有差距，在编辑器中修改。


    image.png
    image.png
    image.png
    image.png

    相关文章

      网友评论

          本文标题：4. 四方演化博弈代码

          本文链接：https://www.haomeiwen.com/subject/fxdldltx.html

          延伸阅读

          深度阅读

          • 您也可以注册成为美文阅读网的作者，发表您的原创作品、分享您的心情！

          栏目导航

          热点阅读

          关于我们|服务条款|联系我们|4. 四方演化博弈代码|投稿指南|网站地图|RSS订阅|排版工具|手机版

          提供经典美文摘抄,优美散文欣赏,现代诗歌精选,短篇小说,心情随笔,表白情书范文,故事会在线阅读欣赏

          Copyright © 2014-2023 Haomeiwen.com All Rights Reserved. 好美文阅读网 版权所有

          备案信息：桂公网安备 45052102000051号 · 桂ICP备13007215号-3

          本站所收录作品、热点评论等信息部分来源互联网，目的只是为了系统归纳学习和传递资讯

          所有作品版权归原创作者所有，与本站立场无关，如不慎侵犯了你的权益，请联系我们告知，我们将做删除处理！