最優(yōu)捕魚策略

1、最優(yōu)捕魚策略模型摘要本文涉及的問題是漁業(yè)資源可持續(xù)發(fā)展，即在我國(guó)一定漁場(chǎng)，在一段時(shí)間內(nèi)，如何實(shí)現(xiàn)最大收獲量的問題；同時(shí)保證漁場(chǎng)能實(shí)現(xiàn)穩(wěn)定生產(chǎn)。我們的解題思路就是：以漁場(chǎng)生產(chǎn)過程中的兩個(gè)相互制約的因素，分別是年捕撈量與再生產(chǎn)能力，從而確定最優(yōu)發(fā)展策略 : 用微分方程分別描述各齡魚1 群數(shù)量雖時(shí)間變化的規(guī)律，并在此基礎(chǔ)上確定總效益即總收獲量為目標(biāo)函數(shù)，以漁場(chǎng)的可持續(xù)捕獲為約束條件，分別對(duì)長(zhǎng)期生產(chǎn)和固定期生產(chǎn)建立規(guī)劃模型。針對(duì)問題一：通過對(duì) 4 齡魚在年末的不同狀態(tài)（ 全部死亡；仍為 4 齡魚 ） 的考慮，從可可持續(xù)捕獲條件出發(fā)，分別建立2 個(gè)模型。最后求解在計(jì)算機(jī)上實(shí)現(xiàn)。針對(duì)問題二：確定一個(gè)整體效

2、益函數(shù)，綜合考慮年捕撈能力和年再生產(chǎn)能力，用計(jì)算機(jī)數(shù)值解法進(jìn)行搜索逐年確定各年的最優(yōu)策略，從而得出五年的總最優(yōu)策略。1. 先假設(shè)每年捕撈量強(qiáng)度相等，建立一個(gè)簡(jiǎn)單模型；2. 再假設(shè)每年捕撈強(qiáng)度不相等，建立一個(gè)復(fù)雜模型；3. 最后給出魚群生產(chǎn)能力破還不大的含義（即魚群減少率的上限）在它的約束之下再建立一個(gè)模型。關(guān)鍵詞 ：微分方程；捕撈強(qiáng)度；再生產(chǎn)能力；規(guī)劃模型1、 問題的重述為了保護(hù)人類賴以生存的自然環(huán)境，可再生資源（如漁業(yè)、林業(yè)資源）的 開發(fā)必須適度以可持續(xù)發(fā)展。那么一種合理、簡(jiǎn)化的策略是在可持續(xù)收獲的前 提下，追求最大產(chǎn)量或最大效益。要求研究的問題是：在一段期間中，對(duì)某種 魚的最優(yōu)捕撈策略。1

3、.1 魚的情況假設(shè)這種魚分為4個(gè)年齡組：1, 2, 3, 4齡魚。各年齡組每條魚的平均重 量分別為5.07,11.55,17.86,22.99（克）；各年齡組的魚自然死亡率均為 0.8（1/年）；這種魚為季節(jié)性集中產(chǎn)卵繁殖，平均每條 4齡魚的產(chǎn)量為L(zhǎng)109 乂 107個(gè)），3齡魚的產(chǎn)卵量為這個(gè)數(shù)的一半，2齡魚和1齡魚不產(chǎn)卵， 產(chǎn)卵和孵化期為每年的最后4個(gè)月；卵孵化并成活為1齡魚，成活率（1齡魚 條數(shù)與產(chǎn)卵總量n之比）為L(zhǎng)22 乂 10仙/ （1.22 x 1011 +n）0具體數(shù)據(jù)如下表:i加（9）r（1/ 年）增（個(gè)/條）15.070.80211.550.80317.860.80.5545

4、x 105422.990.80.5545 x 105其中，i表示i齡魚，明表示i齡魚的重量，r表示i齡魚的自然死亡率，叫 表示平均每條i齡魚的產(chǎn)卵量。又漁業(yè)管理部門規(guī)定，每年只允許在產(chǎn)卵孵化期前的8個(gè)月內(nèi)進(jìn)行捕撈作業(yè)。如果每年投入的捕撈能力（如漁船數(shù)、下網(wǎng)次數(shù)等）固定不變，這時(shí)單位 時(shí)間捕撈量將于i成正比，比例系數(shù)稱捕撈強(qiáng)度系數(shù) 乜，通常使用13mnW眼的 拉網(wǎng)，這種網(wǎng)只能捕撈3、4齡魚，其兩個(gè)捕撈強(qiáng)度系數(shù)之比為：%也=。鉉1 ki = = 0漁業(yè)上稱這種方式為固定努力量捕撈。1.2 問題1）建立數(shù)學(xué)模型分析如何實(shí)現(xiàn)可持續(xù)捕撈（即每年開始捕撈時(shí)漁場(chǎng)中 各年齡組魚條數(shù)不變），且在此前提下得到最高

5、的年收獲量（捕撈總量）2）某漁業(yè)公司承保這種魚的捕撈業(yè)務(wù) 5年，合同要求5年后魚群的生 產(chǎn)能力不能受太大的破壞。已知承保時(shí)各年齡組魚群的數(shù)量分別為，如果122,29.7,10.1,3.29 （*109條）,如果仍用固定努力量的捕撈方式， 該公司 采取怎樣的策略才能使總收獲量最高。2、 模型假設(shè)1 .假設(shè)只考慮某種魚群的繁殖和捕撈，期間不考慮魚群的遷入和遷出，即假設(shè) 大規(guī)模魚群數(shù)量增長(zhǎng)隨時(shí)間是連續(xù)變化的。2 .根據(jù)模型已知條件，假設(shè)魚群在一年的任何時(shí)間都會(huì)發(fā)生自然死亡，在一年的后四個(gè)月都會(huì)發(fā)生產(chǎn)卵。3 . 4齡魚的數(shù)量在第4年仍存活的數(shù)量占總的數(shù)量比例很小， 這里可以假設(shè)其 全部死亡，令其退出完

6、畢。4 .不考慮環(huán)境的影響，各年齡組的自然死亡率為0.8 （1/年）。符號(hào)說明a年平均固定夕匕亡率，單位：1/年%Q)：t時(shí)刻j齡魚的捕撈總重量，j=3 或 4t時(shí)間，t 0,1 H年總收獲量，即捕撈總重量jj 齡魚，j=1,2,3,4i年數(shù) i C 0,5j齡魚每條魚的平均重量，單 位：克T時(shí)刻j齡魚第i年的捕撈重量，j=3,4k年平均捕撈率，單位：1/年第i年的產(chǎn)卵量n母年產(chǎn)卵量為?。﹋齡魚在第i年時(shí)刻t的數(shù)量勺j齡魚在時(shí)刻t的數(shù)量第i年總收獲量，即捕撈總重量功0)j齡魚在年初的數(shù)量卜式八j齡魚在年末的數(shù)量問題的分析四、（一） 對(duì)于問題一的分析1 .對(duì)死亡率a的理解：我們定義平均死亡率a是

7、單位時(shí)間魚群死亡數(shù)量與現(xiàn)有魚群數(shù)量的正比例系數(shù)，由假設(shè)可知，它是一個(gè)環(huán)境等其他因素?zé)o關(guān)的常數(shù)。由于魚群的數(shù)量是連續(xù)變化的，而1、2齡魚全年以上及3、4齡魚在后4個(gè)月的數(shù)量只與死亡率有關(guān)，與其 他因素?zé)o關(guān)，設(shè)魚群量為 x,則在時(shí)間t,t+ ?t內(nèi)，魚群數(shù)量的減少二魚群死亡數(shù) 量，即L * 八 八 八工(t +?t) 一工雙t + ?t) -x(t)=-必t)?t, :=-ax(t).* L當(dāng)“一 0時(shí)得：=-ax(t),(1.1)2 .對(duì)于捕撈強(qiáng)度系數(shù)k的理解：題目告訴我們，捕撈強(qiáng)度系數(shù) K一定，且只在捕撈期內(nèi)(即每年的前8個(gè)月)捕撈3、4齡魚，因此只會(huì)影響3、4齡魚群數(shù)量，而不會(huì)影響其他的魚群

8、數(shù)量。我們可 以看3、4齡魚魚群的數(shù)量在捕撈期內(nèi)不僅與 k有關(guān)，而且還與死亡率a有關(guān)，類似第一點(diǎn)的分析，可以得到(1.2)(1.3)-；=。4小：；7=1 12其中x(t)表示3、4齡魚的魚群數(shù)量，t表示每年的前8個(gè)月，即tC 0 3 .對(duì)于持續(xù)捕獲的理解：隨著時(shí)間的推移，各年齡組的魚數(shù)量必將發(fā)生變化，但持續(xù)捕獲要求每年開始捕撈時(shí)漁場(chǎng)中各年齡組魚群條數(shù)不變， 再根據(jù)魚群的生長(zhǎng)規(guī)律，我們可以得到關(guān)系式:上一年齡組魚群年底的數(shù)量=下一年齡組魚群年初的數(shù)量(1齡魚除外)，即(1.4)X+ 1)4 .對(duì)成活率m的應(yīng)用又假設(shè)知，此種魚在每年的8月底一次產(chǎn)卵完畢，又已知3、4齡魚每條產(chǎn)卵的個(gè)數(shù), 因此可將

9、每年的產(chǎn)卵量n表示為：n = 1*109 x 10 x + 勾(I)(1.5)又已知成活率(1.6)122 x IO111.22 x 10n產(chǎn)卵量乂成活率=1齡魚每年年初的數(shù)量，即n x m =(1.7)5 .對(duì)最高收獲量的描述根據(jù)第2點(diǎn)的分析，在t時(shí)刻的捕撈重量s(t) =3齡魚捕撈重量#3)+4齡魚捕撈重量回(力.而53(E)=0.42k43)g3,54(t)=k#4(t)，則5(。= 0.42必3)外+ 心 4(。94(1.8)由于捕撈被看成連續(xù)的作業(yè)，因此捕撈總收獲量即年收獲量可以用t時(shí)刻的捕撈重 量s(t)關(guān)于t在捕撈期內(nèi)的積分來刻畫2 H=/s(t)dt,(1.9)要求最高的年收獲

10、量即求H的最大值。6 .對(duì)4齡魚在年末情形的兩個(gè)假設(shè)：(1)認(rèn)為4齡魚在年末與魚群總數(shù)量相比十分微小，它們既不產(chǎn)卵又不會(huì)被捕撈，可 以將它們忽略不計(jì)，令其退出系統(tǒng)。(2)近似的認(rèn)為在年末未死亡的 4齡魚的各個(gè)特征(例如：重量、產(chǎn)卵個(gè)數(shù)等)均不 發(fā)生改變，即仍回到4齡魚組中。7 .模型建立大綱(1)以第6點(diǎn)的第一個(gè)假設(shè)為基礎(chǔ)，建立一個(gè)簡(jiǎn)單的模型I ,其實(shí)實(shí)質(zhì)上是聯(lián)立以上分析的幾個(gè)方程成為一個(gè)方程組，其中；U ：一 ：.(2)以第6點(diǎn)的第二個(gè)假設(shè)為基礎(chǔ)，改進(jìn)模型 II ,使得方程組中的一個(gè)方程xo(4)二右被方程Xo=前(3) + xi(4)代替即可。(二)對(duì)于問題二的分析1 .問題一的相似之處由

11、于對(duì)各年齡組魚群數(shù)量起到影響作用的各因素(例如：平均死亡率、成活率、捕撈期等)不變，因此，在每年各年齡組的魚群數(shù)量變化情況與問題一類似。2 .與問題一的不同之處，主要區(qū)別在于(1)問題一要求持續(xù)捕獲，問題二只要求魚量不受到太大的破壞，不限制各年齡組年初魚群的數(shù)量，因此作為約束條件的方程組中各年齡組的魚量肯定與年數(shù)有關(guān)，而不像 問題一是常量。(2)問題一中的各變量呈周期變化，因此，只要考慮一個(gè)周期的變化情況即可，而問題二則不同，其各年的初始值在變化，因此，要考慮每一年的捕獲量，再講5年求和，得到一個(gè)目標(biāo)函數(shù)。綜合以上兩點(diǎn)，可以得到一個(gè)優(yōu)化問題。3.根據(jù)優(yōu)化問題我們又提出了三個(gè)模型模型I :簡(jiǎn)化使

12、得每一年的捕撈強(qiáng)度系數(shù)相同，化為一元函數(shù)最優(yōu)值的求解問題。模型II :考慮每一年的捕撈強(qiáng)度系數(shù)不同，得到一個(gè)多元函數(shù)最優(yōu)值的求解問題。五、模型的建立問題一模型I :假設(shè)4齡魚在年底推出系統(tǒng)和連續(xù)捕獲前提| tE|OJ下如何得到高年收獲量。由問題分析可知，可得以下優(yōu)化問題:maxH=Js(t)dt，s(t)見問題分析,(2.1)山白draft)丁2(。J- =- ax3(t) - 0.42 攵acdX3(t)dt =。刈dX4(t)t =_ ax4(t) - kx4 atdx4(t)dt =-叫t e 0,1r 21t E 0,-2，于 1R fl211.22 x io11-2 x 0A2k-

13、2-2x D. 21.109X100.53 +C 3(2.2)求解過程見附錄。將（2.2）式代入目標(biāo)函婁X再利用一維搜索法求元函數(shù)的最小值的方法上機(jī)求得H的最大值，k=17.4664, H=3.886e+011 （噸），即k=17.4664時(shí)，H取最大。模型II :進(jìn)一步假設(shè)4齡魚在年末的特征不變，仍當(dāng)做4齡魚，則在持續(xù)捕撈的情況下, 求得最大捕獲量。此模型類似模型I也可得優(yōu)化問題，區(qū)別僅在式中，應(yīng)改為Xo(4) =xi(3) + xi (4)同理得1.22 x IO111.22 x IO11X3 - 2 - 2 X 0.42F723鏟- -(0.42k - a) e1.109 x IO5 0

14、.5e 3+2k一 k1-e 3 J重復(fù)模型II的步驟解得k=17.4664, H=3.886e+011 (噸)。 問題二：根據(jù)題意，我們既要保證五年后魚群的生產(chǎn)能力不能受到太大破壞，又要使公司總 收獲量最高。但是如果同時(shí)滿足這兩個(gè)條件，總的來說很難把握兩者之間的度。因此， 我們采用首先使捕魚收獲量最大，其次再看其破壞程度是否影響魚的生長(zhǎng)，如果破壞程 度能保證不影響魚的生長(zhǎng)，則使用該方法。即五年合同到期后魚群盡可能接近可持續(xù)魚 的情況下使捕撈量達(dá)到最大。根據(jù)問題的分析可以得到以下優(yōu)化方程:maxH=maxS；= 1（?國(guó) + “國(guó)）4式。I 出一式。dtd 看 3。）dt=一%1(。=-叫,式

15、。=- a均,3) - 042吊修,3d43dt=dXt.4（0at- =一。丸4（。一 年笛,4 心”。21 t e o.q te 1 ，e 閱 唱1我們將此問題分成了兩種情況進(jìn)行處理情況一：假定五年中捕撈強(qiáng)度k 一直保持不變，即i=k,且第6年初的各年齡魚數(shù) 與平衡捕撈時(shí)的年初值發(fā)。0）0=1,2,3,4）的偏離不超過的的e倍。17.86* 042Max f(k)=t= 1 0.42k + 0.8s.t!X2(0 = XiG- 1)9一。.式f) = .2(jTk另2kX4(i) = X4(i - l)e + 3 + Xj(i - l)e1.5* 1,109* 10s-口叫 _竺Na =白

16、。一1)2 ! 口另+3 *1,109 * 10s T一標(biāo)一4。一 1月一、 1.22 * IORM + M。 的 三22*1臚+ (；+ 心 |箝-劭(。| S Xfl(0用Matlab進(jìn)行求解(代碼見附錄),結(jié)果可得：K=17.1446M=0.9989 N2=0.9987 3=0.9997 =4 =0.9500此時(shí)，最高總收獲量為1.6011K噸情況二：假設(shè)每年的捕撈強(qiáng)度系數(shù)都不相同，即 明-與i有關(guān)，且年相互獨(dú)立。第6年初的各年齡魚數(shù)與平衡捕撈時(shí)的年初值刈(j)(j=1,2,3,4)的偏離不超過工。的乜倍1186*042Max f(k尸& = ij 也 42r + 0.8工式i T) 1

17、 一 eI 0.42匕) Tx2(0 = xi(i-l)e-0i8 x3(i) = q。- 1)6一0總s.tx4(0 = X4(i- l)e15* 1109 * 10(yO 工3(1 - )C 34-2k3* L109* 10JM1 =而Kt(iT)e 0.84k-58 + - L1.6尸_b081.61.22 * 10n(/V + /V4)*1 - L22 * IO11 + (Na +/V4)1(5)-%o(O I 0(D用Matlab進(jìn)行求解(代碼見附錄),結(jié)果可得：口=13.0018 七=13.9406%=21.6540 儲(chǔ)=21.5971=17.2248M=0.9989N 匕=0.9

18、500=3 = 0.9500=0.9500此時(shí)，最高總收獲量為1.6299“噸六、模型的評(píng)價(jià)與改進(jìn)6.1模型的優(yōu)點(diǎn)1）、每個(gè)模型的數(shù)值都可以通過 Matlab計(jì)算解得，且解出來的值與實(shí)際情況相吻合， 都可以用一般常識(shí)來解釋，解的穩(wěn)定性較好。2）、本題用到的一些數(shù)學(xué)方法及知識(shí)，比如說常微分方程求解、一元函數(shù)極值求解等都 比較簡(jiǎn)單。3）、本模型通過了解魚群數(shù)量關(guān)系，運(yùn)用連續(xù)性隨機(jī)變量的方法并建立微分方程解決問 題，使魚群死亡率在模型中更接近實(shí)際情況，實(shí)用性較強(qiáng)。4）、本模型的建立實(shí)現(xiàn)了可持續(xù)捕撈，在獲得最高收獲量時(shí)還能有效地保持生產(chǎn)平衡。5）、沒有單一性，可以預(yù)測(cè)情況相差不大的其他物種之間的關(guān)系和

19、發(fā)展趨勢(shì)。6.2模型的缺點(diǎn)1）本題沒有考慮社會(huì)資源的限制，這與實(shí)際情況不相符。2）沒有考慮禁捕期內(nèi)的3、4齡魚的產(chǎn)卵過程的連續(xù)性，造成了一定的誤差。3）忽略了一些其他因素對(duì)漁場(chǎng)的影響，得到的結(jié)果與實(shí)際有一定的差距。5.2模型的改進(jìn)1）該模型雖然提出了產(chǎn)卵與孵化過程應(yīng)該是連續(xù)的過程，但在求解過程中實(shí)質(zhì)上還是利用了離散的過程代替了連續(xù)過程。因此可以在此方向上進(jìn)行改進(jìn)，可以分成兩種情況： 1、每個(gè)時(shí)刻t, 3、4齡魚都可以固定的概率產(chǎn)卵；2、每個(gè)時(shí)期t, 3、4齡魚都以變概 率產(chǎn)卵，這樣又可以引入兩個(gè)模型。2）對(duì)4齡魚在每年年末的情況還可以再分一種情況，設(shè)其為與前面無關(guān)的第5類魚，再根據(jù)有關(guān)的資料查

20、出產(chǎn)卵量與年齡的關(guān)系，得到一個(gè)數(shù)學(xué)表達(dá)式，對(duì)模型進(jìn)行適當(dāng)?shù)?改進(jìn)。七、參考文獻(xiàn)1孫祥,徐流美，吳濤編，MATLAB7.0礎(chǔ)教程清華大學(xué)出版社 2005.2吳振奎 王全文主編運(yùn)籌學(xué)中國(guó)人民大學(xué)出版社2005.3羅君，劉鵬周，鳴煒，最優(yōu)捕魚策略模型，數(shù)學(xué)的時(shí)間與認(rèn)識(shí)， 27卷： 51-27,1997.4朱元鼎，王文浜，中國(guó)動(dòng)物圖譜（魚類） ，科學(xué)出版社，北京， 1973.5E.C皮洛，數(shù)學(xué)生態(tài)學(xué)引論，科學(xué)出版社，北京，1978.附錄微分方程求解過程：d#i（）lr=-axdc）dm丁二”式。dxj（t）dt = 3（。- 042 也打I dx3（t）t G OJdX4（t）心 =一%3 一如hdt

21、二一 ax4(t)t e 04 r 2 t e o,-T到 ，e l ,g I1由（1）式兩邊對(duì)t積分J：。篇：一乳匕得K|（i）=a）u.同理對(duì)（2）式求積分，結(jié)果得乂3=式2）（3）式對(duì)t從0到積分O；o翳=w（-0.42女-己）比,得心）=3（0）廣出二對(duì)t E10。積分，得到&=xo（3）e 口向5（4）式對(duì)t從1到i積分，得到及=如已 s ）同上分別對(duì)（5）、（6）式積分，得I 4 I 一）=:x4（t） = xo（4）e+X(4)=孫(4)e 一式由模型假設(shè)，魚產(chǎn)卵數(shù)n=1.109*l。*酋（9+ MJ.2(0 42k f a)EI ,=1.109*1.22* 10”因?yàn)榈诙昴瓿?/p>

22、1齡魚的數(shù)量n五.內(nèi)1.因?yàn)?；=.II，所以得到1=1.109* *-；.-)問題一代碼(第一種情況)：%最優(yōu)捕魚策略ch431%件名：ch431.mx=sym( x);E3=0.42*x;% 捕撈率d=1.22*10A11;r=0.8;q=d*exp(-(3*r+2/3*E3)*(32529.55*exp(r)+65059.1*exp(-2/3*x)/(1-exp(-(r+2/3*x)；N10=d*q/(d+q);N40=d*q/(d+q)*exp(-(3*r+2/3*E3)/(1-exp(-(r+2/3*x);a3=E3/(r+E3)*(1-exp(-2/3*(r+E3)*d*q/(d+

23、q)*exp(-2*r);a4=x/(r+x)*(1-exp(-2/3*(r+x)*d*q/(d+q)*exp(-(3*r+2/3*E3)/(1-exp(-(r+2/3*x)；M3=17.86;M4=22.99;M=M3*a3+M4*a4;M1=-M;M10=char(M);M11=char(M1);fplot(M10,0,100)E4=fminbnd(M11,0,100);E3=0.42*E4;d=1.22*10A11;r=0.8;q=d*exp(-(3*r+2/3*E3)*(32529.55*exp(r)+65059.1*exp(-2/3*E4)/(1-exp(-(r+2/3*E4);N1

24、0=d*q/(d+q);N40=d*q/(d+q)*exp(-(3*r+2/3*E3)/(1-exp(-(r+2/3*E4);a3=E3/(r+E3)*(1-exp(-2/3*(r+E3)*d*q/(d+q)*exp(-2*r);a4=E4/(r+E4)*(1-exp(-2/3*(r+E4)*d*q/(d+q)*exp(-(3*r+2/3*E3)/(1-exp(-(r+2/3*E4);M3=17.86;M4=22.99;M=M3*a3+M4*a4;E3E4Max=M%第一種情況，四齡魚過一年全掛了N1=N10 %各年齡組魚的數(shù)量N2=N1*exp(-r)N3=N2*exp(-r)N4=N40問

25、題二代碼(情況一、情況二) ：function fishf()global x0x=zeros(6,4);r1=0;r2=100;r,v1=fminbnd(fish1,r1,r2);% 解第一問ee=exp(-0.8);x10=f1(r);x0(1)=x10;x0(2)=x10*ee;x0(3)=x10*ee*ee;x0(4)=x10/(1-exp(-0.8-2*r/ 3)*exp(-2.4-0.84*r/ 3);rw1=-v1x0pause;5 年的 r 不變時(shí)5年 r 不同時(shí)r2,v2=fmincon(fish2,10,0,100,fish2con);% 解第二問， r2w2=-v2x=f

26、2(r2);x_x0=x(6,:)./x0pause;r0=10 10 10 10 10;vlb=0 0 0 0 0;vub=150 150 150 150 150;r3,v3=fmincon(fish3,r0,vlb,vub,fish3con);% 解第二問，r3w3=-v3x=f3(r3);x_x0=x(6,:)./x0 function f=fish1(r)c1=0.8/(3*1.109e5*(exp(-2.4-2.84*r/ 3)/(1-exp(-0.8-2*r/ 3); c2=0.5*exp(-1.6-0.84*r/ 3)+(exp(-2.4-2.84*r/ 3)/(1-exp(-0

27、.8-2*r/ 3);c3=17.86*0.42/(0.42*r+0.8)*exp(-1.6)*(1-exp(-0.8-2*r/ 3)*(1-exp(-1.6-2*r/ 3); c4=22.99/(r+0.8)*(exp(-2.4-0.84*r/ 3)/(1-exp(-0.8-2*r/ 3)*(1-exp(-1.6-2*r/ 3); f=-1.22e9*r*(1-c1 /c2)*(c3+c4);returnfunction f=f1(r)c1=0.8/(3*1.109e5*(exp(-2.4-2.84*r/ 3)/(1-exp(-0.8-2*r/ 3);c2=0.5*exp(-1.6-0.84

28、*r/ 3)+(exp(-2.4-2.84*r/ 3)/(1-exp(-0.8-2*r/ 3);f=1.22e11*(1-c1 /c2);return function x=f2(r)global x0ee=exp(-0.8);x(1,:)=122 29.7 10.1 3.29*1e9;for i=1:5,n=3*1.109e5/0.8*(exp(-1.6/3)-exp(-0.8);n3=n*x(i,3)*exp(-0.84*r/ 3);n4=n*x(i,4)*exp(-2*r/ 3);x(i+1,1)=1.22e11*(0.5*n3+n4)/(1.22e11+0.5*n3+n4);x(i+1

29、,2)=x(i,1)*ee;x(i+1,3)=x(i,2)*ee;x(i+1,4)=x(i,4)*exp(-0.8-2*r/ 3)+x(i,3)*exp(-0.8-0.84*r/ 3);endreturn function g,eg=fish2con(r)global x0ee=exp(-0.8);x(1,:)=122 29.7 10.1 3.29*1e9;for i=1:5,n=3*1.109e5/0.8*(exp(-1.6/3)-exp(-0.8);n3=n*x(i,3)*exp(-0.84*r/ 3);n4=n*x(i,4)*exp(-2*r/ 3);x(i+1,1)=1.22e11*(

30、0.5*n3+n4)/(1.22e11+0.5*n3+n4);x(i+1,2)=x(i,1)*ee;x(i+1,3)=x(i,2)*ee;x(i+1,4)=x(i,4)*exp(-0.8-2*r/ 3)+x(i,3)*exp(-0.8-0.84*r/ 3); enda=0.05g=(1-a)*x0-x(6,:)-(1+a)*x0+x(6,:);eg=;return function g,eg=fish3con(r)global x0ee=exp(-0.8);x(1,:)=122 29.7 10.1 3.29.*1e9;for i=1:5,n=3*1.109e5/0.8*(exp(-1.6/3)

31、-exp(-0.8);n3=n*x(i,3)*exp(-0.84*r(i)/ 3);n4=n*x(i,4)*exp(-2*r(i)/ 3);x(i+1,1)=1.22e11*(0.5*n3+n4)/(1.22e11+0.5*n3+n4);x(i+1,2)=x(i,1)*ee;x(i+1,3)=x(i,2)*ee;x(i+1,4)=x(i,4)*exp(-0.8-2*r(i)/ 3)+x(i,3)*exp(-0.8-0.84*r(i)/ 3); enda=0.05g=(1-a)*x0-x(6,:)-(1+a)*x0+x(6,:);eg=;return function x=f3(r)global x0ee=exp(-0.8);x(1,:)=122 29.7 10.1 3.29.*1e9;for i=1:5,n=3*1.109e5/0.8*(exp(-1.6/3)-exp(-0.8);n3=n*x(i,3)*exp(-0.84*r(i)/ 3);n4=n*x(i,4)*exp(-2*r(i)/ 3);x(i+1,1)=1.22e11*(0.5*n3+n4)/(1.22e11+0.5*n3+n4);x