波前法及matlab實現

1、有限元二維熱傳導波前法MATLAB程序 · 二維熱傳導有限元· 使用高斯消去法解線性方程組的二維熱傳導有限元程序· 波前法的基本概念與算法· 使用波前法解線性方程組的二維熱傳導有限元程序· 消元過程· 波前法與高斯消去法的效率之比較· 小結：波前法的過去、現在和未來波前法是求解線性方程組的一種方法，廣泛用于有限元程序。它最初由英國人（？）B.M. Irons于1970在“國際工程計算方法雜志”上發(fā)表。30多年來，波前法有了不少變種。本文所用算法，采于法國人Pascal JOLY所著Mise en Oeuvre de

2、 la Méthode des Eléments Finis。這本書是我1993年在比利時一家書店買的，書中有一節(jié)"波前法"，六頁紙，解釋了基本概念和算法，但沒有程序，也沒有細節(jié)討論。我曾花了兩個半天的時間，在網上尋找波前法程序，或更詳細的資料，沒有找到（需要花錢才能看的文獻不算）。倒是看到不少中國人，也在尋找。一些人說，波前法程序太難懂了。 通過自己編寫程序，我同意這些人的說法，確實難。我還真很少編如此耗費腦力的程序。完工之后，我曾對朋友老王說，有了算法，編程序還這么難，當初想出算法的人，真是了不起。 現將我對波前法的理解和編程體會

3、解說如下，供感興趣的網友參考，也為填補網絡上波前法空白。二維熱傳導有限元波前法和有限元密不可分。因而，在編寫波前法程序之前，必須有個有限元程序。為了簡化問題，最好是能解算一個節(jié)點上只有一個自由度的問題的有限元程序，而且要盡可能地簡單。我手邊現有的有限元程序都不符合這個要求。就決定先開發(fā)一個解算二維熱傳導問題的MATLAB有限元程序。 二維熱傳導問題的微分方程是其中 T 是溫度，Kx 和 Ky 分別是 x 和 y 方向上的熱傳導系數，q 是熱源。 對于這樣的比較經典的二階微分方程，如何導出有限元表達式？ 這個問題，是有限元的首要問題！我相信，許許多多學過有限元的人，

4、甚至每天都在用有限元的人，并不真的十分明了。 我自己曾經就是這樣。在我于2005年3月到巴西教書之前，我搞過20年有限元，其中有六年，還是在比利時一個專門開發(fā)有限元程序的公司里度過的。但是，如果您那時問我，任給一個二階偏微分方程，例如上述熱傳導方程，如何導出有限元表達？說老實話，不看書，我還真就不會！ 直到2006年，我教了一遍有限元后，才弄明白（惟教書才是最好的學習）。 其實簡單極了：只需將那二階偏微分方程，乘上一個任意標量函數，然后在某個有限單元上積分。 請看下列推導：即其中， e 是單元面積； 是任意標量函數。  注意在以上積分中，溫度要

5、遭受二階偏導，這很不好。有限元的精髓，在于通過分步積分，將其中一階偏導轉移到那個任意標量函數 上，這樣就可降低一階溫度偏導，改善對它的苛刻待遇。這得用到您在高等數學最后一章里學過的散度定理（Theorem of divergence）： 其中， 是面積的邊界； （反） 是梯度算子；F 和 G 是任意兩個矢量函數。   對于二維問題，上述散度定理可寫為：將此散度定理應用于（2）式中的第一項積分，令得：將此積分結果帶入（3）式，得到熱傳導偏微分方程的弱化表達（weak form）：所謂“弱化”，是指對溫度函數的可導階數要求降低了。在原熱傳導偏微分方程（1）中，溫度函數必須是二階

6、可偏導函數，而在弱化表達（6）中，溫度函數只要一階偏導就行了。 有限單元法就是以偏微分方程的弱化表達式為出發(fā)點的。 現在，到了說明那個任意標量函數，是何方神圣的時候了。 有限單元法有各種各樣的變種，而最常見的，應用最廣泛的，是所謂迦遼金法（Galerkin），即將這個任意標量函數，等同于，有限元的插值函數。迦遼金法的優(yōu)點是可以最終形成對稱勁度矩陣，從而具有較好的數值穩(wěn)定性。 我們知道，在一個有限單元中，任意一點的值，例如溫度，是用節(jié)點上的溫度表達的。對于一個四節(jié)點雙線性單元來說，設四個節(jié)點的溫度分別是T1,T2,T3,T4，則單元內任意一點的溫度 T 可

7、表達為 其中 1，2，3，4，為插值函數，也稱為形函數，定義如下：其中 和 稱為形坐標，取值區(qū)間為 -1,1。 基于式（7），對溫度的偏導數可寫為：將此二偏導數代入弱化表達式（6），該式就轉化為以節(jié)點溫度為變量的代數方程：到此，為得到原偏微分方程的有限元表達，我們只需將任意標量函數，依次取為四個插值函數，1-4，就得到四個代數方程： 注意到式(9)-(12)中下標的規(guī)律，我們可將這四個方程合并寫成矩陣的形式：使用下標表達，式（13）可進一步縮寫為：式中對應于下標 i = 14 的每一個值，下標 j 取遍14。 式（14）就是熱傳導偏微分方程（1）的四節(jié)點

8、雙線性有限元表達，也是我們接下來編寫有限元程序的出發(fā)點。  使用高斯消去法解線性方程組的二維熱傳導有限元程序這個程序是專為解一個特殊的熱傳導問題而設計的。所解問題是：已知一個無限長圓筒，內徑100毫米，外徑200毫米，筒內壁表面溫度 1°C，外壁絕熱，求圓筒截面上的溫度分布。圓筒材料各向均質，熱傳導系數為 1 （單位還得查一下，但也無所謂）。問題的解答很簡單：均布，截面各點溫度都是 1°C。為什么？因為外部絕熱，沒有熱量損失。溫度只能是均布。而內壁溫度為 1°C，所以到處都是1°C。 因為問題的幾何圖形簡單，有限元網格便容易自動生成。

9、在以下程序中，第12行至第51行，生成四節(jié)點單元的有限元網格。控制變量有兩個，Cdiv 和 Tdiv。 前者定義沿圓周分成多少單元；后者定義沿徑向即筒壁厚度方向分成多少單元。如果 Cdiv = 8，Tdiv = 2， 則所生成有限元網格如下（由第52行子程序DrawFEM畫出）：第64行使用MATLAB命令 syms 定義了兩個符號變量 ksi(即前面公式中的 ),eta()。在MATLAB中，可對符號變量進行代數運算，例如定義公式，求導，積分等。第72行就利用代數運算定義了本文前面給出的四節(jié)點單元的形函數。第76和77行分別對形函數關于 ksi 和 eta 求導。第78至第99行，計算這些導

10、數在高斯積分點的數值。本程序中，每個單元有四個高斯積分點，也就是說，在 ksi 和 eta 兩個方向上，各有二個高斯積分點。 第101至124行，根據式（14）計算單元勁度矩陣，Kelem，并將其裝配入總勁度矩陣 K。 本問題沒有熱源，所以在單元循環(huán)水平上，不需計算式（14）中的熱源積分項。 第127至139行，施加邊界條件。圓筒外壁是絕熱條件，即法向熱流等于零。在有限元中，這是自然滿足的，所以式（14）中的邊界熱流積分項為零，不必考慮。唯一需要考慮的，是圓筒內壁溫度等于 1°C 。這種溫度給定的邊界條件，在數學上叫第一類邊界條件。在有限元技術中，有各種

11、各樣的方法施加第一類邊界條件。主要是考慮提高內存效率。鑒于本程序的目的是進一步開發(fā)波前法，不需要仔細考慮如何更有效地施加這種溫度給定的邊界條件，因而所用的方法是最簡單的一種：即將內壁邊界節(jié)點的各行方程，全部換為 T = Tin (Tin = 1°C)。相應地，將這些行的主元素所占據的列，乘以Tin后，移到等號右邊。這樣處理后，就得到待解的線性方程組：K x = F。 第141行，使用本人自編的高斯消去法函數，egauss，解出為未知量 x， 也就是個節(jié)點的溫度，都等于 1。這一行，也可直接用MATLAB命令，x = K F, 取代。我使用egauss的目的，是為下一步與波前

12、法解方程比較效率。程序一：使用高斯消去法解線性方程組的二維熱傳導有限元程序  波前法的基本概念與算法有限元形成的線性方程組的系數矩陣，即剛度矩陣，是稀疏矩陣，也就是說，矩陣里含有許許多多的零元素。有限元網格節(jié)點數目越巨大，非零元素與零元素的比值越小，剛度矩陣越稀疏。用普通高斯消去法求解這樣的線性方程組，完全不考慮矩陣的稀疏性，對大量的零元素也進行加減乘除，結果浪費了大量時間。不僅如此，應用高斯消去法，因為需要將整個剛度矩陣存在計算機內存里，所需計算機內存量與有限元網格節(jié)點未知量總數成平方的關系。以熱傳導問題為例。一千個節(jié)點，光存剛度矩陣，就需內存1000 x 1000 x 8 / （

13、1024 x 1024）  = 7.8 Mb。 這還沒問題。但若要計算一萬個節(jié)點的問題，就需要 10 x 10 x 7.8 = 780 Mb 來存剛度矩陣。內存為 1 Gb的計算機已經不能計算這樣的問題了，因為微軟視窗等其它系統程序還需要內存。您當然可以開始這樣的計算。微軟視窗發(fā)現內存不夠時，會自動啟動虛擬內存，也就是把硬盤上的某一塊地方，當作內存來使用。您的計算仍然能進行。但是，您一定看不見那最終的計算結果，除非幾個月內不斷電，計算機不出毛病。原因在于，與內存相比，虛擬內存的存取時間可認為是無限長！在這種情況下，因為普通高斯消去法需要極其頻繁地使用虛擬內存，它的解算時間也就無限地延

14、長了。 但如果您在這樣的計算機上運行ANSYS，或任何需要花錢買的有限元程序，就會發(fā)現，解算同樣的問題，只需幾分鐘。您甚至可以毫無困難地解算十萬個節(jié)點的熱傳導問題。 秘密就在于，這些商業(yè)有限元軟件，在求解線性方程組時，都盡可能地利用剛度矩陣的稀疏性。波前法就是這樣一種充分考慮了剛度矩陣的稀疏性求解方程的方法。 剛度矩陣為什么會稀疏？因為在有限元中，一個節(jié)點的影響，只限于它所連接的那些單元。為什么？就是因為在前面，我們推導有限元時，在式（2）中，將熱傳導偏微分方程乘以的那個神奇函數。我們說過，是任意標量函數。既然是任意的，當然可以任意選取。然而我們沒有“任意”地、胡

15、亂地選取，而是居心叵測地，讓它恰恰等于有限元的插值函數！而這些插值函數，恰恰只在給定單元內部非零，在單元外邊一律為零！換句話說，插值函數只是些局部函數。我們讓等于插值函數，它也就具有了這種局部性。正是的這種局部性，使得一個節(jié)點的影響，只限于它所連接的單元。有限元方法，之所以能夠在計算力學領域里令人眼花繚亂的各種各樣的計算方法中，獨領風騷，傲視群雄，鶴立雞群，至今幾達50年之久，其全部奧妙，皆出于此！ 為進一步考察這些影響到底有多少，我們來看下面的例子。使用前面的MATLAB有限元程序，給 Cdiv 的值輸入 8, Tdiv 輸入 2 ，就會生成以下網格。它將圓周分成8份，厚度分成2份

16、。圖中括弧內是單元號碼，其余數字為節(jié)點號碼?？梢钥闯觯?節(jié)點只與第1和第8單元相連，其影響也就只限于這兩個單元。這里所說的影響，就是在剛度矩陣中，第1和第8單元的所有節(jié)點，即第1、2、5、4、22、23節(jié)點，要發(fā)生關系。也就是說，在總剛度矩陣（高斯消去法中的K矩陣）中，相應的行與列上的元素非零。例如在第1行當中，第1、2、5、4、22、23列的元素非零，其余元素為零。我們知道，總剛度矩陣的列數與行數是相等的，在本列情況下，列數等于24。在第1行當中，只有6個元素非零，其余18個元素都是零。同理，第4節(jié)點只與第1和第2單元相連，其影響也就只限于第1和第2單元。因而，在總剛度矩陣第4行當中，第1

17、、2、5、4、8、7列的元素非零，其余18個元素為零。第2節(jié)點影響4個單元，分別是第1、9、8、16單元，因而在總剛度矩陣第2行當中，非零元素最多，達到9個，即第1、2、3、4、5、6、22、23、24列的元素非零，其余15個元素為零。如此，可想而知，總剛度矩陣的每一行的大部分元素都是零。 現在我們要考慮怎樣利用這些零元素了。在以上網格中，共有16個單元，24個節(jié)點。使用高斯消去法，會生成24 x 24 的總剛度矩陣，即有24行，24列。而使用波前法，總剛度矩陣的行數雖然不變，也是 24， 但列數僅為11（至于為什么是11 ，下面要詳細討論）。最重要的，在本例情況下，是波前法根本就不

18、需要將 24 x 11 的總剛度矩陣存在內存中，而是存在硬盤上的。內存里，波前法只需要放一個 11 x 11 的所謂波前矩陣就行。那么，什么是波前矩陣呢？就是在某一時刻，彼此發(fā)生關系的節(jié)點的剛度系數組成的矩陣。這個矩陣是方的，其中含有最多非零元素的那一行就叫波前。什么叫某一時刻？就是某一單元！如前面MATLAB程序所示，計算有限元剛度矩陣有兩重循環(huán)，最外面那層循環(huán)，是對單元循環(huán)，即從編號為第一的單元，到編號最大的單元。在波前法中，當循環(huán)到某一單元時，在計算該單元剛度矩陣以后，還要看看哪一個節(jié)點可以消去了，也就是消元。被消元的節(jié)點，對以后其它單元剛度矩陣就不再有影響了，該節(jié)點的剛度系數就可以存入

19、硬盤指定文件中，而這些系數就可以從波前矩陣中清除掉，以便空出位置來，存儲其它節(jié)點信息。因此，一個節(jié)點可以被消元的時間，是可以用單元循環(huán)的進度來度量的。 那么，一個節(jié)點，什么時候可以消元了？就是與該節(jié)點相連的所有單元都循環(huán)到了的時候。例如，若循環(huán)從第1單元開始，當循環(huán)完了第2單元（計算單元矩陣并裝配到波前矩陣中）時，第4節(jié)點就可以消元了，因為第4節(jié)點所連接的2個單元都循環(huán)到了。同理，當循環(huán)完了第3單元時，第7節(jié)點也可以消元了。而第1節(jié)點的消元要等到循環(huán)到第8單元。而第2、3節(jié)點的消元時間最遲，要循環(huán)到最后一個單元，第16單元。 據此，可以編制一個節(jié)點消元時間表，也就是循環(huán)到了

20、哪個單元，該節(jié)點便可以被消元了。算法很簡單，就是查找每一個節(jié)點所連接的最大單元編號。程序如下：程序二：計算節(jié)點消元時間以上程序中，第三行的ICO變量是個兩維數組，18行，4列。它的每一行代表一個單元，該行的4列給出該單元的四個節(jié)點。這段程序執(zhí)行的結果，是在一維數組變量NodeDeactiveTime中定義每個節(jié)點可以消元的時間（即單元號）。此時，NodeDeactiveTime 的值是：8    16    16     2    10  

21、60; 10     3    11    11     4    12    12     5    13    13     6    14    14   

22、0; 7    15    15     8    16    16。第1個數“8”代表第1節(jié)點的消元時間是8（單元）；第2個數“16”代表第2節(jié)點的消元時間是16（單元）；余類推。請注意，消元最早的節(jié)點是第4節(jié)點，時間是“2”（單元）。其次是第7節(jié)點，時間是“3”（單元）。我們下面介紹在波前矩陣里如何消元時要用到這兩個信息。 知道了各節(jié)點的消元時間，就可以計算波前矩陣的最大階數了。程序如下：程序三：計算波前矩陣的最大階數在以

23、上程序中，第1行開了一維輔助數組，NodeInFront，用于記錄每一個節(jié)點是不是在波前矩陣中，“1”表示在，“0”表示不在。第2行將我們要計算的波前矩陣的最大階數maxFrontWidth的值初始為零。第3行開始對單元循環(huán)。對編號為 i 的單元，第4行從單元總表ICO中取出該單元的節(jié)點（本例為四個節(jié)點）；第5行將這些（四個）節(jié)點在NodeInFront中的值定為“1”，代表它們進入了波前矩陣。注意，此時，有的節(jié)點可能已經在波前矩陣中了，即它們在NodeInFront中的值已經是“1”。但這沒有關系，現在只是重新再植一次“1”，再一次表示該節(jié)點在波前矩陣中。第6行計算 maxFrontWidt

24、h。就是將NodeInFront中的“1”相加，再與當前maxFrontWidth比較，誰的值更大，maxFrontWidth就取誰的值。也就是說，maxFrontWidth的值只增不減。第8行對 i 單元的（四個）節(jié)點循環(huán)，查找其中每個節(jié)點是不是到了消元時間（由NodeDeactiveTime給出）。 如果是的，第10行的邏輯變量deactive的值為“1”，并在第12行將該節(jié)點在NodeInFront中的值改為零。這表示該節(jié)點被清除出波前矩陣。 這段程序結束全部循環(huán)后，便得到maxFrontWidth的值為11，就是波前矩陣的最大階數。前面說的波前法總剛度矩陣是24 行 x 11

25、列。其中的11，就是這樣得出的。它的含義，就是在整個計算過程中，某一時刻，同時存在于波前矩陣的節(jié)點數，其值最大為11。  以上程序，實際上模擬了節(jié)點進入和離開波前矩陣的裝配和消元過程。下表給出波前矩陣中的節(jié)點隨著單元循環(huán)的整個變化過程。第1列是單元號碼。第2列是將該單元的剛度矩陣裝入波前矩陣后，波前矩陣中的節(jié)點。第3列是這些節(jié)點的數目。第4列是此時刻可以消元的節(jié)點。第5列是消元后，將消元節(jié)點清除之后，波前矩陣里剩下的節(jié)點數，即第3列減去第4列?？梢钥闯觯瑑纱芜_到最大波前，分別當循環(huán)到第7單元和第10單元時。讀到這里，對照前面那個有限元網格，您要是能夠理解這個表中所有數字的來龍去脈，您

26、就理解了波前法！剩下的，是一些技術上的細節(jié)，可以通過閱讀下列全部程序來最后徹底“摳”懂。  這些技術細節(jié)里，比較難懂的，是如何將單元剛度矩陣裝配入波前矩陣。單元剛度矩陣是個4階矩陣，即4x4矩陣，代表了單元4個節(jié)點之間的相互影響。例如，第1行里的4列元素，是單元4個節(jié)點對單元第1節(jié)點的影響；第2行里的4列元素，是那4個節(jié)點對單元第2節(jié)點的影響；余類推。 每個單元節(jié)點的局部編號都是1、2、3、4，但它們的總節(jié)點編號是不同的，而且是唯一的。單元節(jié)點的局部編號與其總節(jié)點編號的對應關系，在本文的程序中，由二維數組ICO給出。在用普通高斯消去法解方程時，任給一個單元剛度矩陣，我們可以

27、通過ICO表，查到該單元4個節(jié)點的總編號。而節(jié)點總編號與總剛度矩陣 K 的行與列之間具有一一對應的關系。第1節(jié)點占據第1行第1列，第2節(jié)點占據第2行第2列，第n節(jié)點占據第n行第n列。所以，使用普通高斯消去法解方程時，把單元剛度矩陣裝配進總剛度矩陣的算法很簡單，如本文第一個程序的第122行所示，只一行程序便可實現。打個比方，這種情況下的總剛度矩陣 K，就像占地廣闊的大觀園，有許多房子，紅樓夢里的所有人等，例如金陵十二釵，各有各的住所，還有許多空地（相當于K里的零元素），可供黛玉葬花，姐妹們賞月吟詩。 波前法的裝配算法則要復雜得多。因為波前矩陣小，不能同時裝下所有節(jié)點的所有元素。這種情況，就如同旅

28、店。世界上所有旅店是住不下裝得下世界上所有居民的，旅客必須有進有出才行。節(jié)點進出波前矩陣，就如同旅客進出旅館。有的旅客住的時間長，有的短。節(jié)點在波前矩陣呆的時間也一樣，有長有短。隨著節(jié)點的進進出出，波前矩陣的每一行和列都可能先后被不同節(jié)點占用，就像旅館里的每個房間都會被不同旅客先后住過。所以，波前矩陣的行與列，與總節(jié)點編號之間，不再有像使用普通高斯消去法解方程時那樣的一一對應的固定關系。單元矩陣的某一行、某一列的元素，應該放到波前矩陣的哪一行、哪一列，是動態(tài)分配的。有兩種可能：如果某節(jié)點已經存在于波前矩陣中，那么就把該節(jié)點在單元矩陣中的元素加到波前矩陣的相應元素上（因而需要知道它在哪里）；反之

29、，如果某節(jié)點還沒進入過波前矩陣，那么就在波前矩陣給它分配一個自由的行與列。在下面的程序中，我們使用一維數組freeLines來記錄波前矩陣每一行（同時也是每一列。我們假定行數等于列數，也就是說，一個節(jié)點占據了第 n 行，它也就占據了第 n 列）的狀態(tài)：0 表示自由；大于 0 表示已占用（即占用該行的節(jié)點) 。我們還使用一維數組GlobalID2FrontId來記錄每一節(jié)點在波前矩陣占據的行數（=列數）。也就是說，給出一個節(jié)點的總編號，就能在GlobalID2FrontId中查到它在波前矩陣中的位置。這個位置如果是零，就表示該節(jié)點還沒進入過波前矩陣。這些算法的具體實施可見下列程序中的第

30、148-168行。程序四：使用波前法解線性方程組的二維熱傳導有限元程序 消元過程消元也是波前法里比較難于理解的技術細節(jié)。我在調試波前法程序過程中，所花時間最多的，就是消元。而現在，當我想解釋消元過程時，發(fā)現遇到了更大的困難。因為我不知道應該從哪里說起。消元是個細活兒，必得有研究針尖上能站立多少天使的牛角尖精神才能弄懂。不是所有人都有這種精神。而為寫此文，已經花了如此多工夫，幾乎超過了開發(fā)程序的時間！ 所以，還是算了吧。一簡對三繁，我就把消去第一個可消節(jié)點（第4節(jié)點）的過程列表于下。您若看得明白，就算徹底懂了波前法。 為方便說明問題，我將原問題的內壁邊界條件換為絕熱

31、，外壁換為溫度給定（程序四第63至64行須作相應修改）。這樣，第4節(jié)點就不經過邊界處理那一段程序（第183至184行），而是進入消元那一段（第186至190行）。 表1，將第1單元的單元矩陣裝配入波前矩陣后的波前矩陣然后計算第2單元的單元剛度矩陣，結果如下：注意，第2單元的單元剛度矩陣恰好與第1單元相等，這純屬巧合。巧合的原因是：1. 兩單元節(jié)點關于45°線對稱；2. 介質各向同性。  表3，將第2單元的單元矩陣裝配入波前矩陣后的波前矩陣  表4，第 4 節(jié)點消元后的波前矩陣紅字顯示因為消元而出現的非零元素。例如，第7、8兩節(jié)點，本來與第1節(jié)點不共單元，

32、對第1節(jié)點沒有影響，但現在有了。在消元中，我們用第1行減去第4行（乘除兩個系數后），因而第4行的元素全部出現在第1行中。 消元的結果，是第4行主元素所在的列，除主元素外，其余元素全被消成零。此時，就可以將第4行存入指定文件中，如程序第203行所示。程序第205行將此時的波前節(jié)點號寫入另一指定文件。以后回代需要。程序第204行寫方程Ax = B的等號右端自由項。 有關數據存入文件后，就可以將第4行清零了（程序第206至208行）。此時波前矩陣如下表所示。表5，第 4 節(jié)點清除后的波前矩陣波前法與高斯消去法的效率之比較 比較結果見下表。第1列是比較所用的6種網格。Cd

33、iv 代表沿圓周的分割，Tdiv代表沿厚度方向的分割。第2列是節(jié)點數。第3列是波前矩陣的最大寬度。第4列和第5列分別是使用程序四（波前法）和程序一（高斯消去法）解題的時鐘時間（elapsed time）。由上表可看出，當節(jié)點數很少（24和80）時，兩種方法算題的時間差不多，但當節(jié)點數達到288時，差別開始顯示出來，高斯消去法需時比波前法將近多一倍。64 x 16的網格（1088節(jié)點），時間差別達到15倍。96 x 24 （2400節(jié)電），差別到了60倍。到了最后一個網格，128 x 32（4224節(jié)點），運轉程序一（高斯消去法）需內存563MB，這已經超出了我的計算機內存（512MB），因而所需時間可認為無窮大。  小結：波前法的過去、現在和未來在人類計算技術發(fā)展史上