Kuramoto模型在Python和MATLAB中的简单实现

haidaowang

事情是这样的，我最近在研究团队编组及内部模式对发挥团队能力的影响，以及如何正确编组让团队能力发挥实现最大化，别问我为什么研究这个，反正稀里糊涂就研究上了。我发现在描述团队编组间及内部同步能力的时候，人们对Kuramoto模型（藏本模型）作了大量的研究,其中包括模型达到完全相位同步的充分条件、耦合强度对于同步的影响、一定条件下振子的收敛速率等。但具体实现一般都在MATLAB中，且网上代码过于复杂（我运行了一遍一堆报错），这里我使用Python和MATLAB对Kuramoto模型简单模拟。
模拟的话还是一遍举个栗子，边分析边测试效果最好，百度学术上有一篇关于Kuramoto模型的简单论文，我们就用它来实现模拟。
; W- P  C, j, V, o2 Q5 a  S
好吧，那我们开始，首先是Python实现：
/ V( U4 J3 a% j1 A
N维Kuramoto模型的数学描述如下：
7 Q6 t* G! ^" V1 O% [7 U$ A6 _5 e
# B, a# C6 W1 C
8 x* ]# e7 n& h7 ^* t. r0 q没错，是讨厌的数学公式，没事，它可以改写成这样：
" S; [8 L4 ^" N; k" @/ r
  P* n" O  s$ k6 i9 v: l
9 x! A) \! O; P好像还是有点长，那我们在改写一下：
3 h3 i' P9 J: a, {

( M1 G; P3 ]% v% d9 r& p& P看着好多了，那我就来说说式子中参数的意义，Kij为耦合矩阵，是为了便于描述不同振子间耦合程度不同的情形。最下面那个式子的r就是我们的目标，反应振子间的相关性，这个相关性就可以描述我们想要的编组内部同步能力。

2 h1 r3 u8 q/ t- G1 S+ l$ D- ~哎呦，这个式子看起来好简单，这里要补充一下知识点：同步能力可不是一下子各组该怎么同步直接确定的了，它是一个从开始到稳定的阶段，也就是说随时间变化，最终反映在各组的同步能力才会确定，那么最后图像是什么样子才算同步能力好呢？

  G# w" `0 M. {/ D+ H6 N$ {( _* O同步能力好，是指随着时间的推移，各组的同步能力r逐渐稳定，波动现象消失或固定在某一个小范围内。需要注意的是这和各组r值之间的差距没有关系，我们要的是一个平稳的状态。那怎么办找r和t的关系呢？
/ T/ i) G( S) [% V# O/ K% ]
注意看最上面那两个式子，相位（第一项，等号左边那个）上面有个点，这样他可就不简简单单是个相位了，它代表的是相位的变化值，是一个微小的微分值，好吧具体意思就是，那个式子左边展开之后是这样的：


哎呀，t出现了，其实\theta 与t有关，这里你可能有点绕，因为它们之间的关系是一一对应的，就是说每个时间的t对应了一个\theta ，我下面带入具体数值的时候你就知道了。

组间同步能力与时间t的关系出现了！

/ \  S7 w" T2 A, U% J) @0 c也就是说我先用上面的那个公式4计算出来\theta 的值，在带入到公式5，那么t-r关系就可以明确下来了，那现在我们再回过头来看看文章中已经给的例子，看看还有没有未知量。
* E: ~: w# g* h- }
栗子是这样的栗子：
5 O9 @6 ]3 x: O* b& d
假设某机构内部有 4 个编组，每个编组包括 5 个节点(其中 1 个节点为领导节点) 。另外，将上级领导作为一个独立的编组，且只包含一个节点。假设在领导机关增加4名信息传递人员。当以独立编组模式编组时，指定1名信息传递人员为指挥者，其指挥关系与其他编组一致; 当分散编组时，信息传递力量节点的关系与所在编组其他节点指挥关系一 致。其中，完全分散编组模式时，各信息传递力量节点之间无信息共享通道; 不完全分散编组时，在各信息传递人员节点之间建立一条信息共享通道。各编组模式及其拓扑结构如下图所示：
; g4 z5 F' p& e* g9 h# p. G- d/ P
6 b! B: J) z- _独立编组模式
/ [8 B" W! B% [
6 j( y$ b7 ]0 i) c
! S" n9 D2 n& I完全分散编组

* A: X! Z% q/ u' I% [% X6 o
1 a9 m9 L2 H2 d* W
. a; z; w: r0 s  Q5 \不完全分散编组
7 E  A) ?4 k4 D* Z. K
* @& U' n) E7 I; A! ]) U

参数数据


' {& ~; t! {5 [3 \& A5 R8 }3 {9 H0 o9 \

1 v; f2 D1 L- G参数确定一下有没有未知量：
# O4 d8 ?. v+ b0 K: {4 P
; S* c) L5 e' T. l$ J首先N，数据数目已知，这个有了。
: \* Z0 L# T9 P( r3 d
K值是分组内的连接强度，这个是看实际情况，由甲方提供或者自己看着给的，这里就是甲方给的编组图，i与j点的链接强度一目了然，这个有了。
/ H/ Q( P  t1 D7 ^, z5 t8 W
\omega _{i} 是振子i的固有频率，也称自然频率，甲方会给，没法自己估计，这个有了。
: A& s. D/ j6 A: }9 r% m( m" o
\theta ，\theta 怎么办，初始的\theta 会给，自己也能测的出来，但那么多\theta {j} -\theta{i} 得多少不知道啊，这里通过翻看文章，我发现其实文章是有一个特殊条件的，不然的话是需要研究耦合因子求三种约束条件解情况的，特殊条件就在这：

假设编组内节点的初始相位差为π/2，且编号最小者为0，随编号增大而增大。
% S+ S2 ~( W6 i
9 K! l# n3 E# C) ]3 |. A哦，初始相位差知道了，你还告诉了我各个初始相位，那么\theta _{j} -\theta _{i} 的值就在一个范围内的几个固定值里面啊！
4 R% @1 Z3 }. P
' D# w% l/ M: \' G+ x好的，没有未知量了，就是找K的时候麻烦点，没办法，这个决定了编组的不同，写脚本算一下吧：

#codingutf-8

5 r7 }8 I9 r: @) Z# q4 C##ScriptName:KuramotoSimulation.py
0 s* \5 B( N2 c; N7 [& L
import matplotlib.pyplot as plt

from pylab import

from sympy import
5 V' D4 o$ P1 X9 L! j
* r& Z) w/ i: A- Y) i+ efrom matplotlib.ticker import MultipleLocator, FormatStRFormatter

import math
9 ?5 v/ k: Z, Z6 ^, g  }
. n% S4 E, i, X$ t" kimport numpy as np

N = 31      #总节点数

  w2 X9 J1 Z+ v% v$ U+ L% t8 {c=[0,0,math.pi  2,math.pi,3  math.pi 2,0,0,0,math.pi  2,math.pi,3  math.pi 2,0,0,0,math.pi  2,math.pi,3  math.pi 2,0,0,0,math.pi  2,math.pi,3  math.pi 2,0,0,0,math.pi  2,math.pi,3  math.pi 2,0,0]
$ \# H/ R8 d) p( O' V( A. s' I
! c! e* V/ B& j( @( _3 L% ?w = [4,3,3,2,2,1,4,3,3,2,2,1,4,3,3,2,2,1,4,3,3,2,2,1,4,3,3,2,2,1,5]
4 X  m% r! k1 U1 \% K: q
2 _& v9 R/ h  D) P& L0 Dk1 = [

- d3 X3 X( s/ |[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],

[1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2],
7 p4 l! @% Y5 w7 y
[1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
1 K$ v% _# k; W* |
[1,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
% J3 K- A, N0 @
[1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],

[1,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
3 {5 [8 h/ m. a: n& Z7 D
[1,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],

[0,0,0,0,0,0,1,0.9,0.9,0.9,0.9,0.9,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2],

[0,0,0,0,0,0,0.9,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],

$ O) |: b4 F: Y& F: Q1 [[0,0,0,0,0,0,0.9,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],

$ \3 P! m$ p8 K+ G7 d[0,0,0,0,0,0,0.9,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],

6 j* W1 l' K9 D# v: m[0,0,0,0,0,0,0.9,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],

[0,0,0,0,0,0,0.9,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],

[0,0,0,0,0,0,0,0,0,0,0,0,1,0.8,0.8,0.8,0.8,0.8,0,0,0,0,0,0,0,0,0,0,0,0,2],
) m) g* V) T" F& Z6 H" v" y
( v( w( s9 a3 t: n. |& G2 F/ Z[0,0,0,0,0,0,0,0,0,0,0,0,0.8,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
0 c% Z( ]3 d& P; ^, T. O, ]7 {
[0,0,0,0,0,0,0,0,0,0,0,0,0.8,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
7 R2 b  F3 W! n
5 b% d; A0 K5 R[0,0,0,0,0,0,0,0,0,0,0,0,0.8,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
) ~8 Q3 m$ y% G* g! N# I6 B
; g& f3 i& n: d[0,0,0,0,0,0,0,0,0,0,0,0,0.8,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0],

9 h' r! E+ z$ N  U7 H8 o[0,0,0,0,0,0,0,0,0,0,0,0,0.8,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0],

[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0.7,0.7,0.7,0.7,0.7,0,0,0,0,0,0,2],
# w) C+ e+ b  t1 F* x" F' {8 k
# d2 {& Y9 B) t[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.7,1,0,0,0,0,0,0,0,0,0,0,0],
$ ?1 h( w  ^# H
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.7,0,1,0,0,0,0,0,0,0,0,0,0],
4 V# h+ l  g4 F/ y* P3 g! P) [
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.7,0,0,1,0,0,0,0,0,0,0,0,0],
7 Y) x; o6 w+ K* }/ X& A
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.7,0,0,0,1,0,0,0,0,0,0,0,0],

[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.7,0,0,0,0,1,0,0,0,0,0,0,0],

[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0.6,0.6,0.6,0.6,0.6,2],
$ b( s: l0 [2 }2 E. \* \
! [( j7 ?% h# X# E! w[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.6,1,0,0,0,0,0],
- X% K5 o/ g% ^! C) l% o  @; r
6 u& V  t) F+ T# k) K  s3 x; I[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.6,0,1,0,0,0,0],
  L" _3 M* c, b6 N# e
6 h! g# V8 d, X[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.6,0,0,1,0,0,0],

$ |. v& H, l3 H' `[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.6,0,0,0,1,0,0],
" b) E7 w; `6 O0 k& h( L/ P
" }0 Q0 ~; k& r* z, h9 y' X; s- q7 x8 f[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.6,0,0,0,0,1,0],
+ Q' F# z3 F* q: U
% F0 F: ]" O; N9 _! \$ @) j; `1 g' `[2,0,0,0,0,0,2,0,0,0,0,0,2,0,0,0,0,0,2,0,0,0,0,0,2,0,0,0,0,0,1]]   
7 k0 y) H# `& z
n = [i + 1 for i in range(22)]                #目标划分,24个值,1-24

" R  B# ]2 q! P" J! y$ t* ?- Q" vt = [j  for j in range(1000)]

' J6 ?2 d9 K$ l% d* B% S4 Zci = 0

$ B, @  S4 a/ V0 XC = []

C.append(c)
# _6 a. O6 R! z
for d in range(1100)

8 {8 n. l5 ~" B5 i6 a, |    for i in n

& ~7 |3 J4 Z+ @' K        for j in range(22)

            cj = c[j + 1] - c

4 ?. ^' ]9 Q$ U4 U  n5 l' R            ci += k1[j + 1]  math.sin(cj)
2 J4 ~$ p$ ^  q: D
        cii = ci + w

        h = [0.01  cii + z for z in c]
; z' H- J$ E7 k# u, z
        C.append(h)
) ~$ q( {6 S4 |
        c = h
# s. w% ]: q. l
0 g1 v3 k; G1 i0 i* p' g4 a& D  ~def r_function(u)

1 f9 U! h# e2 @. O! e' i0 H    y1 = 0

    y2 = 0

    y12 = 0

    y22 = 0

    r = []
# g2 K5 n. ~' H; S  i
    for x in range(1000)

        for ul in u

  G% ]3 `+ D' m- w  X, L/ F* {# q            y1 += math.cos(C[x][ul])
) ^; j. e2 K# ]  e1 W7 i- Q
            y2 += math.sin(C[x][ul])

        y12 = y1  2
; M( y0 \# {: H; o, E# N& U
        y22 = y2  2
% x2 N" ^) l' w5 i3 X' Z+ w% i9 h
        r.append((float(1)  N )  ((y12 + y22)  0.5))
. t6 e" u6 m0 U, E' T' A
    return r
$ Y: Z+ Q& N6 k% G
  }0 l) S. C" v+ p) R; M  ir1 = r_function([1,2,3,4,5,6])
2 e9 U! U' Y( x  X+ q5 v
r2 = r_function([7,8,9,10,11,12])

r3 = r_function([13,14,15,16,17,18])

( C: \- d# U$ e  l4 _r4 = r_function([19,20,21,22,23,24])
+ l2 Q2 `9 y( ~0 c( p5 L! `3 M
$ _, I. N8 {: u% B1 x' Wr5 = r_function([25,26,27,28,29,30])
7 l3 H6 x4 N* V7 \5 b% Z
. ?1 b2 {* E, X! d+ V+ b& R#r6 = r_function([31])
+ Y) A; o. W% e; m" v: a
. Q7 V/ ?6 Y8 |2 ?$ M

2 i7 A- U+ b6 D1 yax = subplot(111) #注意一般都在ax中设置,不再plot中设置

% r4 u1 T7 L3 r  v7 j5 G& ~ymajorLocator  = MultipleLocator(0.1) #将y轴主刻度标签设置为0.5的倍数
  P5 P+ r) I8 P2 m0 \
* ?- L4 O- H, V* t" nax.yaxis.set_major_locator(ymajorLocator)

plt.plot(t, r1, marker='o', color='green', label='1')
1 @, e. b+ y9 }
& H$ B. q5 A; W: _  A* Bplt.plot(t, r2, marker='D',color='red', label='2')
% T! P8 ~" V, }: K
plt.plot(t, r3, marker='+',color='skyblue', label='3')

( T/ q7 f3 u+ O, U/ c- nplt.plot(t, r4, marker='h', color='blue', label='4')

plt.plot(t, r5, marker='_',color='yellow', label='5')

' ~" G$ S- N2 \/ m) r5 E$ z  a#plt.plot(t, r6, color='red', label='6')

) \5 [5 R8 k  E) L* ]4 a1 A

% A* \' j5 c. W- p% L$ nplt.legend() # 显示图例
6 w9 v5 B. \1 Y6 n$ W
5 c6 G# b* I3 @" B1 g5 s4 C7 q+ u# ~plt.xlabel('iteration times')
2 [/ \+ ]+ z& _' A
plt.ylabel('r')

plt.show()

7 y; u- P3 ?' y
独立编组结果如图：

, L7 `5 J) D! p$ \
: ]1 H; P6 E9 `- S, g) P8 L" p
k1独立编组
' s% k1 E# q& z9 s3 C5 O好的，从图像我们来看看Kuramoto模型在描述这个编组的时候，5组最终稳定，我们说这个团队编组还算科学，但我们改变一下K的值，换成分散编组：
4 g! H) T/ r) {' L# O" Q& z2 t
3 j% `# i8 }( _6 D" w. e3 D* c( Yk2 = [

  a, a7 @( d( l5 E' N! J[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],

% a5 k" X5 O2 S# i[1,1,0.9,0.8,0.7,0.6,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2],
: m& q$ _5 b1 i8 K7 V5 o' R- T# _
  I( Q, Q! K9 D' ~7 q& l[1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],

[0.9,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
. A+ C" L2 S% j- _" s6 x' Z
3 O& b' o+ t2 k8 l: ^. N7 J# z[0.8,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
  T' [) F4 W  K- U& @3 D
[0.7,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
! y/ w8 F2 Q) U. J
0 W8 X9 r4 d4 _0 [[0.6,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],

[0,0,0,0,0,0,1,1,0.9,0.8,0.7,0.6,0,0,0,0,0,0,0,0,0,0,0,2],
; Y7 a! ?/ `( J6 y+ W
[0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],

[0,0,0,0,0,0,0.9,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
5 n) M6 ~) T0 A" i
[0,0,0,0,0,0,0.8,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],

[0,0,0,0,0,0,0.7,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
: a) U; K4 D: v
[0,0,0,0,0,0,0.6,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0],

[0,0,0,0,0,0,0,0,0,0,0,0,1,0.9,0.8,0.7,0.6,0,0,0,0,0,0,2],
/ n4 u/ `1 a  h& m9 ]+ x9 U) j) }
[0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0],
2 ?8 v. ?3 G: P8 ?, t/ @, p
[0,0,0,0,0,0,0,0,0,0,0,0,0.9,0,1,0,0,0,0,0,0,0,0,0],
0 T6 Y; r: D/ f7 F* `% \% e
[0,0,0,0,0,0,0,0,0,0,0,0,0.8,0,0,1,0,0,0,0,0,0,0,0,0],
% }+ `* [2 C! M6 m
( l7 U. R* `. A1 Z/ J+ w# _[0,0,0,0,0,0,0,0,0,0,0,0,0.7,0,0,0,1,0,0,0,0,0,0,0,0],

[0,0,0,0,0,0,0,0,0,0,0,0,0.6,0,0,0,0,1,0,0,0,0,0,0,0],

[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0.9,0.8,0.7,0.6,2],

1 X$ A5 _1 a. n& o[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,0,0,0,0],

[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.9,0,1,0,0,0,0],

# x& x# D' T2 _1 J0 @6 F7 O[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.8,0,0,1,0,0,0],

[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.7,0,0,0,1,0,0],

[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.6,0,0,0,0,1,0],

1 U- p  ]7 O6 |( B) I[2,0,0,0,0,0,2,0,0,0,0,0,2,0,0,0,0,0,2,0,0,0,0,0,1]
( J& D; m+ u5 ?8 ^, H; S


k2独立编组
4 u( K' T+ b, e* x4 U4 f5 {

1 F6 e2 |) L6 k" S
. Z  l! ]* O/ N# v1 U2 s这两幅图像都在开始阶段大幅波动，而后在一定范围内趋于稳定，那么到底哪个分组模式最符合实际，最能突出编组能力呢？

% `( U2 _' u% G+ a, ~/ e) p这里还有一个公式，来解决这个问题，编组同步能力的量化：
6 {" s5 E2 M; [7 I7 G$ a
. h# q/ d: B5 g5 E
M_{s} 就可以描述某编组的同步效果，r_{s} 是达到稳定状态后序参量的均值，β∈(0,1)是调节因子。我们可以用M_{s} 比较编组内部的好坏。那编组间能力的好坏怎样比较呢？
, M/ t. i6 j; y* X& u; w6 f
5 \' A$ S8 e& k- u, B8 ~; r这个Kuramoto模型同样有所考虑，它有一个描述整个系统编组能力的公式：

- _) d, J/ u; }9 o; P
/ Z0 T6 u! X# g其中，P是编组的数量，M_{i} 是第i个编组的同步能力，\omega _{i} 是编组在整个系统中的权重，\psi {e} 是各编组平均 相位的均值，\Delta \psi 是各编组平均相位的标准差。具体的计算不是这篇文章的重点，就不在计算M{s} 和M的值来比较上述例子独立编组和分散编组的好坏了，本篇文章主要是讲下Kuramoto模型的解决思路，尤其是上面解决\theta 值的方法可以套用在其他Kuramoto模型中，做一个目标估计绰绰有余的。

6 Z9 N# k# e. A' R" m# r下面是解决Kuramoto模型常用的MATLAB编程方法，具体思路与上述基本一致，这里不再赘述，K的值我们给另一种编组模式：不完全分散编组模式，也是现在实际上最长用的编组方式，直接上代码：

  e) p! [& B* j( Dclc;
* U! u" r+ p4 X7 U* K' K* a3 o
+ C: K4 P3 Y. ?) K+ {+ Uclear all;

% Y# q  I' ~, d! ~; lk=[0  1  1  1  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0.5;
, Y, X6 M; X( @3 g% `) ]% ^
  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0;   
$ V3 u5 {- S1 K5 c$ R& l$ B/ j
6 e$ X1 K0 J% c/ c% Z  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0;   
6 R4 P! T# I! w
& g8 Q* c- R* @% o! b5 v& G- ~6 O- e  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0;   
5 V* \1 F# K2 [) |
  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0;   
1 ]0 z* ~- Q$ C6 Z& ^0 L0 A& ?
5 c$ }7 e4 J8 j8 I7 f! F# ^  0  0  0  0  0  0  1  1  1  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0.5;
2 K2 E2 Y- r7 v
9 v9 i' m; v# Q# m  0  0  0  0  0  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0;   
, `7 A2 V" |5 ~; l5 U; G
  0  0  0  0  0  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0;   

- }& P4 w" t( ^  0  0  0  0  0  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0;   
7 r3 s) ^$ z& ^9 J
) I* S/ j. p" c; _' P) B! s  0  0  0  0  0  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0;   

  0  0  0  0  0  0  0  0  0  0  0  1  1  1  1  0  0  0  0  0  0  0  0  0  0.5;

3 n# q/ {; v* c; Y  0  0  0  0  0  0  0  0  0  0  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0;   
& H/ J# I" e# S! _+ {/ d( [
  0  0  0  0  0  0  0  0  0  0  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0;   
% q$ o5 M# G6 `
  0  0  0  0  0  0  0  0  0  0  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0;   
" e& U" c/ w. z7 P! k+ F/ L% l
; {' t, M$ z8 \  0  0  0  0  0  0  0  0  0  0  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0;   
& P' q4 j7 `" ?' `% e
  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  1  1  1  1  0  0  0  0  0.5;
: q2 d( y6 t! B
  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  1  0  0  0  0  0  0  0  0  0;   

( Z- r7 A7 X4 g  q$ k" ~" r9 Z  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  1  0  0  0  0  0  0  0  0  0;   
1 O, v% ^% }6 |& U0 M, w
  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  1  0  0  0  0  0  0  0  0  0;   

  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  1  0  0  0  0  0  0  0  0  0;   
8 _: ]  F3 m0 f0 d: N
  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  2  2  2  0.5;
4 m0 K4 m/ F& y$ l
( `' Z. T; _" [  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  2  0  0  0  0;   
, I5 h! N' ]3 R% M! J
  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  2  0  0  0  0;   
9 u) V/ E, t1 A8 R( }: w' `2 [6 W. i
  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  2  0  0  0  0;   
; t) r: m- [: N2 F+ O& v9 A6 K
2 `( e! k/ h7 x  0.5 0 0  0  0  0.5 0 0  0  0 0.5 0  0  0  0  0.5 0 0  0  0  0.5 0 0  0  0;   
" }) }- I% W& X
, c( r5 ?2 K. I# q& x  ]     
& q1 W& i. z/ b4 r3 H
$ e, x: @6 l. _, T0 p- o+ GN = 25;
6 T: w  m8 Z/ W! |, \
4 u. }1 u) O  s+ Z! T$ ^! eStep= 10000;
: w" W4 j" P/ e0 ^& V
Theta=zeros(Step,N);
, j5 o+ p" f  \  [- t. l% A
Omega =zeros(Step,N);
# \3 z$ S# H: X+ o# [2 Z
% U8 R; a+ X0 `% [& KDeltaT=0.01

& s& S  C9 P. \2 F8 S( m. eTheta(1,：）=[0 pi/2 pi 3*pi/2 0 pi/2 pi 3*pi/2 0 pi/2 pi 3*pi/2 0 pi/2 pi 3*pi/2 0 pi/2 pi 3*pi/2 0 pi/2 pi 3*pi/2 0];
# g6 ~$ E  ]6 F5 Q" B+ @
4 @6 z4 d1 u$ ~, y) r9 VOmega(1,：D=2*pi*[2  3  3  4      4 2    3  3      4 4    2  3      3 4    4  2      3 3    4  4      2 2    3  4      1];

%  求解微分方程
! O; m$ e' v! H- m
6 v+ ^' B5 c% xfor n = 1:Step
* P6 K; O1 V1 c2 q
    for i = 1:N

1 v/ a7 x; m6 l1 y8 M+ n* M            Sigma = 0;
5 Y8 {+ `3 \  n6 P6 P- S, o
        for j = 1:N
$ `. P  [) u/ ]! l
. m. S$ M' I5 j* K4 j! d1 J0 }            %%判断k对同步效果的影响
4 P& K2 G* Y, _) I4 h# V% Y
= Sigma + k(i,j) * sin( Theta(n,i) - Theta(n,j) );
3 r* w5 h7 u2 O' w, f
' \& J4 m" l; M' x7 M        end

            Omega(n+1,i) = Omega(n,i) - Sigma;

            Theta(n+1,i) = Omega(n+1,i)*DeltaT + Theta(n,i);
" ^8 Y* P4 w" i: T8 u; D
. `; Y3 {3 h* @    end

end
$ K% d% F9 J7 U2 _& K* V6 z
! l. [: ^% o0 w) r%  求解序参量r(t)
- D$ b! |5 R& q. R8 A
groupN = 5;
  ~1 d5 {8 ~, o' c! V# S
Step = 1000;

for i = 1:Step

- I& j2 @" j( r& l    sigma1 = 0;

    sigma2 = 0;

    sumTheta = 0;
6 k2 e) C+ A7 l6 y2 `+ q
    for j =6:10  %表示相应的群组
0 I! k) s. b- M4 n: H1 t# P2 D
- z. J/ S$ v7 S% r+ v: W        sigma1=sigma1+cos(Theta(i,j));
/ K5 @- u; i9 a
        sigma2=sigma2+sin(Theta(i,j));

    end

    r(i)= sqrt( sigma1^2+sigma2^2 )/groupN;

) {3 T4 F6 N1 t% X. \0 _8 lend
! T9 A+ @2 T: c+ k  K6 T
x= 0.01：DeltaT：DeltaT*Step;
* ?  Y' \6 [: D2 D
plot(x,r);
MATLAB版不完全分散编组结果如下图：
3 M/ [# v; i2 X' `

regngfpcb

Kuramoto模型在Python和MATLAB中的简单实现