模擬工具及方法

1 1 1 (1 ) 1

b I b b I

S c a S

I c I I I c I

    

    

，

a33= _{2 2 2 2} ^{2 2}

2 1 1

( )( 1 )

1 1

I a b I

I c I

 

   

 

，

a =34 _{2 2 2} ^{2 2 2} _{2 1 2 1} ^{1 1 2}

2 2 1 1 1 2 2

1 1

[( )( )] [ ( )( ) ]

1 (1 ) 1 1 1

b b I b I

S S S c a

I I c I I c I



     

    

，

a41=0，

a42=0，

a43=



_{2 2}I ， a44=



₂S₂

 

₂。

所求得的平衡點代入如上 Jacobian matrix，可以得到一 4X4 的矩陣，接著 可以求得特徵值，從特徵值分析平衡解在該情形的局部穩定性。

第三節 模擬工具及方法

我們使用數學軟體 Maxima 計算出平衡解，再用 Matlab 程式語言執行數值模 擬，分析模型(3.2)是否因為係數 , ,a b c 不同而影響動態情形。為便於了解確切_{i i i} 影響模型的係數，須先選擇一項係數變動，其餘係數固定不變。根據動態情形來 分析該係數對模型(3.2)的影響。再固定模型(3.2)中的部分係數，只變動c ，根₂ 據觀察平衡解的變化，判斷c 對模型(3.2)的影響。 ₂

24

例子 1 ：

比較c (_i

i  1, 2

)值不相同時，兩棲息地上可能被感染生物最終情形，固定 5

.

1



0

A ，



₁



0.2，



₁



0.7，



₁



0.5，A₂



0.4，



₂



0.2，



₂



0.8，



₂



0.5，

1



6

a ，b₁



0，c₁



0，a₂



6，b₂



0：

(1)當c₂



0時，可以得到穩定的平衡解E₂



(0.685, 0, 0.625, 1.2758)，此種平 衡解屬於 DF+E(disease free+ epidemic)。其動態演化圖如【圖 4.1】所 示：

0 50 100 150 200 250 300 350 400 450 500

0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6

time t

S1 I1 S2 I2

【圖 4.1】E₂



(0.685, 0, 0.625, 1.2758)的模擬。

另外，檢驗是否有其他平衡解存在，如圖【圖 4.2】，經由計算得到有邊界 平衡E₁



(0.714, 1.212, 0.7557, 0)，

【圖 4.2】c₂



0時，模型另外一平衡解E 位置。 ₂

25

經由計算得到兩平衡解的特徵值(eigenvalue)，做成如下【表 4.1】:

【表 4.1】c₂



0時，兩平衡解的特徵值及穩定性。

2 0

c



E ₁ E ₂

座標

(0.714,1.212, 0.7557, 0) (0.685, 0, 0.625,1.2758)

特徵值

1 12.6212

  

2 0.3137 0.3317i

   

3 0.3137 0.3317i

   

4 0.1046

 

1 12.7102

  

2 0.3552 0.3504i

   

3 0.3552 0.3504i

   

4 0.0205

  

穩定性 不穩定 穩定

2 0

c



時，兩平衡解中，E 的特徵值皆為負數，故₂ E 為穩定平衡解；而₂ E 的特₁ 徵值



₄



0.1046為正數，故E 為不穩定平衡解。 ₁

將E₂



(0.685, 0, 0.625, 1.2758)，E₁



(0.714, 1.212, 0.7557, 0)兩平衡解標 記在三維圖上，並模擬平衡解周圍的動態情形，得到如下圖型【圖 4.3】。

0.6 0.65 0.7 0.75 0.8 0

0.5 1

1.5

0 0.5 1 1.5 2

I1 S1

E1 E2

I2

【圖 4.3】c₂



0時，兩平衡解及其周圍動態情形。

在三維圖型中，由上圖【圖 4.3】可看出在E ，₁ E 兩平衡解周圍的起始點₂ 最終都會趨向E 。由此可知₂ c₂



0時，平衡解E 是模型(3.2)中的穩定平衡點。 ₂

26

(2)當c₂



1時，可以得到穩定平衡解E^*



(0.714,1.065,0.625,0.2)，此種平衡解屬 於 E+E(epidemic+epidemic)。其動態演化圖如【圖 4.4】所示：

0 50 100 150 200 250 300 350 400 450 500

0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6

time t

S1 I1 S2 I2

【圖 4.4】E^*



(0.714,1.065,0.625,0.2)的模擬。

另外，檢驗是否有其他平衡解存在，如圖【圖 4.5】，經由計算得兩個邊界 平衡E₁



(0.713, 1.212, 0.755, 0)和E₂



(1.337, 0, 0.625, 1.015)。

【圖 4.5】 c₂



1時，模型另外兩平衡解E 和₁ E 位置。 ₂

27

經由計算得到兩平衡解的特徵值(eigenvalue)，做成如下【表 4.2】:

【表 4.2】c₂



1時，兩平衡解的特徵值及穩定性。

2 1

c



E ^* E 1 E ₂

座標

(0.714,1.065, 0.625, 0.2) (0.713,1.212, 0.755, 0) (1.337, 0, 0.625,1.015)

特徵值



₁

 

11.5787

2 0.3268+0.3549i

  

3 0.3268 0.3549i

   

4 0.0767

  

1 12.6212

  

2 0.314 0.3316i

   

3 0.314 0.3316i

   

4 0.104

 

1 9.5333

  

2

= 0.3282 0.2476i

  

3

= 0.3282 0.2476i

  

28

(3)當c₂



100時，可以得到穩定平衡E^*



(0.714, 1.262, 0.625, 0.0023)，此種平 衡解屬於 E+E(epidemic+ epidemic)。其動態演化圖如【圖 4.7】所示：

0 50 100 150 200 250 300 350 400 450 500

0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6

time t

S1 I1 S2 I2

【圖 4.7】E^*



(0.714, 1.262, 0.625, 0.0023)的模擬。

另外，檢驗是否有其他平衡解存在，如圖【圖 4.8】，經由計算得到兩個邊 界平衡E₁



(0.714, 1.212, 0.755, 0)和E₂



(3.739, 0, 0.625, 0.054)。

【圖 4.8】 c₂



100時，模型另外兩平衡解E 和₁ E 位置。 ₂

29

經由計算得到兩平衡解的特徵值(eigenvalue)，做成如下【表 4.3】:

【表 4.3】c₂



100時，兩平衡解的特徵值及穩定性。

2 100

c



E ^* E 1 E ₂

座標

(0.714,1.262, 0.625, 0.0023) (0.714,1.212, 0.755, 0) (3.739, 0, 0.625, 0.054)

特徵值



₁

 

11.43

2 0.0871

  

3 0.323 0.356i

   

4 0.323 0.356i

   

1 12.6212

  

2 0.3137 + 0.3317i

  

3 0.3137 0.3317i

   

30

(4) 當c₂



200時，可以得到穩定平衡E^*(0.714,1.263,0.625,0.00116)，此種平衡 解屬於 E+E(epidemic+ epidemic)。其動態演化圖如【圖 4.10】所示：

0 50 100 150 200 250 300 350 400 450 500

0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6

time t

S1 I1 S2 I2

【圖 4.10】E^*(0.714,1.263,0.625,0.00116)的模擬。

另外，檢驗是否有其他平衡解存在，如圖【圖 4.11】，經由計算得到兩個 邊界平衡E₁(0.714,1.2119,0.755,0)和E₂(0.714,0,0.755,0.0277)。

【圖 4.11】 c₂



200時，模型另外兩平衡解E 和₁ E 位置。 ₂

31

經由計算得到兩平衡解的特徵值(eigenvalue)，做成如下【表 4.4】:

【表 4.4】c₂



200時，兩平衡解的特徵值及穩定性。

2 200

c



E ^* E 1 E ₂

座標

(0.714,1.263,0.625,0.00116) (0.714,1.2119,0.755,0) (0.714,0,0.755,0.0277)

特徵值



₁

 

11.4215

2 0.323 0.3567i

   

3 0.323 0.356i

   

3 0.3137 0.3317i

   

32

從本節中，我們可以指出c 的變化可以影響系統最終的平衡解變化。在其它₂ 因素固定下，當c₂



0時，會有 DF+E 的平衡解；當c₂



1時，會有 E+E 的平衡解；

當c₂



100時，會有 DF+E 的平衡解。 下圖【圖 4.13】為c 影響模型(3.2)的過 ₂ 程。 v

【圖 4.13】c 影響模型(3.2)的過程。 ₂

第四節 結論

從上章節中，可以得知模型(3.2)會因c 的改變而影響模型動態，使得模型₂ 的穩定平衡解變動。當c₂



0時，穩定平衡解為E ，在周圍的起始值被吸到₂ E 。₂ 當c₂



1時，穩定平衡解為E ，在周圍的起始值被吸到^* E 。當^* c₂



100時，穩定 平衡解亦為E ，在周圍的起始值被吸到^* E ，但^* E^*



0.714,1.262,0.625,0.0023



十分 接近邊界值



0.714,1.262,0.625,0



。當c₂



200時，穩定平衡解亦為E ，在周圍的^* 起 始 值 被 吸 到 E ， 但 此 時^* E^*



0.714,1.263,0.625,0.00116



更 加 接 近 邊 界 值



0.714,1.263,0.625,0



。由此可知，當c 越大，₂ I 會越接近 0。₂ c 是為第二棲息地₂ 的媒體影響效力，而I 為第二棲息地上確定受感染的生物數量，當₂ c 越大，第二₂ 棲息地上受感染生物數量就會越少。

也可以固定其他係數，只變動a 、_i b 或其他係數，來分析該係數對於模型的_i 穩定性影響。或是修改模型，來對用模型對流行疾病做更進一步的了解。

33

在文檔中 非常速遷徙SIR疾病傳染模型的穩定性分析 (頁 33-43)