公共財ゲーム関連

scilabのプログラム例

ポピュレーションダイナミクス＋公共財進化ゲーム

Hauertらのモデルの軌跡描画プログラム例

funcprot(0);

function ap = avg_payoff(x,y) // 報酬の計算、引数：協力、非協力の割合 x,y
    ap = zeros(2,1);
    z = 1.0 -x - y;
    ap(2) = r*x/(x+y) * (1.0 - (1 - z**n)/(1-z)/n);
    ap(1) = ap(2) - (1.0 + (r-1.0)*z**(n-1) - r*(1-z**n)/n/(1-z));
endfunction

function xdot=f(t,x,payoff_func) // リプリケータダイナミクスの式
    z = 1.0 - x(1)-x(2);
    d = 0.5;
    payoff = payoff_func(x(1),x(2));
    t_payoff = (x(1)*payoff(1) + x(2)*payoff(2));
    a_payoff = t_payoff/(x(1)+x(2));
    xdot=zeros(2,1);
    xdot(1) = (z*payoff(1)-d)*x(1);
    xdot(2) = (z*payoff(2)-d)*x(2);
endfunction

//ここからメイン

r=input('r=');
n=input('N=');
x0=[0.24;0.01] // 初期値
t=linspace(0,2000.,10000);
y=ode(x0,0, t,list(f,avg_payoff)); // リプリケータダイナミクス方程式を解く

//　出力のための変数変換
ratio_x=zeros(length(y(1,:)),1);
xy = ratio_x;
for i=1:length(y(1,:))
    ratio_x(i) = y(1,i)/(y(1,i)+y(2,i));
    xy(i) = y(1,i)+y(2,i);
end

plot2d(xy,ratio_x)

軌跡をアニメーションとして表示する例

function ap = avg_payoff(x,y)
    ap = zeros(2,1);
    z = 1.0 -x - y;
    ap(2) = r*x/(x+y) * (1.0 - (1 - z**n)/(1-z)/n);
    ap(1) = ap(2) - (1.0 + (r-1.0)*z**(n-1) - r*(1-z**n)/n/(1-z));
endfunction

function xdot=f(t,x,payoff_func)
    z = 1.0 - x(1)-x(2);
    d = 0.5;
    payoff = payoff_func(x(1),x(2));
   // a_payoff = (x(1)*payoff(1) + x(2)*payoff(2))/(x(1)+x(2));
  xdot=zeros(2,1);
  xdot(1) = (z*payoff(1)-d)*x(1);
  xdot(2) = (z*payoff(2)-d)*x(2);
endfunction

//ここからメイン
funcprot(0);
r=input('r=');
n=input('N=');
x_plus_y = input('x+y=');
ratio_x = input('x/(x+y)=');

x_init = x_plus_y*ratio_x;
x0=[x_init; x_plus_y - x_init]; // 初期値
t=linspace(0,2000.,10000);

y=ode(x0,0, t,list(f,avg_payoff)); // 微分方程式の計算

//　出力のための変数変換
xv=zeros(y);
for i=1:length(y(1,:))
    xv(2,i) = y(1,i)/(y(1,i)+y(2,i));
    xv(1,i) = y(1,i)+y(2,i);
  //  printf("%f %f : ",ratio_x(i),xy(i));
end

// アニメーション

my_handle = scf(100001); // 引数は整数で図のid 返り値は図のハンドル
clf(my_handle,"reset");
title("evolution game model");

function h = poly2d(x,y)  // 2点間をグラフに書かせる関数定義
    xpoly(x,y); h = gce();
endfunction

curAxe = gca();
drawlater();
curAxe.title.font_size =3
curAxe.data_bounds = [min(xv,'c')';max(xv,'c')']
curAxe.x_label.text = 'x+y'
curAxe.x_label.font_size = 3
curAxe.y_label.text = 'x/(x+y)'
curAxe.y_label.font_size = 3
curAxe.axes_visible = 'on'
curAxe.labels_font_size = 3

//the trace
p = poly2d(xv(1,1),xv(2,1));
drawnow()

// 線を次々書かせる繰り返し
for k=1:size(xv,2)
   if ~is_handle_valid(my_handle) then
      break;
    end

    if is_handle_valid(p) then
      p.data=[p.data;
        xv(1:2,k)'];
    else
      break;
    end
end