社会を＜モデル＞でみる：数理社会学への招待
10章：なぜ好きなのにケンカするのか
x:xさんの行為、正の値：２人にとって良い行為;負の値:２人もしくはyさんにとって悪い行為
y:yさんの行為、正の値：２人にとって良い行為;負の値:２人もしくはxさんにとって悪い行為
x(t+1)=x(t)+a*x(t)+b*y(t)
y(t+1)=y(t)+c*x(t)+d*y(t)

# なぜ好きなのにけんかをするのか
fighting_couple<-function(a,b,c,d) {
  timeSep=0.05;ts=seq(1,50,timeSep);n_ts=length(ts)
  x=matrix(0,nrow=n_ts,ncol=1)
  y=matrix(0,nrow=n_ts,,ncol=1)
  initX=c(rep(-40,5),rep(-20,5),rep(0,5),rep(20,5),rep(40,5))
  initY=c(rep(c(-40,-20,0,20,40),5))
  initX=initX[-13]
  initY=initY[-13]
  lengthINI=length(initX)
  for (i_ini in 1:lengthINI) {  
    x[1]=initX[i_ini];y[1]=initY[i_ini];
    for (i_gen in 2:n_ts) {
       x[i_gen]=x[i_gen-1]+(a*x[i_gen-1]+b*y[i_gen-1])*timeSep
       y[i_gen]=y[i_gen-1]+(c*x[i_gen-1]+d*y[i_gen-1])*timeSep
    }  
    if (i_ini==1) {  
       plot(x,y,xlim=c(-50,50),ylim=c(-50,50),col=4,type='l',lwd=2,
         xlab="X's Action",ylab="Y's Action")
       arrows(x[2],y[2],x[3],y[3],col=4,lwd=2,length=0.15)} else {  
       lines(x, y, col=4, lwd=2)
       arrows(x[2], y[2], x[3], y[3], col=4,lwd=2, length=0.15)
    }  
  }  
}  

par(mfrow=c(2,3))
fighting_couple(-1,0.0,0.5,-1)
fighting_couple(-1,2,-1,-1)
fighting_couple(0,-1,1,0)
fighting_couple(1,-2,2,0)
fighting_couple(1,0,0.5,1)
fighting_couple(1,-4,-4,0)

fun1<-function(a,b,c,d,N,initX,initY) {
  dt=0.05;
  x=matrix(0,nrow=N,ncol=1);y=matrix(0,nrow=N,ncol=1);
  x[1]=initX;y[1]=initY;
  for (i_rep in 2:N) {
    x[i_rep]=x[i_rep-1]+(a*x[i_rep-1]+b*y[i_rep-1])*dt
    y[i_rep]=y[i_rep-1]+(c*x[i_rep-1]+d*y[i_rep-1])*dt
  }
  return(data.frame(x,y))
}
fun2<-function(a,b,c,d,N,initX, initY) { 
  res<-fun1(a,b,c,d,N,initX[1],initY[1]);
  nXYs=length(initX)
  plot(res$x,res$y,xlim=c(-50,100),ylim=c(-50,100),type='l')
  for (i_loop in 2:nXYs){
    res<-fun1(a,b,c,d,N,initX[i_loop],initY[i_loop]);
    lines(res$x,res$y,type='l')
  }
}
# 実行例
initX=rep(seq(10,50,10),5)
initY=sort(rep(seq(10,50,10),5))
fun2(0,-1,1,0,500,initX,initY)