工科數學
——祖沖之思路

鐘萬勰

大連理工大學

工科數學
——祖沖之思路

鐘萬勰

大連理工大學

前 言

中國的大學工科數學教材，講的全部是國外數學家成就，“言必稱希臘”。作者常常懷疑，難道歷古以來光輝燦爛的中國文化，竟然在數學方面一無所成？然而，我們是在工程力學方面做研究與應用的，畢竟是數學外行。教育部邀請許多數學家所編寫的教材寫成為：全部由洋人所貢獻，也是無可如何。有疑問：中國工科數學教材，難道中國人的貢獻竟可被不屑一顧地完全忽略，連立錐之地也不容嗎？

對此難以認同，不平則鳴么！作者不是數學家，難免有所舛誤，有錯就請批評，不用客氣。

作者這些年在動力學方面努力，尤其是有約束的動力學計算分析，發現作為機器人動力學基礎的約束動力系統微分-代數方程(DAE, Differential-Algebraic Equation)的求解，例如現代著名著作

E.Hairer G.Wanner:Solving ordinary differential equations II-stiff and differential-algebraic problems 2nd ed.ch.7.[M],Springer,Berlin,1996. E.Hairer,Ch.Lubich and G.Wanner:Geometric-Preserving Algorithms for Ordinary Differential Equations.Springer,2006.

大力論述DAE的求解，效果不理想。一些著名洋人軟件，例如廣泛應用的ADAMS，其數值結果也不理想。因為這些著作的求解方法，是先進行微商，將約束方程歸化到微分方程，其微商次數稱為Index，可稱為Index法。看起來約束條件處處滿足，而實際上數值結果的約束條件滿足不行。按“實踐是檢驗真理的唯一標準”，不行。

鐘萬勰，中國科學院院士，大連理工大學工程力學研究所所長、教授。工程力學、計算力學專家。1956年同濟大學橋梁與隧道系畢業。1993年當選為中國科學院院士。

上世紀60年代發現潛艇耐壓錐、柱結合殼失穩的不利構造形式。70年代與小組基于群論研制了大量工程應用軟件，并主持研制了三維大型有限元系統JIGFEX／DDJ。80年代提出了基于序列二次規劃的結構優化算法及DDDU程序系統提出結構極限分析新的上、下限定理，繼而又提出了參變量變分原理及相應的參變量二次規劃算法用于彈－塑性變形及接觸問題，是中國計算力學發展的奠基人之一。1989年以來，發現了結構力學與最優控制相模擬據此又提出了彈性力學求解新體系與精細積分的方法論。

中國著名南北朝數學家祖沖之(429～500)在計算圓周率π時，已經達到π＝3.1415926～3.1415927之間，唐朝魏征等撰寫的史書《隋書》中有記載∶

（古之九數，圓周率三，圓徑率一，其術疏舛。自劉歆、張衡、劉徽、王蕃、皮延宗之徒，各設新率，未臻折衷。宋末，南徐州從事史祖沖之，更開密法，以圓徑一億為一丈，圓周盈數三丈一尺四寸一分五厘九毫二秒七忽，朒數三丈一尺四寸一分五厘九毫二秒六忽，正數在盈朒二限之間。密率，圓徑一百一十三，圓周三百五十五。約率，圓徑七，周二十二。又設開差冪，開差立，兼以正圓參之。指要精密，算氏之最者也。所著之書，名為《綴術》，學官莫能究其深奧，是故廢而不理。）

學官厲害，我不懂就廢而不理。當然這是古代官僚機制。找官告狀，就得先跪下。古代數學受到嚴重打擊。祖沖之的《綴術》為之失傳；然而今天畢竟還有劉徽（三國魏人）的《九章算術》傳世。今天我們應將祖沖之算法挖掘出來，“繼絕世”，融合現代數學，解決今天的問題，發揚光大。貫通古今，融合中西。

如所周知，π是現代數學不可回避的基礎，看到中國祖師爺的重要貢獻，欣喜不已。于是我們就探討祖師爺祖沖之在當年條件下是怎么計算的。由此引申下去，其實許多數學基本概念，例如無窮序列，極限等等，中國祖師爺也是明白的。并且在實際上應用了，按此路線探究下去就接地氣。尤其是接中國人的地氣。

所以，中國數學并非一無所成，而是挖掘不夠。竊以為，既然在中國大學講數學，就不應將中國數學祖師爺的工作忽略，而應有所傳承。我們不可“光說不練”，本文就從祖沖之如何計算圓周率π講起，加以傳承，貫通古今，融合中西。中國數學也應占有一席之地的。可先從中國工科數學教學做起。

2014年5月，鐘萬勰參加了在廣東工學院的工科數學教材會議，介紹了祖沖之類算法。希冀大學數學教材要加入中國元素。本文是以實際工作貫徹該意圖。不僅是講講而已，而要有實際行動的。期待官方能任領導改革之職。下面就由探討祖沖之的工作開始。

著作[10]的序言中講：“1999年5月，教育部委托上海大學在錢偉長教授主持下召開了一次應用力學教改的會議。該會議使我下決心寫出這本書，為此花費了大量的精力。”表明鐘萬勰很早就有志于教學。雖然曾努力推動，但思路不得法，效果不理想。現在再一次努力，希望能發揮應有作用。但鐘萬勰困惑：教育改革沒有SCI，得不到評價，又怎能鼓勵青年教師呢。官僚機制厲害呀！當官只要數數SCI就可以了，懶政！此“官不任事”者也。面對這種評價機制，只能長嘆：“此天之亡我，非戰之罪也！”，對此唏噓不已。

一，祖沖之是如何計算圓周率π的估計

中國人古代早就關注圓周率了，所謂“周三徑一”，那只是一個約略的估計。設平面上有一個直徑為2的圓，半徑r＝1，要計算圓周長度，當然是2π了。但當年祖沖之又是如何計算的呢？

中國在東周時期，已經發明了算盤，并且申了遺。算盤是古代的計算機。并且中國勾、股、弦的定理也早已發現，周公時期(公元前1000多年)已經有所記載，稱商高定理。古代希臘有畢達哥拉斯者，對此有系統論述，所以大家說這是畢達哥拉斯定理。我們可仍稱呼為商高定理。因為祖沖之未必知道有畢達哥拉斯者，他實際可使用中國傳承下來的商高定理。

平面上兩點之間的最短距離是其連接直線的長度，今天說法是平面歐幾里得幾何的短程線。這些基本概念是祖沖之算法的基礎。

一個圓有內接正多角形，有如切西瓜(古代稱為割圓術)，1分為2,2分為4,4分為8,…。每次劃分全部是n=2i，i＝0,1,2,3,4,…。圓周生成了n條邊的內接正多角形。當切了i=2次，就有內接4邊形，總是半徑r=1，此時每塊等腰三角形的張角是90o。切下的內接正多角形圓弧，不計算圓弧長度而用短程線的長度替代，成為等腰三角形。i=2時兩點連直線的長度是l2=。等腰三角形有中垂線將三角形劃分為2個直角三角形。其勾的長度是a=l2/2，弦的長度就是r=1。中垂線的長度（股b）可根據商高定理a2+b2=r2=1而計算。股的長度b當然小于半徑r=1。于是得到1-b，這是將中垂線延伸到圓周點p的長度，見圖1。

圖1.割圓術

以延伸長度1-b為勾，以前面的a=l2/2為股，再根據商高定理可計算其弦，就是n=2i+1正內接多邊形的邊長li=3。于是計算圓周長度2π的近似，內接正多角形邊的總長是2ixl。i

圖1為半徑為r=1的圓，計算圓周率，采用割圓術計算，用內接正多邊形的總周長近似為圓的周長。實際計算的可以是半個圓周。假設正n邊形的總邊長為li，n=2i，用li作為弧長的近似。按上文所述，勾a=li/2，弦為1，而按商高定理得b=[1-(li/2)2]1/2的股。

再用商高定理。以伸長線1-b為勾

1-b=1-[1-(li/2)2]1/2=1-[4-li2]1/2/2=[2-(4-li2)1/2]2，以a=li/2為股，則給出2li+1=2×[(1-b)2+(li/2)2]1/2=[4(1-b)2+li2]1/2。將4(1-b)2=[2-(4-li2)1/2]2=8-4(4-li2)1/2-li2，代入有

這樣，就形成了一次等腰三角形的分割，式(1.1)形成了遞推回歸形式。無非是將2條邊的長度si+1=2li+1相加在一起而已。

具體計算可從半個圓開始，取為i=1,s1=l1=2。即計算半個圓周的弧長，半徑1，張開角是1800。結果弧長應該是π，而與i＝1，s=l1＝2相差很多。再分解一次則有l2，而張開角成為900的2段圓弧，近似為2個三角形其邊長之和是s2=2l2=2=2.8284；與π逼近了些但相差依然很大。再繼續劃分2個等腰三角形，900的一段圓弧，成為張開角各為450的2段圓弧，近似長度s3=2l3=22l2；如此繼續下去…。從li計算li+1的公式是(1.1)。其數值結果羅列于表1之中。

內接正多角形的總邊長，永遠趕不上2π，因為連接線是直線短程線；然而可不斷逼近2π，得到序列。所以，得到的無窮序列是無限逼近于2π的，2π就是其極限。具體數值計算給出

表1：不同正多邊形對應的圓周率

表1為用SiPESC集成平臺調用的Python編程計算得到的表格，可見當取正16384邊形時，近似計算得到圓周率為.14159263，此過程不斷進行，得到了無窮序列，其極限就是圓周率π。

所以說，祖沖之肯定有無窮序列和極限的概念。通過割圓術這樣的具體問題來講述無窮序列和極限，接地氣，容易理解。教學么，本來應當如此。

中國大學工科數學教材可在祖沖之的基礎上講授下去。中國特色、中國元素么。直觀，而又接地氣，畢竟工科與理科純數學的要求不同。形成切合中國自己的那一套，讓學生容易懂，那才是今天應關注的。

科學網2015-5-24傳達韓啟德說：“中國科學家缺位科普”。在工科數學中，中國數學家也是缺位的。讓我們來改變這個不合理的現狀吧。

感慨：人云：“如欲亡其國，應先滅其史”！滅其史干什么，為了喪其志！連爸、媽是誰都不知道了，喪了志，渾渾噩噩，任人擺布，那就什么也沒有了。警惕呀。

中國數學如自動缺位，前途堪憂。

應從基本概念方面提煉出祖沖之算法的特點：圓弧全部用短程線的直線代替。除節點外，短程線的點全部不在圓周上。將圓周看成約束，則除節點外，只要短程線的要求得到滿足，不必再考慮其約束。這是在二維歐幾里得幾何條件下得到的。我們當然希望將祖沖之思路推廣到更廣泛領域中去。事實上這是成功的。在動力學DAE的推廣，比國外許多著作用的Index法好多了。見后文。

雖然祖沖之當年未必知道有短程線之說，但直觀上一定知道直線是最短距離。直觀沒有什么不好，不要苛求。馮·諾伊曼的文章《數學家》[1]應該讀一讀，不要追求絕對的嚴格性。

數學大師馮?諾伊曼指出：“能以更好的形式推廣為更有效的新理論的理論將戰勝另一理論。…必須強調的是，這并不是一個接受正確理論、拋棄錯誤理論的問題。而是一個是否接受為了正確的推廣而表現出更大的形式適應性的問題”，見文獻[1]。

上文推測了祖沖之計算圓周率的方法。如僅僅只是用在二維歐幾里得空間，那未免太局限了。祖沖之算法的思路非常寶貴，應尋求將其使用于廣大領域。有中國特色成果，應予以收復失地。

對于中國大學工科數學教育來說，本文不免有些“說三道四”。但我在大學工科工作了一輩子，難道連這些話也不該說嗎！

黨中央推動改革。官員有公權力，希望能有擔當，領導我們改掉不合理的教材吧。希望有關部門能予以關注，改革3的路就在腳下，難道工科數學就不用“認祖歸宗”嗎？ 振作起來干吧！

二，無約束動力學的時間積分

牛頓的奇跡年是1665～1666，后來1687年牛頓出版了巨著《自然哲學的數學原理》是用希臘文寫就的，簡稱principia，是近代科學開始的標記，同時也創立了微積分，這是劃時代的成就。書名表明是自然哲學，經典動力學的數學理論由此開始。牛頓主要考慮天體力學方面的問題。隨后分析動力學則得到當年眾多大數學家的關注，Bernoulli,Eule r,Lagrange,Hamilton,Jacobi,Poisson,Poincare, …，經典動力學就是他們奠基奉獻的。Euler-Lagrange方程，Lagrange函數，Hamilton變分原理，正則方程，Hamilton-Jacobi偏微分方程，正則變換，作用量，攝動法，…，成為優美的數學分析體系。牛頓提出微分方程求解的方法；Euler-Lagrange提出能量方法，Hamilton則提出了正則方程，其自變量擴展到了2n，其中n是位移自由度數。位移法的Euler-Lagrange方程，是n個二階微分方程；而Hamilton正則方程則是狀態空間的2n個一階對偶微分方程。

經典力學眾多著作是許多著名數學家努力耕耘的成果。內容豐富、數學理論優美，難以窮盡。這里不過是取其中心部分內容講解，畢竟是教數學么。然而對于Hamilton體系的認識，卻有曲折過程。

馮康說[4]：從歷史上考察人們對Hamilton體系的評價，Hamilton本人是從幾何光學著手創建他的理論模式的。1834年Hamilton曾說：“這套思想與方法業已應用到光學與力學，看來還有其他方面的應用，通過數學家的努力還將發展成為一門獨立的學問”，這僅僅是他的期望。19世紀同時代人對其反應則很冷淡，認為這套理論“漂亮而無用”。著名數學家F.Klein(International Congress of Mathematicians,ICM，1897年第一屆蘇黎士大會主席)在對Hamilton體系的理論給予高度評價的同時，對其實用價值亦持懷疑態度，他說“這套理論對于物理學家是難望有用的，而對工程師則根本無用”。這種懷疑，至少就物理學的范疇而言，是被隨后的歷史所完全否定了。到了20世紀量子力學的創始人之一Schroedinger曾說：“Hamilton原理已經成為現代物理的基石…，如果您要用現代理論解決任何物理問題，首先得把它表示為Hamilton形式”。說明了Hamilton狀態空間法的重要性。

量子力學當年面臨的一個關鍵問題是光譜分析。觀察到的譜線分裂現象，被Hilbert的大弟子，大數學家Hermann Weyl用群論分析清楚[2]。進入Hamilton狀態空間后，H.Weyl在1939年又提出Symplectic對稱[3]。華羅庚先生音譯為‘辛’。H.Weyl提出了Hamilton正則方程的辛對稱，表明了動力系統運行的基本性質，表征了動力學的辛階段。

辛，是大數學家Hermann Weyl在名著[3]中提出的。作者寫道：“The name ‘complex group’ formerly advocated by me in allusion to line complexes, …has become more and more embarrassing through collision with the word ‘complex’ in the connotation of complex number.I therefore propose to replace it by the Greek adjective ‘symplectic’ ”

表達了為避免“complex”容易產生的混淆，特地引入的希臘形容詞‘symplectic’。

H.Weyl的工作吸引了純數學家們高度關注。從微分幾何的角度出發，給出了“辛幾何”的定義；但其理論過分高深以致讓工程師難以接受。今天已經是信息時代，計算機的發展使人們認識到“計算科學”的重要性。“計算科學與理論、實驗共同構成現代科學的3大支柱”的論點得到了廣泛認同。計算科學是當代數學的大發展。計算科學的主要特點是離散。計算機本身就是離散的么。

“世界潮流，浩浩湯湯，順之者昌，逆之者亡。”

辛幾何：數學家從純數學微分幾何出發，定義：微分形式，切叢、余切叢(tangent,cotangent bundle)，交叉外乘積(exterior product)，Cartan幾何等。艱澀難懂。離散后就成為：

辛代數：容易理解、掌握。回歸到力學，離散系統。可從虎克定律講清楚，見[9]。自然科學基金1993年支持了課題(19372011)“結構動力學及其辛代數方法”。辛代數的提法不是現在才有的。

我國數學家馮康率先提出[4]，動力學積分的差分格式應保辛，差分格式就是離散積分。國外著作[5]也接受了。

祖沖之算法首先是用于離散系統的，切實符合計算科學的需要，所以前程遠大。自然可用于離散后的動力學方程。牛頓是動力學和微積分的祖師爺(德國的Leibnitz同年代也從幾何方面提出了微積分)；同年代的虎克(Robert Hooke)是結構力學的祖師爺，彈性力學不能回避虎克定律。然而當年兩人不和[6]，因此分道揚鑣，各自獨立發展。Hamilton狀態空間方法未能進入彈性力學領域，直至著作[7～8]的出版，才改變了求解思路。

作者隨后在著作[9～12]，建立了分析動力學與分析結構力學的模擬理論。會同以前建立的計算結構力學與最優控制間的模擬理論[13,14]，將多個領域皆奠基于狀態空間的數學理論。于是各個領域的數學計算方法可以互通有無了。基礎就是狀態空間法的辛對稱。

1900年在第二次世界數學大會上，D.Hilbert做了一個著名報告《數學問題》，深刻影響了20世紀數學的發展，指出：“只要一門科學分支能提出大量的問題，它就充滿著生命力；而問題缺乏則預示著這門學科獨立發展的衰亡或終止。”能看到問題就好，表明可繼續發展。辛指出Hamilton正則方程有對稱性，是經典力學繼續發展的機會，有大量的新問題。《數學問題》中提出23個問題，其中第23號是變分法的進一步發展。Hilbert說：“我已經廣泛地涉及了盡可能是確定的和特殊的問題…，用一個一般的問題來做結束…我指的是變分法”，見文獻[15]。變分法不單純是一個數學問題，而是一個方向，是大師的遠見卓識。

變分原理的提出，由來已久。“大自然總是走最容易和最可能的途徑”，這是費爾馬(Fermat)著名的自然哲學原理。1744年，J.Bernoulli提出了“最速下降線”問題，大體上可認為是數學變分法的開始。以后蓬勃發展，Euler-Lagrange方程，繼而總結為Hamilton變分原理。分析動力學與相應的常微分方程理論的成功，自然要發展到偏微分方程。到位勢理論，有Laplace方程的求解，Green,Gauss,Dirichlet等先后指出，可將其轉化為變分原理。Riemann(1826～1866)將其命名為Dirichlet原理，屬于橢圓型偏微分方程理論，本來在蓬勃發展。然而在1870年發生了曲折，強調數學嚴格性的維爾斯特拉斯(Weierstrass)否定了Dirichlet原理；然而數學物理中許多重要結果都依賴于此原理而建立。1899年,Hilbert用邊界條件的光滑化保證了極小化函數的存在，從而挽救了Dirichlet原理。表明變分原理經歷了“鳳凰涅槃，浴火重生”的過程。

此后一個多世紀以來，變分原理有了巨大進展，從有限元法發展出來的“計算科學”也是變分法的發展；而辛也屬于變分法的發展。我國在變分原理的研究方面是有成就的，胡海昌的三類變量變分原理蜚聲世界。有限元法離散取得了巨大成功，有限元單元推導的基礎就是變分原理。動力學變分原理稱為最小作用量原理。做一點解釋如下。

設物體有n個自由度。Lagrange函數L(q, )=T-U是動能T減勢能U，是n維位移q(t)和速度(t)向量的純量函數。在數值求解之前，介紹一下國外發展的許多方法是有益的。

隨著計算力學、有限元，以及航空航天等的推進，動力系統的求解成為不可缺少。絕大部分課題分析求解無望，只能數值求解。國外數學家大力發展Runge-Kutta，Newmark法等許多差分算法，可謂五花八門；繼而國外進行了許多研究，一批保辛差分算法相繼出現[5]。動力學系統的積分特別講究所謂首次積分(fi rst integral)，其實就是積分不變量。其中能量守恒特別重要。

然而，繼續的研究卻走入了歧途。1988年，[17]提出了“保辛則能量不能守恒”的誤判。文獻[18]用參變量達成了離散近似系統在保辛的同時依然保守的算法，并在著作[19]中給出了系統的論述。表明計算科學的發展，催生了辛對稱的深入發展，并且還應當突破矩陣群論的等維數限制，還有時間滯后等問題，前途廣闊。著作[19]只是初步試探而已。

既然分析求解困難，則只能離散數值求解，差分法是常用的。我國數學家馮康提出，差分格式應保辛，棋高一著。設時間劃分為步長η的一系列節點 to=0,…,tk=kη,…。逐步積分要在積分出tk-1的狀態后，進一步積分時間區段k#:[tk-1,tk]，以計算tk的狀態。此時可用最小作用量原理。時間區段k#:[tk-1,tk]的作用量是

作用量可從例如[16]找到。動力系統通常是沒有約束條件的。

采用時間有限元法[21]積分。時間tk-1的狀態qk-1,pk-1已經計算出來，以下要求解時間tk的狀態qk,pk。在k#:[tk-1,tk]內位移用兩端位移向量qk-1,qk插值，例如簡單些線性插值。時間有限元法完成積分得到有限元近似的作用量

然后，按一般理論，有

其中qk-1,pk-1為已知。至于子矩陣Kaa,k就簡單地寫成Kaa了。于是

該表達式的右側全部是已知的qk-1,pk-1，用狀態向量表達則是vk-1；而右側是待求的qk,pk，用狀態向量表達則是vk。組成傳遞矩陣

其中子矩陣Kaa,k簡單地寫成Kaa了，等。其實這就是按[9]所述，2nx2n對稱矩陣Kk變換到傳遞辛矩陣Sk的公式，數值積分可按此進行。的公式可自行驗證。然而，此種傳遞辛矩陣的積分，卻不能保證每步的能量保守。如果要能量保守，則還應按文獻[23]加參變量的處理。

一般的動力學積分是沒有約束的，相對比較簡單些。而祖沖之類算法并不要求區段k#:[tk-1,tk]內滿足約束，而只要求在離散時間點處滿足約束條件。

雖然依然有許多問題，例如剛性(stiff)等不同尺度的問題。但畢竟沒有約束，相對比較好對付些。如果還有約束條件，同時要進行時間積分，就會出現微分-代數方程了。求解DAE是機器人動力學必要的基礎，祖沖之類算法可以發揮基本作用，下一節會講的。首先要積分無約束的動力學問題。

這方面有許多數值例題，不再重復。本節是初步的分析動力學的內容，時間有限元么。對于祖沖之類算法，不管約束條件而取“動力學意義下的短程線”來說，是不可缺少的基礎。有了基礎，就可講動力學的祖沖之類算法了。

三，祖沖之類算法

現代化推進使機器人成為發展的熱點。機器人要運動，當然需要動力學分析。機器人起碼是機構，機械工程離不開機構。機器各個部件皆可運動，但相互間是有位移幾何約束的，稱為完整約束。因此我們面臨受約束的動力學系統分析問題。國外發展了稱為微分-代數方程(DAE,Differential-Algebraic Equation)的求解方法。DAE的求解成為一個重要問題，祖沖之算法可推廣用于此類重要問題。著作[20]的副標題就直接寫上了stiff and differential-algebraic problems。論述了DAE的求解，但效果不理想。因為這些著作的求解方法，是先進行微商，將約束方程歸化到微分方程，其微商次數稱為Index，可稱為Index法。看起來約束條件處處滿足，而實際上數值結果的約束條件滿足不行。

上節，介紹了的是傳統的分析力學，無約束。哈密頓體系理論是很一般的，并不限于線性體系，而是對于保守體系的[16]。動力系統的數值積分吸引了眾多研究，有代表性的著作見[20,5]。許多實際課題引導到有約束的Hamilton系統(Constrained Hamiltonian system)，總是推導到微分-代數方程(DAE)。國外通常求解方法采用對于代數約束進行微商，Index法，并歸化到聯立常微分方程組再進行數值求解的方法。但數學理論方便，并不代表數值方法有效。從數值求解的角度看，反而使問題變復雜了。

代數方程本來就是微分方程組的積分，反而將它化成微分方程再進行數值積分，多余的“虛功”么，坐回頭車了。歸化到微分方程組再進行離散數值積分，前面的歸化步是連續的精確的，而離散數值積分是近似的，往返不等價么。方法論已經不合適了，怎么可能比原來的代數方程約束更精確、更方便呢。走偏路了！

采用祖沖之類算法的思路，時間積分的格點處，要求位移約束嚴格滿足；而在時間區段內則不管位移約束了，要求“動力學意義下的短程線”即可。現在是動力學，要求的短程線應是2n維狀態空間的短程線，因為在時間積分的區段內不管位移約束條件了。

格點處的位移代數方程顯然比微分方程容易數值處理，可在數值求解時作為等式約束同時迭代滿足。等式約束比不等式約束的處理容易多了。求解DAE的這方面可見文[24]，從其中的簡單例題，可看到約束條件的滿足是非常好的，表明index積分方法不可取。雖然這是國外眾多著作采用的，但數值例題表明不夠理想。

完整約束是直接對位移的，在每個時間節點可以方便地表達而不涉及位移的時間微分或動量的。容易在每個節點單獨處理的。理論上講，約束條件是完全在位移空間(Configuration)的。節點位移滿足約束條件，意味著位移和動量都只有n-nc個獨立分量，其中n,nc分別是位移和約束的維數。代數約束方程(位移空間的約束)不包含動量，故離散后的約束只與本節點的位移有關，而與相鄰節點無關。節點k的相鄰時間區段是k#:(k-1,k)以及(k+1)#:(k,k+1)。k#區段積分時得到的，受到的約束與(k+1)#區段的k點同，因為總是與k點的位移約束方程的切面相垂直(理想約束)。因此可直接用于。即動量在節點處是連續的。這給DAE的祖沖之類算法求解帶來了方便。

設q(t)是n維廣義位移向量，無約束時動力系統的Lagrange函數為，T,U分別為動能，勢能。但當運動時有nc維的理想約束，約束方程為

運動只能在約束下的超曲面(流形,Manifold)上運動，即q(t)已不是獨立的廣義位移。其動力方程的推導是引入nc維的Lagrange乘子函數λ(t)，組成擴展的Lagrange函數

導出的方程組是

對應的DAE為(這些是國外的方法論)

H(q,p)=E是保守的。此即在狀態空間q,p下的微分-代數方程[5,20]。國外求解方法論的基礎是微分方程，要推出聯立微分方程引起指標(Index)問題等，給理論與計算上帶來了復雜的附加因素。

微分—代數DAE方程的求解有廣泛應用，但非線性系統的DAE方程并不是可輕易地求解的，多種常見的差分近似已經有了較深入的探討。但差分的插值本來不滿足約束條件的。反而難以保證軌道在格點處滿足位移約束條件。只能在每步積分后再采用投影等修補手段，這種修補成了一種干擾。因其方法論是先用差分法離散約束產生的微分方程，然后再考慮約束。差分近似本身的精度就不夠，哪怕注意到了差分近似的保辛，但其投影修補是否保辛也是問題。隨著該思路有許多研究，但方法論不行。走偏路了。

“行成于思，毀于隨”。既然現有求解方法論不理想，就應重新探討。孫子兵法曰：“出其所不趨，趨其所不意”，本文要改變其求解方法論。其實中國數學祖師爺早已有世界首創的工作了。這就是引申為祖沖之類的算法。看看我們中國人是怎么求解的。

在計算圓周率π時，祖沖之給出的思路是，約束條件不必處處滿足，只要在節點處嚴格滿足就可以了，而相鄰節點間則可用短程線代替，而不用管nc維的約束條件。

同樣的思路也可運用于DAE的方程求解。先進行離散，約束條件只要求在節點處嚴格滿足，而軌道則可按無約束動力系統積分，(動力學意義下的短程線)，例如可用時間有限元近似積分。只要節點劃分足夠密，則軌道定會逼近真實解的。這樣的思路下導出的算法，可稱之為祖沖之方法論、祖沖之類算法。我們理應接住1500多年前祖師爺傳來的球，挖掘繼承，發揚光大，與近代數學融合。“古為今用，洋為中用”，向前沖之。

回顧Lagrange力學的基本思想，是先引入廣義位移以滿足全部約束條件的，但祖沖之類算法改變了其思路。Lagrange體系要處處滿足約束條件的廣義位移難以找到，但本文運用分析結構力學的基本概念，進入離散系統來分析。首先保證離散積分點的nc維約束條件嚴格滿足；再在時間區段內，用時間有限元[23]離散代替差分離散，然后運用作用量的變分原理以代替微分方程，離散的時間步內就不再考慮約束了，因插值函數本來不能滿足。離散在先，分析結構力學理論保證可達到每步積分的自動保辛，繼承了祖沖之方法論的思路，稱祖沖之類算法，特色思路么。祖沖之的古代，沒有微積分、動力學等，但其思路可以繼承(貫通古今)。融合近代科學理論(融合中西)，得到的算法可稱之為祖沖之類的算法。中國在現代數學發展中，不可缺位。祖沖之類算法，可提供一個例證。

基于分析結構力學理論的方法，可從最簡單的例題著手闡明。

例題1：質量m的單擺振動，取q(t)={x,y}T為未知數，而不用θ(t)。約束條件g(q)r2-(x2+y2)=0,r=1。初始條件x(0)=0.99,。

解：y坐標以向下為正，勢能U(x,y)=-mgry, gr=10m/sec2；動能，Lagrange函數，其變分原理

其中gr是重力加速度。

首先將時間劃分為步η的一系列節點to=0,…,tk=kη。按祖沖之方法論，只考慮時間格點處的約束條件，是離散的。對應地，其nc維約束的Lagrange參數向量λ也是離散的，不再有對nc維向量λ微商之說。對應于離散格點有λk,k=1,2,…的nc維向量待求。

擴展Lagrange函數的作用量，也要相應地離散。區段k#:(tk-1,tk)內按祖沖之思路，不用考慮約束，因此位移仍是n維。但到k號時間節點tk時，要將nc維向量λk考慮進去。

分別對區段位移分量進行簡單的線性插值，單元內部的約束條件則放松掉。運用有限元法近似計算作用量。線性插值有。但還要加入時間節點tk處的約束對偶向量λk，于是與完全無約束的(2.2)對比，應修改為

它相當于結構力學的區段變形能。整體變形能無非是全部單元變形能之和

由于插值公式，相當于作用量(區段變形能)用的有限元近似與真實作用量略有不同。Lagrange原理的位移約束條件已經在節點處嚴格滿足，區段內部的約束條件則由有限元插值自然就近似滿足了。其逐步積分求解是積分區段k#:(tk-1,tk)時，將qk,λk同時求解，然后再積分下一個時間步。

將節點約束條件(1.1)推遲到數值積分時一起處理。利用節點的nc維Lagrange參變向量λk。參變量λk則用于處理tk處的約束條件。單元內位移則用有限元插值

積分就得區段作用量sk(qk-1,qk,λk)。因未曾滿足節點約束條件，故位移qk仍是原來n維的，不是約束后的獨立位移。該作用量對參變量λk是線性的，可表達為

根據

與無約束系統相比，多了線性參變量 。組成各站的2n維狀態向量

按分析結構力學方法進行，從vk-i可遞推vk，該變換保辛，但帶有參變量λk。確定λk要根據節點約束條件g(qk)=0。

建議運用線性的有限元插值函數 N(t)，因區段內的約束條件并未嚴格滿足，用線性函數插值近似來滿足，可能效果好些。

本方法用參變量線性λk滿足節點約束條件g(qk)=0。未知數包含線性參變量λk與全部節點狀態vk，有2n+nc個未知數要求解，仍是非線性聯立代數方程，保辛。能提供的方程是帶有參變量λk的傳遞矩陣2n個方程，以及節點位移約束條件g(qk)=0，也是共2n+nc個非線性方程。

可歸納分析結構力學(帶參變量λk)的算法如下：

1)形成Lagrange函數(動能—勢能)，并引入約束的Lagrange參數，形成擴展Lagrange函數。

2)將時間坐標離散，得一系列的時間點to=0,t1,…,tk…，以各點的n維位移qk當作未知數。

3)按分析結構力學，計算區段(tk-1,tk)的作用量sk(qk-1,qk,λk)，可用有限元線性插值得作用量。

4)生成對偶向量

組成狀態向量。

5)根據

以及初始條件，對1號單元，用插值公式計算初始p0并組成初始狀態。

6)各單元的辛矩陣可從方程

會同約束條件g(qk)=0解出。成為根據qk-1,pk-1計qk,pk,λk，完成逐步積分的保辛遞推。

分析結構力學的算法運用時間有限元，不具體考慮約束對偶λ的微分方程，而是將λk逐點當作待定的常參數。因不涉及其微商，故不必考慮其微分方程，也不會產生DAE微分方程理論的Index問題。通過數值例題可看到祖沖之方法論，祖沖之類算法特色思路的效果。本課題的數值結果展示了優越性。

可提供更多的數值例題以供比較。本文展示約束保守體系的分析結構力學DAE有限元保辛算法的解。

例題2：選擇3維球面擺[5]p.210的例題。

解：按[5]取重力加速度取gr=1，擺的長度為1，質點質量為1，z坐標以向下為正。Hamilton函數為

初始位置為

初始動量為

積分步長分別取為0.03秒和0.1秒。圖2給出了球面擺質點的軌跡圖。因滿足約束條件是本算法的基本點，故不再展示約束。圖3給出了系統的Hamilton函數隨時間變化情況。文獻[5]213頁分別用辛歐拉算法和投影辛歐拉算法，給出了同樣問題的Hamilton函數和質點偏離約束面的情況，可以看到兩種算法的Hamilton函數都隨時間線性增長，保辛效果不理想。將本文的圖3與它們比較，可以看到本文算法的保辛效果好，對約束的滿足具有非常高的精度，具有更大的優勢。

例題3：空間雙擺問題：質點1的初始位置是

初始速度是{0.1 0 0}T，質點2的初始位置是

圖2.球面擺質點的軌跡圖(步長分別為0.03秒和0.1秒)

圖3，球面擺Hamilton函數隨時間變化，真實H=0.99680

圖4.質點1的軌跡，約束條件滿足非常好

初始速度是{0.2 0 0}T，積分步長為0.01秒。積分結果如以下。圖4是質點1的軌跡圖，圖5是質點2相對于質點1的軌跡，圖6是質點2的絕對軌跡。圖7是系統的Hamilton函數隨時間變化情況，可以看到Hamilton函數在兩個確定數值之間震蕩，不會線性地偏離，并且這兩個數值和系統真實的Hamilton函數相差很小，這說明保辛效果很好。

通過這些數值例題，讀者可看到如何計算求解，并看到其效果。這些數值結果，其位移曲線出現混沌現象。這是非線性系統的特性。該滿足的約束條件，滿足得非常好；等步長積分，能量Hamilton函數，雖然有所偏離但也滿足得很好了。比國外著名算法好多了。

DAE不但在非線性動力學求解中非常有用，(以上講的例題全部是經典動力學的)；而且在電網的網絡控制問題等課題中也很有用。雖然問題早就存在，好在DAE求解數值問題得到關注的時間還不是很長。以上提出的迭代法逐步積分，與洋人的Index方法完全不同，其效果從這些簡單例題中已經表達。用事實說話。

我們學習了許多著名數學家的名言，有些體會。教材就應當易學易懂，不怕人家說土，說水平不高。我們就是土生土長的么。Hilbert在《數學問題》中說：“清楚的、易于理解的問題吸引著人們的興趣，而復雜的問題卻使我們望而卻步”，又說：“在討論數學問題時，我們相信特殊化比一般化起著更為重要的作用。可能在大多數場合，我們尋找一個問題的答案而未能成功的原因，是在于這樣的事實，即有一些比手頭的問題更簡單、更容易的問題沒有完全解決或是完全沒有解決。這時，一切都有賴于找出這些比較容易的問題并使用盡可能完善的方法和能夠推廣的概念來解決它們。這種方法是克服數學困難的最重要的杠桿之一。”

圖5.質點2相對于質點1的軌跡，約束條件滿足依然非常好

圖6.質點2的軌跡，混沌出現了

圖7.Hamilton函數隨時間變化，真實H=28.85251

D. Hilbert《數學問題》：“數學中每一步真正的進展，都與更有力的工具和更簡單的方法的發現密切聯系著，這些工具和方法同時會有助于理解已有的理論并把陳舊的、復雜的東西拋到一邊。數學科學發展的這種特點是根深蒂固的。因此，對于個別的數學工作者來說，只要掌握了這些有力的工具和簡單的方法，他就有可能在數學的各個分支中比其他科學更容易地找到前進的道路”。

本文涉及了辛矩陣對稱群，讀者可能對群論有些怕。其實本文只用到群的最簡單原理。英國M.Atiyah(1990年當選的皇家學會會長，著名純數學家)說[26]：“群在自然中產生，它們是使事物運動的東西，它們是變換或置換…。理解這些東西的本性，并且使用它們才是目的”。又說：“重要的東西常常不是技術上最困難的即最難證明的東西，而常常是較為初等的部分。因為這些部分與其它領域、分支的相互作用最廣泛，即影響面最大”。“在群論中有許多極端重要的，并且在數學的各個角落到處都出現的東西。這些是較為初等的東西：群及其同態、表示的基本觀點，一般的性質，一般的方法—這些才是真正重要的”。作者也由此得益，這些觀點對讀者也許有益。

作為機器人基礎的DAE，基于祖沖之方法論得到的解，比國外著名算法的解好多了。中國祖師爺的優秀成果應努力予以挖掘繼承，融合近代數學，發揚光大。中華文化博大精深，不只是懷念懷念，而應以實實在在的東西體現出來。當然還應經歷挖掘、品味、提煉、繼承、融合近代數學、然后才能發揚光大。學習大學微積分，讀完后沒見到中國人實實在在的貢獻，遺憾；到現在信息時代，挖掘出祖沖之方法論，理應占有一席之地。我們應為此而站起來大喊大叫。“南海諸島一千多年前就在中國的管轄之下”，我們應理直氣壯地講，因為有事實依據。有人不承認，難道我們就不敢講了嗎！今天看到祖沖之方法論的優越性，可達到貫通古今，融合中西的境界。中國人就應占有這一席之地，有什么不敢講的。難道這也要等著，讓洋人來發揚光大嗎?! 他們會嗎？

當前中國一些人治學，重洋輕華，“言必稱希臘”。保辛離散體系，是新事物，洋人也措手不及，例如他們提出“不可積系統，保辛近似算法不能使能量守恒”的誤判，我們也予以了正名，中國人切不可自卑。所謂的SCI評價體系，同樣的論文，英文發表就比中文發表值錢，真荒唐！民族虛無主義，可恥！

對自己沒有信心，弱者的心態么，要不得。

當年紅軍英勇奮斗，請了李德，就打敗仗。道路決定命運么。評價體系是學術研究的指揮棒。提出一個什么SCI等的洋指標，完全放棄了自主，又請來了一個“李德”，特別不爽。

中國數學界泰斗、中科院院士吳文俊認為：“數學應該是有助于解決實際問題的，這應該是研究數學的主要目的，數學不是什么抽象的理論。…中國的應用數學應該以解決實際問題為出發點的，絕不能只做抽象的理論研究。”

改革！你講你的，我講我的。力求返璞歸真，是本文的特點。走自己的路，一定要自信。要有道路自信、理論自信、體系自信，要敢于擔當，真正做到‘千磨萬擊還堅韌，任爾東西南北風’。特色思路么。真正做到擔當，不容易呀！

[1]馮?諾依曼：數學在科學和社會中的應用。大連理工大學出版社，2009.

[2]H.Weyl∶The theory of groups and quantum mechanics, Dover.1928.

[3]H.Weyl∶The classical groups;Their Invariants and representations, University press, Princeton, 1939.

[4]馮康，秦孟兆：Hamilton體系的辛計算格式。浙江科技出版社，2004.

[5]E.Hairer,Ch.Lubich and G.Wanner∶Geometric-Preserving Algorithms for Ordinary Differential Equations. Springer,2006.

[6]杰克曼(J.Jakeman)：牛頓：上帝、科學、煉金術，大連理工大學出版社，2008。原版 Newton－A beginner’s guide.

[7]鐘萬勰：彈性力學求解新體系。大連理工大學出版社，1995.

[8]姚偉岸，鐘萬勰：辛彈性力學。高等教育出版社，2003.

[9]鐘萬勰：力、功、能量與辛數學，3rd ed.大連理工大學出版社，2012.

[10]鐘萬勰：應用力學對偶體系。科學出版社，2002.

[11]鐘萬勰：應用力學的辛數學方法。高等教育出版社，2006.

[12]鐘萬勰，高強，彭海軍：經典力學-辛講，大連理工大學出版社，2013.

[13]鐘萬勰，歐陽華江，鄧子辰：計算結構力學與最優控制。大連理工大學出版社，1993.

[14]鐘萬勰，吳志剛，譚述君：狀態空間控制理論與計算。科學出版社，2007.

[15]希爾伯特：數學問題。大連理工大學出版社，2009.

[16]H.Goldstein, Classical mechanics, 2nd ed. Addison-Wesley, 1980.

[17]G Zhong,JE Marsden.Lie-Poisson Hamilton-Jacobi theory and Lie-Poisson integrators.Physics Letter A, 113(3)∶ 134-139,1988.

[18]高強，鐘萬勰：Hamilton系統的保辛-守恒積分算法。動力學與控制學報, 2009.

[19]鐘萬勰，高強：辛破繭。大連理工大學出版社，2012.

[20]E.Hairer&G.Wanner∶Solving ordinary differential equations II-stiff and differential- algebraic problems 2nd ed.ch.7.[M],Springer,Berl in,1996.

[21]鐘萬勰，姚證：時間有限元與保辛，機械強度，27(2)，2005.

[22]Q.Gao,S.J.Tan,H.W.Zhang,W.X.Zhong.Symplectic algorithms based on the principle of least action and generating functions. International Journal for Numerical Methods in Engineering, 89(4)∶ 438-508,2012.

[23]高強，鐘萬勰：Hamilton系統的保辛-守恒積分算法。動力學與控制學報, 2009.

[24]鐘萬勰，高強∶ 約束動力系統的分析結構力學積分，動力學與控制學報，4(3), 2006.

[25]杜瑞芝：數學史詞典。山東教育出版社，2000.

[26]阿蒂亞：數學的統一性。大連理工大學出版.

相關鏈接：

祖沖之（公元429年4月20日～公元500年）是我國杰出的數學家，科學家。南北朝時期人，漢族人，字文遠。生于宋文帝元嘉六年，卒于齊昏侯永元二年。祖籍范陽郡遒縣（今河北淶水縣）。為避戰亂，祖沖之的祖父祖昌由河北遷至江南。祖昌曾任劉宋的“大匠卿”，掌管土木工程；祖沖之的父親也在朝中做官。祖沖之從小接受家傳的科學知識。青年時進入華林學省，從事學術活動。一生先后任過南徐州（今鎮江市）從事史、公府參軍、婁縣（今昆山市東北）令、謁者仆射、長水校尉等官職。其主要貢獻在數學、天文歷法和機械三方面。

祖沖之算出π的真值在3.1415926和3.1415927之間，相當于精確到小數第7位，這一紀錄直到1427年才被阿拉伯學者卡西所超越。他提出約率22/7和密率355/113，這一密率值是世界上最早提出的，比歐洲早1100年，所以有人主張叫它“祖率”，直到16世紀年才由荷蘭人奧托得到。

（以上下載自：百度百科http∶//baike.baidu.com/ view/2122.htm）

祖沖之計算出π的真值在3.1415926和3.1415927之間，如果用割圓術計算，需要計算圓內接正49152邊形才可以達到3.14159265，而以當年的計算工具顯然做不到。祖沖之到底用什么方法計算到這么精確的值現在已經是個謎。