楊波 范敏 劉文奇 陳曉松

1)(昆明理工大學數(shù)據科學研究中心,昆明 650500)

2)(昆明理工大學理學院,昆明 650500)

3)(中國科學院理論物理研究所,理論物理前沿重點實驗室,北京 100190)

4)(中國科學院大學物理科學學院,北京 100049)

自我質疑機制下公共物品博弈模型的相變特性?

楊波1)2)?范敏1)2)劉文奇1)2)陳曉松3)4)

1)(昆明理工大學數(shù)據科學研究中心,昆明 650500)

2)(昆明理工大學理學院,昆明 650500)

3)(中國科學院理論物理研究所,理論物理前沿重點實驗室,北京 100190)

4)(中國科學院大學物理科學學院,北京 100049)

公共物品博弈是研究群體相互作用的經典模型,廣泛用于解釋自私個體間合作的涌現(xiàn)和保持.本文從理論分析和蒙特卡羅模擬兩個方面研究了二維正方格子上一個有償懲罰機制下隨自我質疑更新規(guī)則演化的公共物品博弈模型的相變特性.理論分析方面,將公共物品博弈模型轉化為一個外場不為零的鐵磁Ising模型.通過有效能量發(fā)現(xiàn):不存在懲罰時,個體間的耦合強度為零,體系只有外場作用;存在懲罰時,個體間包含最近鄰、次近鄰和第三近鄰相互作用且外場不為零.蒙特卡羅模擬方面,首先驗證了理論分析的正確性,然后對公共物品博弈模型相關的一級相變和二級相變進行了有限尺度標度分析.研究發(fā)現(xiàn):1)蒙特卡羅模擬所得結果與類Ising模型分析結果完全吻合;2)相比二維Ising模型,公共物品博弈的二級相變臨界指數(shù)發(fā)生了變化;3)公共物品博弈的一級相變與二維Ising模型相同.

Ising模型,有限尺度標度理論,蒙特卡羅模擬,自我質疑更新規(guī)則

1 引 言

演化博弈理論(evolutionary game theory)突破了傳統(tǒng)博弈論關于“理性人”和“完全信息”的限制,它強調有限理性的博弈個體在重復博弈過程中根據自身掌握的局部信息,通過自適應學習做出盡可能占優(yōu)的策略［1,2］.在人類社會和自然界中,自私個體之間產生合作是一個“驚人”的現(xiàn)象,得到了許多學者的重視和研究［3,4］.空間演化博弈理論(將多個個體放置在社會網絡上進行重復博弈)的提出從空間結構和社會網絡的角度很好地解釋了合作的涌現(xiàn),同時也開啟了演化博弈研究的新篇章［1,2,5,6］.

公共物品博弈(public goods game)是一種典型的多人空間演化博弈模型,在博弈過程中,N個參與者獨立決定向公共池子投資(合作)或不投資(不合作),初始投資總額通過增益系數(shù)放大r(1＜r＜N)倍后,平均分給群組內所有個體,無論初始投資與否.顯然,采用不合作策略可以在事先不進行任何投資的情況下獲得收益,搭上合作者的便車.公共物品博弈生動再現(xiàn)了個人與集體間的矛盾,有著廣泛的應用,如:在班級大掃除中,每個學生付出的勞動是不對等的,但干凈明亮的學習環(huán)境卻是大家一起共享的.盡管所有成員都合作時,群組可以實現(xiàn)利益最大化,然而對一個自私且理性的個體而言,最優(yōu)策略是選擇背叛.從而在個體最優(yōu)策略與群組最優(yōu)策略之間形成了社會困境,即所謂的公共品悲劇(the tragedy of the commons)［1,7?9］.然而,現(xiàn)實生活中這種悲劇并未出現(xiàn).為了解釋合作的廣泛傳播與個體的自私性之間的矛盾,許多機制被相繼提出,如:期望誘導重連機制、刪邊機制、懲罰機制、獎勵機制、社會多樣性機制、馬太效應、志愿者參與機制、遷徙機制、調整投資機制、自適應有限投資反饋機制和條件策略［10,11］等.本文基于前人對懲罰機制的研究［12?18］,采用一類簡化的有償懲罰機制作用下的公共物品博弈模型進行研究.

空間演化博弈的兩個研究重點是策略更新規(guī)則和網絡的空間結構［1］.網絡演化博弈的策略更新規(guī)則可以劃分為兩類:學習機制（模仿）和自我質疑機制（自?。?學習機制主要是向網絡中的最近鄰學習,包括學習鄰居中收益最大的［5］;以一定的概率學習那些收益比自己高的個體們［19,20］;任意選擇最近鄰中的一個,比較收益差,以較大的概率模仿收益比自己高的個體［21?23］.自我質疑機制的研究相對較少［24?27］,它是指任意博弈個體先計算當前策略所得總收益,然后采用自身反策略進行一次虛擬博弈,并計算虛擬總收益,通過比較當前博弈和虛擬博弈收益變化的多少,來決定自己下一局所采用的策略.自我質疑機制與統(tǒng)計物理中Ising模型的單自旋翻轉蒙特卡羅模擬方法十分類似.基于此,本文研究自我質疑更新規(guī)則下的公共物品博弈模型與Ising模型的等價關系［28］.之前,演化博弈模型的相關研究主要依靠計算機模擬方法,本文提出的類Ising方法將為該領域的研究提供理論基礎.值得一提的是,經濟學中廣泛研究的隨機最優(yōu)反應均衡(quantal response equilibrium,QRE)與自我質疑機制十分類似,它用隨機反應代替?zhèn)鹘y(tǒng)的最優(yōu)反應,將納什均衡作為它的一個特例進行處理［29?31］.

本文提出一個有償懲罰機制下按自我質疑機制演化的公共物品博弈模型,從以下三個方面展開研究.理論分析方面,找到與該演化博弈模型等價的Ising模型并運用相關理論進行分析預測;蒙特卡羅模擬方面,驗證理論分析的結果;相變特性方面,任選一組懲罰參數(shù),研究博弈模型中伴隨的相變現(xiàn)象,分析相變的種類,測定相關臨界指數(shù).

2 自我質疑更新規(guī)則下的公共物品博弈模型

2.1 模型的引入

在由N個個體參與的公共物品博弈中,每個個體都有兩種策略可供選擇,投資或不投資,假設原始投資額為1.投資總額放大r倍后(可認為是投資后獲利),平均分給群組內所有參與博弈的個體,而不管個體是否向公共池內投過資,其中r為強化因子(enhancement factor),滿足1＜r＜N.假設群組內有N+個合作者,則投資者的收益為gC=rN+/N?1,而不投資者的收益為gD=rN+/N.顯然,不投資者的收益永遠大于投資者的收益.所以,不投資者可以在投資者投資的基礎上“搭便車”.從而,公共物品博弈的納什均衡是,所有人都不投資.為了解釋自私個體間仍有合作行為產生,本文提出一種簡化的有償懲罰機制,即懲罰的代價為C(cost),被懲罰者的處罰額度為F( fi ne).此時,投資者的收益為gC=rN+/N?1?C(N?N+),不投資者獲得收益gD=rN+/N?FN+.

采用一個固定的網絡來刻畫個體之間的博弈關系.在網絡模型確定之后,種群中的個體占據網絡中的節(jié)點.每個個體都是純策略者并只能選擇合作或者不合作策略.每個個體與自己的鄰居組成一個群體進行公共品博弈.在每一輪博弈中,任選一個個體i,然后根據自我質疑機制決定是否改變當前策略.即選定個體i后,個體i采用當前策略與周圍鄰居博弈獲得當前收益Gi,采用當前策略的反策略與周圍鄰居進行一次虛擬博弈,獲得虛擬收益?zhèn)€體i通過比較當前收益和虛擬收益,以一定概率決定是否改變當前策略,概率的選擇為

2.2 公共物品博弈模型與Ising模型間的轉化關系

在空間網絡中,任意個體i,不僅與周圍最近鄰個體博弈,還與次近鄰進行博弈.即個體i不僅參與以自己為中心節(jié)點進行的群體博弈,而且參與以鄰居為中心節(jié)點的群體博弈,即總共經歷ki+1次博弈,其中ki為節(jié)點i的度.當個體i參與以自己為中心節(jié)點的群體博弈時,參與博弈的個體總數(shù)為ki+1,當參與鄰居j為中心節(jié)點的群體博弈時,參與人數(shù)為kj+1.圖1顯示了二維格點上個體i參與的五次公共物品博弈,深色格子表示博弈過程中的中心節(jié)點,淺色格子表示中心節(jié)點的最近鄰節(jié)點,每輪博弈在由深色和淺色節(jié)點組成的群體中進行,白色格子表示其余節(jié)點.

圖1 (網刊彩色)二維正方格子上任意個體i每輪博弈中參與以不同節(jié)點為中心的五次公共物品博弈Fig.1.(color online)The public goods games which involve player i on two-dimension square lattice each round.

個體i每一輪博弈的收益來自于ki+1個群體博弈所得收益的總和,為便于推導,用Si表示個體i的策略,有兩種取值,Si=+1表示合作,Si=?1表示不合作,在以i為中心結點的群體博弈中,當Si=+1時,

其中ni+表∑示i的鄰居中采用合作策略的個體總和,ni+=j(1+Sj/2);ni?表示∑i的鄰居中采用不合作∑策略的個體總和,ni?=j(1?Sj)/2=(ki?jSj)/2,j是i的最近鄰.當Si=?1時,

個體i翻轉所帶來的收益差為

將ni+和ni?的表達式代入(4)式,并用新記號?gii描述個體i在以i為中心節(jié)點進行博弈的過程中所帶來的收益變化,

同理可得,?gij,個體i在以j為中心節(jié)點進行博弈時,改變策略帶來的收益差為

其中,j是i的最近鄰,表示對i的最近鄰求和;表示在以j為中心的博弈群體中,對除i以外的節(jié)點進行求和.從而,一輪博弈中i改變策略帶來的總收益差?Gi為

將(5)式和(6)式代入(7)式,并整理可得

從上面的表達式可以看出,總收益的改變由兩部分組成:前面部分只與博弈個體i所處的狀態(tài)和空間結構有關,后面部分則包含個體i與最近鄰j以及j的鄰居k間的相互作用.顯然,當不存在懲罰時,總收益差只與個體i和網絡結構有關.也就是說,懲罰的引入直接導致個體間產生相互作用.為了便于研究和分析,本文以二維正方格子上的公共物品博弈為例,展開進一步討論.在二維方格子中,ki=kj=4且∑j=4,表達式(8)中第一項可化簡為Si［5?r?10(F?C)］,第二項中求和項是對所有與個體i相關的ki+1個集團中所有的參與者進行求和.圖2是對所有求和節(jié)點的求和次數(shù)的示意圖,可以將這樣的求和規(guī)則改寫成:兩次最近鄰求和加上兩次次近鄰求和再加上一次第三近鄰求和,即其中j表示最近鄰,j′表示次近鄰,j′′表示第三近鄰.代入到(8)式的第二項可得:

因此,表達式(8)可寫成

公共物品博弈中追求利益最大化,而Ising模型則希望系統(tǒng)能量達到最小.所以,定義有效能量差為?Ei=??Gi,相應的,

考慮最近鄰、次近鄰和第三近鄰的Ising模型的表達式為

其中J表示最近鄰相互作用強度,J′表示次近鄰相互作用強度,J′′表示第三近鄰相互作用強度,〈ij〉表示對最近鄰求和,〈ij′〉表示對次近鄰求和,〈ij′′〉表示對第三近鄰求和,H表示外場強度.

圖2 (網刊彩色)任意個體i參與的公共物品博弈等價于與最近次、次近鄰分別進行兩次兩兩博弈和與第三近鄰進行一次兩兩博弈Fig.2.(color online)The number of interaction between player i and it’s neighbours.

Ising模型的局部能量變化關系?E=E′?E可表示為

對比表達式(10)和表達式(12)可得公共物品博弈對應類Ising模型的耦合強度和外場為:

顯然,耦合強度是懲罰費用和懲罰代價的函數(shù),因為F＞0且C＞0,所以耦合強度J＞0.由Ising模型相關理論可知,該模型為鐵磁Ising模型.此外,耦合強度的大小與空間結構和強化因子無關,只與懲罰的取值相關.外場由兩部分構成:前一項表示一輪博弈中,均分后的利益是否大于原始投資;后一項來源于懲罰的引入,只要懲罰的費用大于懲罰的代價,外場都為正,即促進合作的產生.特別指出,運用矩陣分解法也可得到相似的等價關系［2］.

綜合寫出有效能量為

2.3 蒙特卡羅模擬

本部分采用Metropolis算法對周期邊界條件下二維空間正方格子上的公共物品博弈模型進行蒙特卡羅模擬.空間格子的大小為N=L×L,用偽隨機數(shù)生成任意初始構型后按照自我質疑機制進行演化.首先,隨機選擇群體中的任意博弈個體i,然后用自我質疑機制決定個體i是否改變當前策略,重復以上操作,直至系統(tǒng)趨于穩(wěn)定狀態(tài)為止.為了保證計算的準確性,扔掉前10000個蒙特卡羅模步(N次隨機試驗稱為一個蒙特卡羅步).計算接下來的2×105個蒙特卡羅步的平均值.為消除初始條件對系統(tǒng)演化的影響,再選擇100個不同的初始構型進行系綜平均.

在演化博弈模型的社會學研究中,通常選擇合作者占比f作為觀測量,其表達式為

其中N代表參與博弈的總人數(shù),N+代表總人數(shù)中采用合作策略的人數(shù).

由表達式(17)可知,當不存在懲罰時,耦合強度為零,體系的有效能量為

顯然,當r＞5時為正向外場,整個體系趨向于合作,而r＜5時為負向外場,體系趨向于不合作.如圖3所示,合作與不合作以r=5為分界線,溫度較低時,在強化因子逐漸增大的過程中,體系在r=5處,從完全不合作跳變到完全合作,隨著溫度的升高,這種跳變逐漸變得平緩.

圖3 (網刊彩色)不存在懲罰時合作頻率隨強化因子的變化Fig.3.(color online)Frequency of cooperators vs.r for various T.

當存在懲罰時,體系表現(xiàn)為鐵磁相互作用,系統(tǒng)狀態(tài)隨外場而變.顯然,臨界外場Hc=0,即1?rc/5?2(F?C)=0.解得rc=5?10(F?C),則當r＞rc時,系統(tǒng)為完全合作狀態(tài),r＜rc時系統(tǒng)為完全不合作狀態(tài).圖4顯示不同的懲罰額度F和懲罰的代價C下,合作頻率f隨強化因子r的變化,固定溫度T=0.1.顯然,蒙特卡羅模擬的結果與Ising模型理論分析的結果很好地吻合.合作頻率以rc為突變點發(fā)生跳變.從社會學角度看,懲罰的引入增強了個體間的聯(lián)系,相應地也增強了體系對不確定性因素T的魯棒性.

圖4 (網刊彩色)不同懲罰條件下合作頻率隨強化因子的變化Fig.4.(color online)Frequency of cooperators vs.r at T=0.1 for various punishment.

3 公共物品博弈模型的相變特性

從與公共物品博弈模型對應的Ising模型中看到,耦合強度J恒大于零,體系表現(xiàn)為鐵磁相互作用.鐵磁Ising的相變和臨界現(xiàn)象在統(tǒng)計物理中已經廣泛研究,基于判斷有效能量推導的正確性和解釋公共物品博弈誘發(fā)的相變現(xiàn)象兩個原因,我們采用蒙特卡羅模擬和有限尺度標度理論詳細分析公共物品博弈模型的相變特性.在Ising模型中,零場下隨著溫度的變化,存在鐵磁序到順磁序的二級相變(即由完全合作或完全不合作態(tài)向混亂狀態(tài)的轉變),相變溫度稱為臨界溫度,在臨界溫度附近表現(xiàn)出標度特性;外場非零時,只要溫度小于臨界溫度,隨著外場的變化會引發(fā)一級相變(由完全合作態(tài)到完全不合作態(tài)的轉變或其逆過程).為了便于與Ising模型進行對比,本部分仍采用Ising模型的相關觀測量展開研究,即選擇磁化強度(magnetization)作為序參量.其表達式為

磁化強度與合作者占比的關系為m=2f?1.此外,定義磁化率為

一級相變和二級相變具有完全不同的相變特性,如二級相變由于關聯(lián)長度的發(fā)散引發(fā)一系列奇異特性而一級相變則表現(xiàn)出δ函數(shù)奇異性.為了討論方便,先簡單介紹一級相變和二級相變的標度函數(shù),然后再進行蒙特卡羅模擬.

二級相變的相變點通常稱為臨界點,在臨界點附近,系統(tǒng)的關聯(lián)長度發(fā)散,從而引起系統(tǒng)相關熱力學量的奇異性,如比熱和磁化率在臨界點處發(fā)散.通常無窮大系統(tǒng)才表現(xiàn)出臨界現(xiàn)象,然而無論是實際系統(tǒng)還是計算機模擬都不可能無窮大.Fisher于20世紀70年代提出的有限尺度標度理論通過對小系統(tǒng)的計算機模擬實現(xiàn)了對臨界現(xiàn)象的研究［32］.具體而言,在臨界點附近磁化強度和磁化率存在以下標度關系:

其中,fm和fχ為普適函數(shù),t=(T?Tc)/Tc和h=(H?Hc)/Hc分別稱為約化溫度和約化外場,L是有限系統(tǒng)的特征長度,ν和νh分別代表關聯(lián)長度對溫度和外場發(fā)散的臨界指數(shù),1/νh=(γ+β)/ν,β和γ是描述熱力學函數(shù)m和χ臨界行為的臨界指數(shù).為確定臨界點的位置,定義四階Binder累積量為

Binder累計量的標度形式滿足

其中U為普適函數(shù),在臨界點處,它與系統(tǒng)尺寸無關,即不同尺度下U的曲線交于臨界點.

在無窮大系統(tǒng)中,一級相變的關聯(lián)長度不發(fā)散且相變表現(xiàn)為δ函數(shù)奇異性,即磁化率可用一個δ函數(shù)表示.但在有限系統(tǒng)模擬中,δ函數(shù)的奇異性將被平滑掉.文獻［33—37］運用熱力學漲落理論獲得了一級相變磁化強度、磁化率和四階累積量的有限尺度標度函數(shù).磁化強度為

其中Msp是熱力學極限下的自發(fā)磁化強度,χ0是系統(tǒng)處于單相時的磁化率,d為系統(tǒng)維度,kB表示玻爾茲曼常數(shù).

磁化率為

上式表明,在L→∞極限下發(fā)生在H=0處的δ函數(shù)奇異性,當L有限時被平滑為一個峰,其高度正比于Ld,其寬度正比于L?d.

四階累計量的表達式為

從磁化強度、磁化率和四階累計量的表達式看出,空間維度d是一級相變的唯一標度指數(shù).

3.1 外場為零時的相變特性

公共物品博弈對應的有效能量包括:最近鄰、次近鄰和第三近鄰相互作用,精確解未知.外場為零時,體系隨溫度變化表現(xiàn)為二級相變,通過四階Binder累積量可以確定臨界溫度的位置,通過有限尺度標度理論可以獲得相應的臨界指數(shù).模擬過程中,固定懲罰額度F=0.2,懲罰費用C=0.1,此時外場h=(r?4)/2.顯然,當r＞4時,h＞0為正向外場,r＜4時,為負向外場,r=4時,外場為零.為研究的方便,先研究外場為零時系統(tǒng)隨溫度變化表現(xiàn)出的二級相變.

圖5和圖6是零外場時磁化強度和磁化率隨溫度變化的曲線,顯然,不同尺寸下存在顯著的有限尺度效應.磁化強度隨著系統(tǒng)尺寸的增大,臨界點附近的曲線變得越來越陡,越來越接近系統(tǒng)無窮大時的變化規(guī)律.熱力學極限下,磁化率在臨界點處發(fā)散,但在有限系統(tǒng)下為有限值,圖6正好說明了這一點,隨著系統(tǒng)尺度的不斷增大,曲線的峰值變得越來越大.

圖5 (網刊彩色)外場為零時不同尺度下磁化強度隨溫度的變化Fig.5.(color online)Magnetization m plotted vs temperature T at h=0 for various L.

圖6 (網刊彩色)外場為零時不同尺度下磁化率隨溫度的變化Fig.6.(color online)Susceptibility χLplotted vs temperature T at h=0 for various L.

不同尺寸下的Binder累積量曲線相交于臨界溫度. 如圖7所示,可確定臨界溫度Tc=1.065±0.0067.通過有限尺度標度理論,可以測定臨界指數(shù)β/ν=0.154±0.0042,γ/ν=1.71±0.023,1/ν=1.085±0.027.將模擬結果與二維Ising模型的精確值對比,發(fā)現(xiàn)二者的臨界指數(shù)存在一定差異,有可能是普適類發(fā)生了變化.為了進一步驗證我們的設想,圖8和圖9用二維Ising模型的臨界指數(shù)去坍塌模擬所得數(shù)據.可以明顯看到臨界點右半邊曲線的坍塌效果較差,導致這一現(xiàn)象的原因是博弈模型中不僅包含最近鄰相互作用,還包括次近鄰和第三近鄰相互作用.

圖7 (網刊彩色)外場為零時不同系統(tǒng)尺度下四階累積量隨溫度的變化Fig.7.(color online)Reduced cumulant U plotted vs temperature T at h=0 for various L.

圖8 (網刊彩色)外場為零時臨界溫度附近磁化強度的有限尺度標度關系Fig.8. (color online)Scaled magnetization mL1/8 plotted vs scaled temperature εL1/1at T=1.065 and various L.

圖9 (網刊彩色)外場為零時臨界溫度附近磁化率的有限尺度標度關系Fig.9.(color online)Scaled susceptibility χL/L7/4 plotted vs scaled temperature εL1/1at T=1.065 and various L.

3.2 外場非零時的相變特性

當溫度小于臨界溫度時,沿外場演化模型表現(xiàn)為一級相變.當溫度等于臨界溫度時,表現(xiàn)為二級相變.為更好地對比和區(qū)分兩種相變,分別討論T=1(T＜Tc,一級相變)和T=1.065(T=Tc,二級相變)兩種情況.

圖10和圖11分別為T=1和T=1.065時磁化強度隨外場的變化曲線.當系統(tǒng)無窮大時,在相變點r=5處,磁化強度不連續(xù),發(fā)生跳變.當系統(tǒng)尺寸有限時,磁化強度曲線的奇異性被平滑掉.隨著系統(tǒng)尺寸的不斷增大,曲線變得越來越陡峭,系統(tǒng)無窮大時變得不連續(xù).對比兩圖發(fā)現(xiàn):在rc附近,T=Tc處磁化強度的變化比T＜Tc處的變化要緩慢一些.

圖10 (網刊彩色)T＜Tc時不同尺度下磁化強度隨強化因子的變化Fig.10.(color online)Magnetization m plotted vs the factor r at T=1 for various L.

圖11 (網刊彩色)T=Tc時不同尺度下磁化強度隨強化因子的變化Fig.11.(color online)Magnetization m plotted vs the factor r at T=1.065 for various L.

圖12和圖13分別展示了T=1和T=1.065時磁化率隨外場的變化曲線.從(27)式知道,磁化率在相變點處隨系統(tǒng)尺寸冪率發(fā)散,在兩圖中可以明顯看到,隨著系統(tǒng)尺度的增大,磁化率曲線的尖峰不斷增大,當系統(tǒng)尺度為無窮大時發(fā)散,T=Tc處的峰值沒有T＜Tc處的峰值高.這是由于系統(tǒng)無窮大時,一級相變在臨界點附近δ函數(shù)發(fā)散而二級相變則按臨界指數(shù)γ冪率發(fā)散.

圖12 (網刊彩色)T＜Tc時不同尺度下磁化率隨強化因子的變化Fig.12.(color online)Susceptibility χLplotted vs the factor r at T=1 for various L.

圖13 (網刊彩色)T=Tc時不同尺度下磁化率隨強化因子的變化Fig.13.(color online)Susceptibility χLplotted vs the factor r at T=1.065 for various L.

圖14 (網刊彩色)T＜Tc時不同尺度下四階累積量隨強化因子的變化Fig.14.(color online)Reduced cumulant U plotted vs the factor r at T=1 for various L.

圖14和圖15得到了公共物品博弈相變點的位置,可明顯看到不同尺度下所有曲線都交于相變點r=4處.

圖15 (網刊彩色)T=Tc時不同尺度下四階累積量隨強化因子的變化Fig.15.(color online)Reduced cumulant U plotted vs the factor r at T=1.065 for various L.

圖16 (網刊彩色)T＜Tc時臨界外場處磁化率隨尺度變化的雙對數(shù)圖Fig.16.(color online)Log-log plot of susceptibility maximum vs linear dimension at T=1.

圖16和圖17為磁化率隨系統(tǒng)變化的雙對數(shù)曲線,擬合曲線的斜率可以獲得磁化率隨尺度變化的臨界指數(shù),當T＜Tc時臨界指數(shù)d=1.935±0.0093;當T=Tc時,臨界指數(shù)γ/ν=1.714±0.0026.二維Ising模型臨界指數(shù)的精確結果為d=2和γ/ν=1.75.顯然二者間仍然有微小的差別,為了進一步驗證模擬結果的正確性,圖18和圖19用二維Ising模型的臨界指數(shù)對模擬所得數(shù)據進行坍塌,其中ε=(r?rc)/rc為約化增益系數(shù),γ/ν的理論值為7/4,(γ+β)/ν的理論值為15/8.圖16基于表達式(27)進行數(shù)據坍塌,在坍塌過程中忽略修正項χ0.從圖中可以看出,公共物品博弈模型隨外場變化的相變特性和二維Ising模型完全相同.

圖17 (網刊彩色)T=Tc時臨界外場處磁化率隨尺度變化的雙對數(shù)圖Fig.17.(color online)Log-log plot of susceptibility maximum vs linear dimension at T=1.065.

圖18 (網刊彩色)T＜Tc時臨界外場附近磁化率的有限尺度標度關系Fig.18. (color online)Scaled susceptibility χL/L2 plotted vs scaled fi eld εL2at T=1 and various L.

圖19 (網刊彩色)T=Tc時臨界外場附近磁化率的有限尺度標度關系Fig.19.(color online)Scaled susceptibility χL/L7/4plotted vs scaled fi eld εL15/8at T=1.065 and various L.

圖20和圖21用臨界指數(shù)的精確解對模擬得到的磁化強度的數(shù)據進行了坍塌.T＜Tc時,圖20顯示,在臨界點rc附近曲線的坍塌效果較好,隨著r逐漸遠離rc,坍塌效果變差.其原因來源于表達式(26),在rc附近m≈hMspL2kBT.圖21顯示了較好的坍塌效果.

圖20 (網刊彩色)T＜Tc時臨界外場附近磁化強度的有限尺度標度關系Fig.20.(color online)Scaled magnetization m plotted vs scaled fi eld εL2at T=1 and various L.

圖21 (網刊彩色)T=Tc時臨界外場附近磁化強度的有限尺度標度關系Fig.21.(color online)Scaled magnetization mL1/8 plotted vs scaled fi eld εL15/8at T=1.065 and various L.

4 結 論

本文研究了有償懲罰機制下隨自我質疑更新規(guī)則演化的公共物品博弈模型.首先采用博弈模型收益差類比Ising模型能量差的方法,獲得了二維正方格子上博弈模型的有效能量.有效能量顯示:沒有懲罰時,個體間不存在相互作用,合作策略的選取只與外場有關,外場則取決于原始投資和最終獲得收益分紅間的關系;存在懲罰時,個體間包括最近鄰、次近鄰和第三近鄰相互作用,同時也減小了外場的作用(懲罰需要付出代價).也就是說,懲罰的引入增強了個體間的關聯(lián),促使博弈群體具有更強的魯棒性,微小的擾動很難引起個體策略的改變.其次,計算機模擬博弈過程獲得了不同理性程度和不同懲罰額度下合作者占比的變化曲線.結果顯示:模擬結果和類Ising分析結果完全吻合.最后,以Ising模型相關的熱力學量為基礎研究了公共物品博弈模型的相變和臨界現(xiàn)象.分別對二級相變和一級相變進行了討論,結果顯示,公共物品博弈模型隨理性參數(shù)演化的二級相變與鐵磁Ising模型隨溫度變化的二級相變不同,而一級相變則具有相同的有限尺度效應.需要特別說明的是,公共物品博弈通過集體分紅獲得收益,而傳統(tǒng)的囚徒困境和雪堆博弈通過兩兩相互博弈獲得收益.雖然二者的博弈形式不同,但都可以轉化為與之等價的Ising模型,囚徒困境和雪堆博弈對應的Ising模型已經在參考文獻［28］中給出.

［1］Szabó G,Fath G 2007Phys.Rep.446 97

［2］Szabó G,Borsos I 2016Phys.Rep.624 1

［3］He D R,Liu Z H,Wang B H 2009Complex Systems and Complex Networks(Beijing:Higher Education Press)p73(in Chinese)［何大韌,劉宗華,汪秉宏 2009復雜系統(tǒng)與復雜網絡(北京:高等教育出版社)第73頁］

［4］Nowak M A 2006Science314 1560

［5］Nowak M A,May R M 1992Nature359 826

［6］Nowak M A(translated by Li Z Q,Wang S C)2010Evolutionary Dynamics:Exploring the Equations of Life(Beijing:Higher Education Press)pp108–125(in Chinese)［Martin A.Nowak著 (李振清,王世暢 譯)2010進化動力學-探索生命的方程(北京：高等教育出版社)第108—125頁］

［7］Hardin G 1968Science162 1243

［8］Feeny D,Berkes F,McCay B J,Acheson J M 1990Hum.Ecol.18 1

［9］Palfrey T R,Rosenthal H 1984J.Public Econ.24 171

［10］Wang X W 2015M.D.Dissertation(Anhui:University of Science and Technology of China)(in Chinese)［王旭文2015博士學位論文(安徽:中國科學技術大學)］

［11］Wang B K 2013M.S.Dissertation(Xi’an:Xidian University)(in Chinese)［王?？?2013碩士學位論文 (西安:西安電子科技大學)］

［12］Gao J 2012Ph.D.Dissertation(Xi’an:Xidian University)(in Chinese)［高佳 2012博士學位論文 (西安:西安電子科技大學)］

［13］Hauert C,Traulsen A,Brandt H,Nowak M A,Sigmund K 2007Science316 1905

［14］Dreber A,Rand D G,Fudenberg D,Nowak M A 2008Nature452 348

［15］Helbing D,Szolnoki A,Perc M,Szabó G 2010New J.Phys.12 083005

［16］Helbing D,Szolnoki A,Perc M,Szabó G 2010Phys.Rev.E81 057104

［17］Ohtsuki H,Iwasa Y,Nowak M A 2009Nature457 79

［18］Szolnoki A,Szabó G,Perc M2011Phys.Rev.E83 036101

［19］Hauert C,Doebeli M 2004Nature428 643

［20］Santos F C,Santos M D,Pacheco J M 2008Nature454 213

［21］Szabó G,T?ke C 1998Phys.Rev.E58 69

［22］Szabó G,Hauert C 2002Phys.Rev.Lett.89 118101

［23］Szabó G,Vukov J,Szolnoki A 2005Phys.Rev.E72 047107

［24］Wang W X,Ren J,Chen G,Wang B H 2006Phys.Rev.E74 056113

［25］Gao K,Wang W X,Wang B H 2007Physica A380 528

［26］Liu Y K,Li Z,Chen X J,Wang L 2009Chin.Phys.Lett.26 377

［27］Liu Y K 2010Ph.D.Dissertation(Xi’an:Xidian University)(in Chinese)［劉永奎 2010博士學位論文 (西安:西安電子科技大學)］

［28］Yang B,Li X T,Chen W,Liu J,Chen X S 2016Commun.Theor.Phys.66 439

［29］McKelvey R D,Palfrey T R 1995Games Econ.Behav.10 6

［30］He J T 2006Behavioral Economics The Economic Behavior of China(Beijing:China Financial Economic Publishing House)p79(in Chinese)［賀京同 2006行為經濟學與中國經濟行為(北京:中國財政經濟出版社)第79頁］

［31］Zhuang Q,Di Z R,Wu J S 2014PloS One9 e105391

［32］Fisher M E 1967Rep.Prog.Phys.30 615

［33］Privman V 1990Finite Size Scaling and Numerical Simulation of Statistical Systems(Singapore:World Scienti fi c)p223

［34］Landau D P,Binder K 2014A Guide to Monte Carlo Simulations in Statistical Physics(United Kingdom:Cambridge University Press)p82

［35］Binder K,Heermann D 2010Monte Carlo Simulation in Statistical Physics An Introduction(Berlin Heidelberg:Springer-Verlag)p56

［36］Binder K,Landau D P 1984Phys.Rev.B30 1477

［37］Lee J,Kosterlitz J M 1991Phys.Rev.B43 3265

Phase transition properties for the spatial public goods game with self-questioning mechanism?

Yang Bo1)2)?Fan Min1)2)Liu Wen-Qi1)2)Chen Xiao-Song3)4)
1)(Data Science Research Center,Kunming University of Science and Technology,Kunming 650500,China)
2)(Faculty of Science,Kunming University of Science and Technology,Kunming 650500,China)
3)(Institute of Theoretical Physics,Key Laboratory of Theoretical Physics,Chinese Academy of Sciences,Beijing 100190,China)
4)(School of Physical Sciences,University of Chinese Academy of Sciences,Beijing 100049,China)

The spatial public goods game is one of the most popular models for studying the emergence and maintenance of cooperation among sel fi sh individuals.A public goods game with costly punishment and self-questioning updating mechanism is studied in this paper.The theoretical analysis and Monte Carlo simulation are involved to analyze this model.This game model can be transformed into Ising model with an external fi eld by theoretical analysis.When the costly punishment exists,the e ff ective Hamiltonian includes the nearest-,the next-nearest-and the third-nearestneighbor interactions and non-zero external fi eld.The interactions are only determined by costly punishment.The sign of the interaction is always greater than zero,so it has the properties of ferromagnetic Ising.The external fi eld is determined by the factorrof the public goods game,the fi neFon each defector within the group,and the relevant punishment costC.The Monte Carlo simulation results are consistent with the theoretical analysis results.In addition,the phase transitions and critical behaviors of the public goods game are also studied using the fi nite size scaling theory.The results show that the discontinuous phase transition has the same fi nite size e ff ects as the two-dimensional Ising model,but the continuous phase transitions is inconsistent with Ising model.

ising model, fi nite size scaling theory,Monte Carlo simulations,self-questioning update rules

27 May 2017;revised manuscript

4 July 2017)

(2017年5月27日收到;2017年7月4日收到修改稿)

10.7498/aps.66.196401

?昆明理工大學引進人才科研啟動基金項目(批準號：KKSY201607047)和國家自然科學基金(批準號:61573173)資助的課題.

?通信作者.E-mail:yangbo@kmust.edu.cn

?2017中國物理學會Chinese Physical Society

PACS:64.60.De,87.55.K–,02.50.Le,87.23.Ge

10.7498/aps.66.196401

*Scienti fi c Research Foundation for Introduced Scholars,Kunming University of Science and Technology(Grant No.KKSY201607047)and the National Natural Science Foundation of China(Grant No.61573173).

?Corresponding author.E-mail:yangbo@kmust.edu.cn