1、2011高教社杯全國大學生數學建模競賽編 號 專 用 頁 賽區評閱編號（由賽區組委會評閱前進行編號）：賽區評閱記錄（可供賽區評閱時使用）：評閱人評分備注全國統一編號（由賽區組委會送交全國前編號）：全國評閱編號（由全國組委會評閱前進行編號）：核反應堆屏蔽層設計摘要核反應堆，又稱為原子反應堆或反應堆，是裝配了核燃料以實現大規模可控制裂變鏈式反應的裝置。為防護中子、射線和熱輻射，必須在反應堆和大多數輔助設備周圍設置屏蔽層。通常設計最外層屏蔽時應將輻射減到人類允許劑量水平以下，常稱為生物屏蔽。核反應堆屏蔽層是用一定厚度的鉛把反應堆四周包圍起來，用以阻檔或減弱反應堆發出的各種射線。在各種射線中，中子對人

2、體傷害極大，因此，屏蔽設計，主要是了解中子穿透屏蔽的百分比(或概率)，這對反應堆的安全運行是至關重要的。問題一在假定屏蔽層厚度已知的情況下，要求求解中子穿透屏蔽層的百分比。采用傳統的的蒙特卡羅法對問題進行分析，利用軟件計算出當屏蔽層時的時的概率。首先假定屏蔽層是理想的均勻鉛質平板；中子在屏蔽層內相繼兩次碰撞之間游動的距離服從指數分布；兩次碰撞之間中子的平均游動距離為；中子在軸上的位置表示為，中子運動方向與軸夾角為，即用描述第次碰撞后的中子的運動狀態；中子碰撞后的彈射角服從的均勻分布；從而得到第次碰撞后在屏蔽層的位置。由假設可能得到式中為兩次碰撞之間中子的平均游動距離；是區間上均勻分布的隨機數。

3、然后采用逆轉法求解平均分布和指數分布的隨機數。問題二用蒙特卡羅模擬要求穿透率為時，需試驗次。為節省時間，減少計算，在問題一上做改進。采用蒙特卡羅模擬收斂法，先用分析法求出穿透率為時屏蔽層厚度，則三倍此厚度即可達到穿透率不大于的要求。再用計算機收索法進行驗證，得出相同結果且符合實際要求。最后根據上述估計為日本核危機善后提出建議。關鍵詞：蒙特卡羅模擬 分析法 計算機收索法一、問題重述1.1問題背景核反應堆，又稱為原子反應堆或反應堆，是裝配了核燃料以實現大規模可控制裂變鏈式反應的裝置。為防護中子、射線和熱輻射，必須在反應堆和大多數輔助設備周圍設置屏蔽層。有些屏蔽材料俘獲中子后放射出射線，因此在中子屏

4、蔽外要有一層射線屏蔽。通常設計最外層屏蔽時應將輻射減到人類允許劑量水平以下，常稱為生物屏蔽。核反應堆屏蔽層是用一定厚度的鉛把反應堆四周包圍起來，用以阻檔或減弱反應堆發出的各種射線。在各種射線中，中子對人體傷害極大，因此，屏蔽設計，主要是了解中子穿透屏蔽的百分比(或概率)，這對反應堆的安全運行是至關重要的。1.2問題提出假設屏蔽層，在大數定理的意義下，中子穿透屏蔽層的百分比多少？在實際應用中，要求中子穿誘屏蔽層的概率極小，一般數量級為，即穿入屏蔽層的中子若為幾百萬個，也只能有幾個中子穿過屏蔽層。問題是多厚的屏蔽層才能使它被穿的概率小于?根據上述估計，并查閱相關文獻，嘗試為日本福島核泄漏事件的核危

5、機善后工作提出約2000字的建議。二、模型假設假定屏蔽層是理想的均勻鉛質平板；假定中子在屏蔽層內相繼兩次碰撞之間游動的距離服從指數分布；兩次碰撞之間中子的平均游動距離為；中子在軸上的位置表示為，中子運動方向與軸夾角為；中子碰撞后的彈射角服從的均勻分布；3、 符號說明：精確度 ：單層屏蔽厚度：試驗次數 ：兩次碰撞間中子游動距離：穿透率 ：單層穿透率 ：隨機數 ：為中子在軸上的位置：中子運動的方向和軸的夾角4、 模型的建立與求解4.1問題一問題分析首先考慮一個中子進入屏蔽層后運動的物理過程：假定屏蔽層是理想的均勻鉛質平板，中子以初速度和方向角射入屏蔽層內，如圖1所示，運動一段距離后，在處與鉛核碰撞

6、之后，中子獲得新的速度及方向，再運動一段距離后，與鉛核第二次碰撞，并獲得新的狀態等等，經若干次碰撞后，發生以下情況之一則終止運動過程：(1)彈回反應堆；(2)穿誘屏蔽層；(3)第i次碰撞后，中子被屏蔽層吸收。圖1蒙特卡羅法作為一種計算方法，是由S.M.烏拉姆和J.馮諾伊曼在20世紀40年代中葉為研制核武器的需要而首先提出來的。在此之前，該方法的基本思想實際上早已被統計學家所采用了。例如，早在17世紀，人們就知道了依頻數來決定概率的方法。采用傳統的的蒙特卡羅法對問題進行分析，利用軟件計算出當屏蔽層時的概率。使用蒙特卡羅方法進行分子模擬計算是按照以下步驟進行的： 使用隨機數發生器產生一個隨機的分子

7、構型。對此分子構型的其中粒子坐標做無規則的改變，產生一個新的分子構型。計算新的分子構型的能量。 比較新的分子構型于改變前的分子構型的能量變化，判斷是否接受該構型。 若新的分子構型能量低于原分子構型的能量，則接受新的構型，使用這個構型重復再做下一次迭代。若新的分子構型能量高于原分子構型的能量，則計算玻爾茲曼常數，同時產生一個隨機數。 若這個隨機數大于所計算出的玻爾茲曼因子，則放棄這個構型，重新計算。 若這個隨機數小于所計算出的玻爾茲曼因子，則接受這個構型，使用這個構型重復再做下一次迭代。 如此進行迭代計算，直至最后搜索出低于所給能量條件的分子構型結束。模型的建立假定屏蔽層平行板得厚度為，其中為兩

8、次碰撞之間中子的平均游動距離；每次碰撞后中子因損失一部分能量而速度下降，假設第十次碰撞后，中子的速度下降到某一很小的數值而終止運動（被吸收）。由于對穿透屏蔽層的中子感興趣，故用描述第次碰撞后的中子的運動狀態，其中為中子在軸上的位置，中子運動的方向和軸的夾角。中子在再屏蔽層內相繼兩次碰撞之間游離的距離服從指數分布，中子碰撞后的彈射角服從上的均勻分布。從而得到第次碰撞后在屏蔽層的位置 式中，是中子第次碰撞后的彈射角度，是中子從第次碰撞到第次碰撞時所游動的距離。由假設可能得到：式中為兩次碰撞之間中子的平均游動距離；是區間上均勻分布的隨機數。式表明了中子運動的位置和方向是隨機的。生成連續隨機變數的方法

9、很多。選擇哪些算法取決于我們所希望生成變數的的分布，需要考慮到隨機變量的準確性、計算和存儲效率以及算法的復雜性等因素。最常用的兩種算法是逆轉法和接受-排除法。下面采用逆轉法求解平均分布和指數分布的隨機數：設指數分布的密度函數為：分布函數為：生成隨機數，設求解得。由此可得重新整理后得：取方程兩邊的自然對數，可得最后求解可得為了簡化計算用代替。由于是隨機數，所以也是隨機數，因此指數分布的過程生成器將是：分析在區間上均勻分布的隨機變量，該分部的由下列函數給出：用法從這個隨機變量生成觀測值，分布函數求得：要使用從均勻分布生成觀測值，首先生成隨機數，然后設以求解。由此可得：求解可得：它是均勻分布的生成器

10、。可以用下述公式對進行模擬：式中：是均勻分布隨機數。中子在屏蔽層里隨機游動，類似布朗運動，第次碰撞后，中子的位置有三種情況：，中子返回反應堆；，中子穿出屏蔽層；若，則中子在屏蔽層內繼續運動；若，則中子被屏蔽層吸收。中子的三種運動模式如圖2。針對中子的三種運動模式，利用程序模擬求解，模擬結果如下表1：中子數（個）返回（%）吸收（%）穿透（%）10082.00%2.00%12.00%100081.60%5.80%12.60%300081.036.73%12.23%500081.28%6.30%12.42%1000080.98%6.47%11.40%1500081.69%6.30%11.97%200

11、0081.50%5.83%12.67%10000081.75%6.18%12.07%中子擴散圖像如下：中子穿透屏蔽層的百分比約為12%。4.2問題二4.2.1模型的建立從表1可知，當要求精確度時，試驗次數。在實際應用中，要求穿透概率很小，比如，即是說，穿入屏蔽層的中子數若為幾百萬個，只允許幾個穿透屏蔽層。此時，要達到千分之一的精度，則試驗次數為，要獲得更高的精度，則問題一的模擬法已行不通。現在，討論加快模擬收斂的方法。把均勻平行板分為厚度相同的層，只取一層進行模擬。設中子在一層中被吸收和彈回的概率之和為，則穿過一層的概率是，因此穿透層的概率是。由于中子穿過一層的平均速度有所下降，因而總的穿透概

12、率。對方法，試驗次數可保證的精度好于，經層后有，若，就可獲得穿透概率。這樣處理后，不必作高達的試驗，只需作次的試驗就可達到的精度。這一改進比直接的模擬方法大大加快了收斂速度，從而大大減少了模擬時間。4.2.2模型求解 分析法求解：由于模型模擬過程中設計較多中子，又每個中子碰撞最多達次，因此模擬萬個中子，則實際運行最多要處理萬次碰撞。通過循環語句逐個處理則費時，因此采用加快模擬收斂的方法，可以利用的功能，編寫成值文件，提高模擬效率。模擬得當時模擬的結果如下表所示：中子數/個返回/%穿透/%吸收/%100080.413.16.550008112.786.221000081.612.116.29表格

13、 2：模擬結果圖表 1：模擬結果圖示 模擬結果表明取屏蔽層厚不合適，中子穿透屏蔽層的百分比超過了，在實際應用中要求穿透屏蔽的概率很小，數量級為，即進入屏蔽層的中子若為百萬個以上，只允許幾個穿過屏蔽層。設計屏蔽層的厚度，可把屏蔽層看成層厚度均為的平行板，由于中子穿過一層后的平均速度有所下降，若記中子穿過厚度屏蔽層的概率為，則穿過整個屏蔽的概率有：利用前面的模擬結果：時，令，即：即：因此取蔽層的厚度時，可使中子穿透屏蔽層的概率。經過三層屏蔽后，則中子穿透概率，整個屏蔽層的厚度為。計算機搜索法驗證：逐步增大屏蔽層厚度，如：屏蔽層厚度開始取，若則增加厚度為，直到穿透概率為止。經過三層屏蔽后，則可使中子

14、穿透屏蔽層的概率。運用模擬結果如下：厚度返回/%穿透/%吸收/%82.2211.656.1381.982.4915.5381.760.817.4482.81017.19表格 3：模擬結果圖表 2：模擬結果圖示驗證結果與分析法相符合，因此，經過三層屏蔽后，整個屏蔽層的厚度為，則中子穿透概率。5、 建議報告2011年3月11日日本本州島附近海域發生強烈地震，強震導致福島第一核電站發生爆炸引發核危機，從這一刻開始就引發了核泄漏。福島第一核電站的一號、三號和二號反應堆分別在當地時間日時分、日時分和日晨時分發生氫氣爆炸。福島發生的核泄漏事件死亡人、失蹤人，輻射泄漏最初一星期內向大氣釋放的輻射物質活度高達

15、至萬萬億貝克勒爾，對海洋領域造成的損失不可估計。福島核電站從選址到設計都有嚴格質量監控，對抗震也有很高的技術要求，福島核電站發生核泄漏的可能原因是日本遭受的特大地震超過了當初的抗震標準。另一方面核反應堆在停堆后仍需要對堆芯進行冷卻，因為核燃料有自衰變余熱，若長時間得不到冷卻，會使得堆芯溫度上升，導致核燃料棒熔化。一般停堆后，冷卻系統泵所需的電力需要從外部輸入，同時備有多臺應急發電機供電。但日本強震后，福島核電廠外電網癱瘓，應急柴油發電機短暫運行后也陷入中斷狀態。為防止核反應堆內大量放射線從堆芯穿透出來，一般反應堆設有三到四道屏障。第一道是把核燃料放在陶瓷芯塊中，使得大部分裂變產物和氣體產物保存

16、在芯塊中；第二道是合金制造的包殼，把陶瓷芯塊密封在其中，使得其在高溫下不與水發生反應；第三道是耐壓力容器，避免放射性物質泄漏到反應堆廠房；第四道是反應堆安全殼，一般是厚近厘米的鋼筋混凝土墻，里面襯有防輻射金屬材料，防止放射性物質進入環境。安全殼是最重要的一道屏障，以防前面的措施失效，它能保證周圍居民遭受輻射在安全范圍內。核反應堆，又稱為原子反應堆或反應堆，是裝配了核燃料以實現大規模可控制裂變鏈式反應的裝置。制裂變鏈式反應的裝置。反應堆主要有活性區，反射層，外壓力殼和屏蔽層組成。活性區又由核燃料，慢化劑，冷卻劑和控制棒等組成；反射層是為減少堆芯的中子泄漏、降低臨界質量和臨界尺寸、提高并展平堆芯中

17、子注量率而布置在堆芯四周的材料或物體。應具有大的中子散射截面和小的中子吸收截面；屏蔽層是為防護中子、射線和熱輻射，必須在反應堆和大多數輔助設備周圍設置屏蔽保護層。其中慢化劑、控制棒、冷卻劑、屏蔽層是核心組件。慢化劑就是用來將快中子能量減少，使之慢化成為中子或中能中子的物質，常采用石墨、重水、輕水；控制棒在反應堆中起補償和調節中子反應性以及緊急停堆的作用；冷卻劑由主循環泵驅動，在一回路中循環，從堆心帶走熱量并傳給二回路中的工質，使蒸汽發生器產生高溫高壓蒸汽，以驅動汽輪發電機發電，是既在堆心中工作又在堆外工作的一種反應堆成分，這就要求冷卻劑必需在高溫和高中子通量場中工作是穩定的。屏蔽層是為防護中子

18、、射線和熱輻射，必須在反應堆和大多數輔助設備周圍設置屏蔽保護層，其設計要力求造價便宜并節省空間。目前日本福島核電站的放射性物質控制在安全殼內，仍未造成真正意義上的核擴散，但一旦安全殼密閉失效將會造成更為嚴重的后果。根據問題一、問題二的求解可知核反應堆的屏蔽層設計的厚度嚴重影響了中子穿透屏蔽層的百分比。當屏蔽層的厚度（為兩次碰撞間中子的平均游動距離）時，中子穿透屏蔽層的百分比為，遠遠大于實際應用中要求穿透屏蔽層的概率。在實際應用中要求穿透屏蔽層的概率的數量級為，即進入屏蔽層的中子若為幾百萬個以上，只允許幾個穿透屏蔽層。于是用蒙特卡羅方法把屏蔽層看成層厚度均為的平行板，得到當屏蔽層的厚度時，可是穿

19、透屏蔽層的概率達到實際應用要求。可見增加屏蔽層的厚度可以在很大程度上減小放射性物質的泄露。一旦發生核泄漏，應及時用人工方式使反應堆冷卻，通過注入海水是堆芯降溫；同時應采取措施防止冷卻水的流失，封閉核污染區，保證人生安全和環境清潔，及時疏散周邊人員。日本政府應加大對1、2、3、4號機組的監測與處理工作，以反應堆安全殼的破壞造成更大規模的核泄漏和更嚴重的損失。年發生的史上最大的切爾諾貝利核電站核泄漏事故導致事故后前個月內有人死亡，之后年內有萬人死亡，萬人遭受各種程度的輻射疾病折磨，方圓公里地區的萬多民眾被迫疏散；年日本福島發生的核泄漏事件死亡人、失蹤人，輻射泄漏最初一星期內向大氣釋放的輻射物質活度

20、高達至萬萬億貝克勒爾，對海洋領域造成的損失不可估計。核泄漏造成的嚴重損失以及對人類生命財產的致命威脅，為我們敲響了要和平利用核能源與預防核泄漏的警鐘。如何有效地預防核泄漏以及在核泄漏發生后應采取怎樣的措施使核泄漏造成的威脅和傷害達到最低，這就需要我們對核反應堆、核泄漏發生的原因和核輻射有深入的了解。通過福島核泄露的原因分析，需要加強核反應堆的防震功能來應對地震對核反應堆的破壞；同時要設有多臺應急發電機，以防地震造成的電網癱瘓而是反應堆堆芯無法冷卻；加強反應堆停堆系統、應急冷卻系統和安全殼等安全措施的發展是其在發生核泄漏時能發揮應有的作用，從而保證核能的安全開發和利用。六、模型評價6.1模型評價

21、6.1.1優點 問題一，蒙特卡羅法作為一種計算方法，是由S.M.烏拉姆和J.馮諾伊曼在20世紀40年代中葉為研制核武器的需要而首先提出來的。在此之前，該方法的基本思想實際上早已被統計學家所采用了。例如，早在17世紀，人們就知道了依頻數來決定概率的方法。屬于比較成熟的分析方法，所得結果比較可靠。問題二，有些系統構造雖不復雜，卻屬內部機理不明的“黑箱”系統，我們無法用機理分析建立模型，系統模擬是很最重要的一種方法。每模擬一次中子的運動就相當于對中子運動過程做一次“試驗”或“觀察”，按理論來講必須模擬無限次，才能求出中子穿出屏蔽層的概率。在實際中常進行較大次數的模擬，用所得頻率作為所求概率的近似解。

22、設計模擬模型是為了“觀察”系統并收集必要信息，這和觀察一個實現系統非常相似，模擬系統實現系統的仿真。6.1.2缺點 模擬模型本質是一種求解問題的試驗方法，是對思維結果的一種驗證。由于有隨機數每次運行結果不同，與實際有偏差。方法比較單一，運行的精確度雖然已經很高，但是對于實際情況來說仍然需要提高，特別是在當今倡導人權的社會中，需要給予核設施附近的工作人員更好的保護，因此需要將精確度進一步提高。七、參考文獻1阮曉青，周義倉，數學建模引論，北京：高等教育出版社，2005年。2袁新生，邵大宏，郁時煉，北京：科學出版社，2007年。3呂金虎，陸君安，陳士華，武漢：武漢大學出版社，2002年。4李乃文，北

23、京：清華大學出版社，2006年。八、附錄 附錄一：值文件function sim_zhongzi2%模擬核反應堆屏蔽層中中子的運動%為了提高程序運行速度，這里沒有逐個動態顯示中子的運動，%所以一開始就產生所有中子的初始狀態；當然，也可以很容易修改為動態顯示程序。%n=input(中子個數：);%中子個數xx=;%存儲第10次碰撞后中子的橫坐標yy=;%存儲第10次碰撞后中子的縱坐標 N=10;%每個中子最多碰撞次數d = 2;D = 3*d;H =10*10;c=zeros(1,3);%c(1)返回反應堆數量，c(2)穿透,c(3)吸收 %產生中子 x=zeros(1,n); y=H*rand

24、(1,n)/10; %中子運動 for j=1:N,%考慮至多10次碰撞 if isempty(x), break; end R = -d*log(rand(1,length(x); seta = 2*pi*rand(1,length(x); x = x + (R.*cos(seta); y = y + (R.*sin(seta); %c(2)穿透 t=find(xD); c(2) = c(2) + length(t); if length(t)0, xx=xx,x(t); yy=yy,y(t); x(t)=; y(t)=; end %c(1)返回反應堆數量 t=find(x0, xx=xx,

25、x(t); yy=yy,y(t); x(t)=; y(t)=; end if j=N, c(3) = c(3) + length(x); end end xx=xx,x; yy=yy,y; check = sum(c) - n %right-0 check2=length(xx)-n%right-0 bili=c/n %return%將此行注釋或刪除后才能執行下列代碼 %重繪所有中子，當n增大時，會花較長時間 plot(xx,yy,r.) hold on line(0,0,-H,H) hold on line(D,D,-H,H) text(-4*d,H-5*d,返回) text(-4*d,H-

26、9*d,sprintf(%6.2f%),c(1)/n*100) text(D/2-d,H-5*d,吸收) text(D/2-2*d,H-9*d,sprintf(%6.2f%),c(3)/n*100) text(D+2*d,H-5*d,穿透) text(D+d,H-9*d,sprintf(%6.2f%),c(2)/n*100) hold off附錄二：程序subplot(2,2,1)n=10000;xx=;%存儲第10次碰撞后中子的橫坐標yy=;%存儲第10次碰撞后中子的縱坐標N=10;%每個中子最多碰撞次數d = 2;D = 3*d;H =10*10;c=zeros(1,3);%c(1)返回反

27、應堆數量，c(2)穿透,c(3)吸收 %產生中子 x=zeros(1,n); y=H*rand(1,n)/10; %中子運動 for j=1:N,%考慮至多10次碰撞 if isempty(x), break; end R = -d*log(rand(1,length(x); seta= 2*pi*rand(1,length(x); x= x + (R.*cos(seta); y= y + (R.*sin(seta); %c(2)穿透 t=find(xD); c(2) = c(2) + length(t); if length(t)0, xx=xx,x(t); yy=yy,y(t); x(t)

28、=; y(t)=; end %c(1)返回反應堆數量 t=find(x0, xx=xx,x(t); yy=yy,y(t); x(t)=; y(t)=; end if j=N, c(3) = c(3) + length(x); end end xx=xx,x; yy=yy,y; check = sum(c) - n %right-0 check2=length(xx)-n%right-0 bili=c/n %return%將此行注釋或刪除后才能執行下列代碼 %重繪所有中子，當n增大時，會花較長時間 plot(xx,yy,r.) hold on line(0,0,-H,H) hold on lin

29、e(D,D,-H,H) text(-4*d,H-5*d,返回) %text(-4*d,H-9*d,sprintf(%6.2f%),c(1)/n*100) text(D/2-d,H-5*d,吸收) %text(D/2-2*d,H-9*d,sprintf(%6.2f%),c(3)/n*100) text(D+2*d,H-5*d,穿透) %text(D+d,H-9*d,sprintf(%6.2f%),c(2)/n*100) hold offsubplot(2,2,2)n=10000;xx=;%存儲第10次碰撞后中子的橫坐標yy=;%存儲第10次碰撞后中子的縱坐標N=10;%每個中子最多碰撞次數d =

30、 2;D = 6*d;H =10*10;c=zeros(1,3);%c(1)返回反應堆數量，c(2)穿透,c(3)吸收 %產生中子 x=zeros(1,n); y=H*rand(1,n)/10; %中子運動 for j=1:N,%考慮至多10次碰撞 if isempty(x), break; end R = -d*log(rand(1,length(x); seta= 2*pi*rand(1,length(x); x= x + (R.*cos(seta); y= y + (R.*sin(seta); %c(2)穿透 t=find(xD); c(2) = c(2) + length(t); if

31、 length(t)0, xx=xx,x(t); yy=yy,y(t); x(t)=; y(t)=; end %c(1)返回反應堆數量 t=find(x0, xx=xx,x(t); yy=yy,y(t); x(t)=; y(t)=; end if j=N, c(3) = c(3) + length(x); end end xx=xx,x; yy=yy,y; check = sum(c) - n %right-0 check2=length(xx)-n%right-0 bili=c/n %return%將此行注釋或刪除后才能執行下列代碼 %重繪所有中子，當n增大時，會花較長時間 plot(xx,

32、yy,r.) hold on line(0,0,-H,H) hold on line(D,D,-H,H) text(-4*d,H-5*d,返回) %text(-4*d,H-9*d,sprintf(%6.2f%),c(1)/n*100) text(D/2-d,H-5*d,吸收) %text(D/2-2*d,H-9*d,sprintf(%6.2f%),c(3)/n*100) text(D+2*d,H-5*d,穿透) %text(D+d,H-9*d,sprintf(%6.2f%),c(2)/n*100) hold offsubplot(2,2,3)n=10000;xx=;%存儲第10次碰撞后中子的橫

33、坐標yy=;%存儲第10次碰撞后中子的縱坐標N=10;%每個中子最多碰撞次數d = 2;D = 8*d;H =10*10;c=zeros(1,3);%c(1)返回反應堆數量，c(2)穿透,c(3)吸收 %產生中子 x=zeros(1,n); y=H*rand(1,n)/10; %中子運動 for j=1:N,%考慮至多10次碰撞 if isempty(x), break; end R = -d*log(rand(1,length(x); seta= 2*pi*rand(1,length(x); x= x + (R.*cos(seta); y= y + (R.*sin(seta); %c(2)穿

34、透 t=find(xD); c(2) = c(2) + length(t); if length(t)0, xx=xx,x(t); yy=yy,y(t); x(t)=; y(t)=; end %c(1)返回反應堆數量 t=find(x0, xx=xx,x(t); yy=yy,y(t); x(t)=; y(t)=; end if j=N, c(3) = c(3) + length(x); end end xx=xx,x; yy=yy,y; check = sum(c) - n %right-0 check2=length(xx)-n%right-0 bili=c/n %return%將此行注釋或刪除后才能執行下列代碼 %重繪所有中子，當n增大時，會花較長時間 plot(xx,yy,r.) hold on line(0,0,-H,H) hold on line(D,D,-H,H) text(