鐘萬(wàn)勰

大連理工大學(xué)

工科數(shù)學(xué)
——祖沖之思路

鐘萬(wàn)勰

大連理工大學(xué)

前 言

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

對(duì)此難以認(rèn)同，不平則鳴么！作者不是數(shù)學(xué)家，難免有所舛誤，有錯(cuò)就請(qǐng)批評(píng)，不用客氣。

作者這些年在動(dòng)力學(xué)方面努力，尤其是有約束的動(dòng)力學(xué)計(jì)算分析，發(fā)現(xiàn)作為機(jī)器人動(dòng)力學(xué)基礎(chǔ)的約束動(dòng)力系統(tǒng)微分-代數(shù)方程(DAE, Differential-Algebraic Equation)的求解，例如現(xiàn)代著名著作

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的求解，效果不理想。一些著名洋人軟件，例如廣泛應(yīng)用的ADAMS，其數(shù)值結(jié)果也不理想。因?yàn)檫@些著作的求解方法，是先進(jìn)行微商，將約束方程歸化到微分方程，其微商次數(shù)稱為Index，可稱為Index法?？雌饋砑s束條件處處滿足，而實(shí)際上數(shù)值結(jié)果的約束條件滿足不行。按“實(shí)踐是檢驗(yàn)真理的唯一標(biāo)準(zhǔn)”，不行。

鐘萬(wàn)勰，中國(guó)科學(xué)院院士，大連理工大學(xué)工程力學(xué)研究所所長(zhǎng)、教授。工程力學(xué)、計(jì)算力學(xué)專家。1956年同濟(jì)大學(xué)橋梁與隧道系畢業(yè)。1993年當(dāng)選為中國(guó)科學(xué)院院士。

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

中國(guó)著名南北朝數(shù)學(xué)家祖沖之(429～500)在計(jì)算圓周率π時(shí)，已經(jīng)達(dá)到π＝3.1415926～3.1415927之間，唐朝魏征等撰寫的史書《隋書》中有記載∶

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

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

如所周知，π是現(xiàn)代數(shù)學(xué)不可回避的基礎(chǔ)，看到中國(guó)祖師爺?shù)闹匾暙I(xiàn)，欣喜不已。于是我們就探討祖師爺祖沖之在當(dāng)年條件下是怎么計(jì)算的。由此引申下去，其實(shí)許多數(shù)學(xué)基本概念，例如無窮序列，極限等等，中國(guó)祖師爺也是明白的。并且在實(shí)際上應(yīng)用了，按此路線探究下去就接地氣。尤其是接中國(guó)人的地氣。

所以，中國(guó)數(shù)學(xué)并非一無所成，而是挖掘不夠。竊以為，既然在中國(guó)大學(xué)講數(shù)學(xué)，就不應(yīng)將中國(guó)數(shù)學(xué)祖師爺?shù)墓ぷ骱雎裕鴳?yīng)有所傳承。我們不可“光說不練”，本文就從祖沖之如何計(jì)算圓周率π講起，加以傳承，貫通古今，融合中西。中國(guó)數(shù)學(xué)也應(yīng)占有一席之地的?？上葟闹袊?guó)工科數(shù)學(xué)教學(xué)做起。

2014年5月，鐘萬(wàn)勰參加了在廣東工學(xué)院的工科數(shù)學(xué)教材會(huì)議，介紹了祖沖之類算法。希冀大學(xué)數(shù)學(xué)教材要加入中國(guó)元素。本文是以實(shí)際工作貫徹該意圖。不僅是講講而已，而要有實(shí)際行動(dòng)的。期待官方能任領(lǐng)導(dǎo)改革之職。下面就由探討祖沖之的工作開始。

著作[10]的序言中講：“1999年5月，教育部委托上海大學(xué)在錢偉長(zhǎng)教授主持下召開了一次應(yīng)用力學(xué)教改的會(huì)議。該會(huì)議使我下決心寫出這本書，為此花費(fèi)了大量的精力?！北砻麋娙f(wàn)勰很早就有志于教學(xué)。雖然曾努力推動(dòng)，但思路不得法，效果不理想?，F(xiàn)在再一次努力，希望能發(fā)揮應(yīng)有作用。但鐘萬(wàn)勰困惑：教育改革沒有SCI，得不到評(píng)價(jià)，又怎能鼓勵(lì)青年教師呢。官僚機(jī)制厲害呀！當(dāng)官只要數(shù)數(shù)SCI就可以了，懶政！此“官不任事”者也。面對(duì)這種評(píng)價(jià)機(jī)制，只能長(zhǎng)嘆：“此天之亡我，非戰(zhàn)之罪也！”，對(duì)此唏噓不已。

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

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

中國(guó)在東周時(shí)期，已經(jīng)發(fā)明了算盤，并且申了遺。算盤是古代的計(jì)算機(jī)。并且中國(guó)勾、股、弦的定理也早已發(fā)現(xiàn)，周公時(shí)期(公元前1000多年)已經(jīng)有所記載，稱商高定理。古代希臘有畢達(dá)哥拉斯者，對(duì)此有系統(tǒng)論述，所以大家說這是畢達(dá)哥拉斯定理。我們可仍稱呼為商高定理。因?yàn)樽鏇_之未必知道有畢達(dá)哥拉斯者，他實(shí)際可使用中國(guó)傳承下來的商高定理。

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

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

圖1.割圓術(shù)

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

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

再用商高定理。以伸長(zhǎng)線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條邊的長(zhǎng)度si+1=2li+1相加在一起而已。

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

內(nèi)接正多角形的總邊長(zhǎng)，永遠(yuǎn)趕不上2π，因?yàn)檫B接線是直線短程線；然而可不斷逼近2π，得到序列。所以，得到的無窮序列是無限逼近于2π的，2π就是其極限。具體數(shù)值計(jì)算給出

表1：不同正多邊形對(duì)應(yīng)的圓周率

表1為用SiPESC集成平臺(tái)調(diào)用的Python編程計(jì)算得到的表格，可見當(dāng)取正16384邊形時(shí)，近似計(jì)算得到圓周率為.14159263，此過程不斷進(jìn)行，得到了無窮序列，其極限就是圓周率π。

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

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

科學(xué)網(wǎng)2015-5-24傳達(dá)韓啟德說：“中國(guó)科學(xué)家缺位科普”。在工科數(shù)學(xué)中，中國(guó)數(shù)學(xué)家也是缺位的。讓我們來改變這個(gè)不合理的現(xiàn)狀吧。

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

中國(guó)數(shù)學(xué)如自動(dòng)缺位，前途堪憂。

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

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

數(shù)學(xué)大師馮?諾伊曼指出：“能以更好的形式推廣為更有效的新理論的理論將戰(zhàn)勝另一理論?！仨殢?qiáng)調(diào)的是，這并不是一個(gè)接受正確理論、拋棄錯(cuò)誤理論的問題。而是一個(gè)是否接受為了正確的推廣而表現(xiàn)出更大的形式適應(yīng)性的問題”，見文獻(xiàn)[1]。

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

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

黨中央推動(dòng)改革。官員有公權(quán)力，希望能有擔(dān)當(dāng)，領(lǐng)導(dǎo)我們改掉不合理的教材吧。希望有關(guān)部門能予以關(guān)注，改革3的路就在腳下，難道工科數(shù)學(xué)就不用“認(rèn)祖歸宗”嗎？ 振作起來干吧！

二，無約束動(dòng)力學(xué)的時(shí)間積分

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

經(jīng)典力學(xué)眾多著作是許多著名數(shù)學(xué)家努力耕耘的成果。內(nèi)容豐富、數(shù)學(xué)理論優(yōu)美，難以窮盡。這里不過是取其中心部分內(nèi)容講解，畢竟是教數(shù)學(xué)么。然而對(duì)于Hamilton體系的認(rèn)識(shí)，卻有曲折過程。

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

量子力學(xué)當(dāng)年面臨的一個(gè)關(guān)鍵問題是光譜分析。觀察到的譜線分裂現(xiàn)象，被Hilbert的大弟子，大數(shù)學(xué)家Hermann Weyl用群論分析清楚[2]。進(jìn)入Hamilton狀態(tài)空間后，H.Weyl在1939年又提出Symplectic對(duì)稱[3]。華羅庚先生音譯為‘辛’。H.Weyl提出了Hamilton正則方程的辛對(duì)稱，表明了動(dòng)力系統(tǒng)運(yùn)行的基本性質(zhì)，表征了動(dòng)力學(xué)的辛階段。

辛，是大數(shù)學(xué)家Hermann Weyl在名著[3]中提出的。作者寫道：“The name ‘complex group’ formerly advocated by me in allusion to line complexes, …h(huán)as 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’ ”

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

H.Weyl的工作吸引了純數(shù)學(xué)家們高度關(guān)注。從微分幾何的角度出發(fā)，給出了“辛幾何”的定義；但其理論過分高深以致讓工程師難以接受。今天已經(jīng)是信息時(shí)代，計(jì)算機(jī)的發(fā)展使人們認(rèn)識(shí)到“計(jì)算科學(xué)”的重要性。“計(jì)算科學(xué)與理論、實(shí)驗(yàn)共同構(gòu)成現(xiàn)代科學(xué)的3大支柱”的論點(diǎn)得到了廣泛認(rèn)同。計(jì)算科學(xué)是當(dāng)代數(shù)學(xué)的大發(fā)展。計(jì)算科學(xué)的主要特點(diǎn)是離散。計(jì)算機(jī)本身就是離散的么。

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

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

辛代數(shù)：容易理解、掌握?；貧w到力學(xué)，離散系統(tǒng)?？蓮幕⒖硕芍v清楚，見[9]。自然科學(xué)基金1993年支持了課題(19372011)“結(jié)構(gòu)動(dòng)力學(xué)及其辛代數(shù)方法”。辛代數(shù)的提法不是現(xiàn)在才有的。

我國(guó)數(shù)學(xué)家馮康率先提出[4]，動(dòng)力學(xué)積分的差分格式應(yīng)保辛，差分格式就是離散積分。國(guó)外著作[5]也接受了。

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

作者隨后在著作[9～12]，建立了分析動(dòng)力學(xué)與分析結(jié)構(gòu)力學(xué)的模擬理論。會(huì)同以前建立的計(jì)算結(jié)構(gòu)力學(xué)與最優(yōu)控制間的模擬理論[13,14]，將多個(gè)領(lǐng)域皆奠基于狀態(tài)空間的數(shù)學(xué)理論。于是各個(gè)領(lǐng)域的數(shù)學(xué)計(jì)算方法可以互通有無了。基礎(chǔ)就是狀態(tài)空間法的辛對(duì)稱。

1900年在第二次世界數(shù)學(xué)大會(huì)上，D.Hilbert做了一個(gè)著名報(bào)告《數(shù)學(xué)問題》，深刻影響了20世紀(jì)數(shù)學(xué)的發(fā)展，指出：“只要一門科學(xué)分支能提出大量的問題，它就充滿著生命力；而問題缺乏則預(yù)示著這門學(xué)科獨(dú)立發(fā)展的衰亡或終止?！蹦芸吹絾栴}就好，表明可繼續(xù)發(fā)展。辛指出Hamilton正則方程有對(duì)稱性，是經(jīng)典力學(xué)繼續(xù)發(fā)展的機(jī)會(huì)，有大量的新問題?！稊?shù)學(xué)問題》中提出23個(gè)問題，其中第23號(hào)是變分法的進(jìn)一步發(fā)展。Hilbert說：“我已經(jīng)廣泛地涉及了盡可能是確定的和特殊的問題…，用一個(gè)一般的問題來做結(jié)束…我指的是變分法”，見文獻(xiàn)[15]。變分法不單純是一個(gè)數(shù)學(xué)問題，而是一個(gè)方向，是大師的遠(yuǎn)見卓識(shí)。

變分原理的提出，由來已久?！按笞匀豢偸亲咦钊菀缀妥羁赡艿耐緩健保@是費(fèi)爾馬(Fermat)著名的自然哲學(xué)原理。1744年，J.Bernoulli提出了“最速下降線”問題，大體上可認(rèn)為是數(shù)學(xué)變分法的開始。以后蓬勃發(fā)展，Euler-Lagrange方程，繼而總結(jié)為Hamilton變分原理。分析動(dòng)力學(xué)與相應(yīng)的常微分方程理論的成功，自然要發(fā)展到偏微分方程。到位勢(shì)理論，有Laplace方程的求解，Green,Gauss,Dirichlet等先后指出，可將其轉(zhuǎn)化為變分原理。Riemann(1826～1866)將其命名為Dirichlet原理，屬于橢圓型偏微分方程理論，本來在蓬勃發(fā)展。然而在1870年發(fā)生了曲折，強(qiáng)調(diào)數(shù)學(xué)嚴(yán)格性的維爾斯特拉斯(Weierstrass)否定了Dirichlet原理；然而數(shù)學(xué)物理中許多重要結(jié)果都依賴于此原理而建立。1899年,Hilbert用邊界條件的光滑化保證了極小化函數(shù)的存在，從而挽救了Dirichlet原理。表明變分原理經(jīng)歷了“鳳凰涅槃，浴火重生”的過程。

此后一個(gè)多世紀(jì)以來，變分原理有了巨大進(jìn)展，從有限元法發(fā)展出來的“計(jì)算科學(xué)”也是變分法的發(fā)展；而辛也屬于變分法的發(fā)展。我國(guó)在變分原理的研究方面是有成就的，胡海昌的三類變量變分原理蜚聲世界。有限元法離散取得了巨大成功，有限元單元推導(dǎo)的基礎(chǔ)就是變分原理。動(dòng)力學(xué)變分原理稱為最小作用量原理。做一點(diǎn)解釋如下。

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

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

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

既然分析求解困難，則只能離散數(shù)值求解，差分法是常用的。我國(guó)數(shù)學(xué)家馮康提出，差分格式應(yīng)保辛，棋高一著。設(shè)時(shí)間劃分為步長(zhǎng)η的一系列節(jié)點(diǎn) to=0,…,tk=kη,…。逐步積分要在積分出tk-1的狀態(tài)后，進(jìn)一步積分時(shí)間區(qū)段k#:[tk-1,tk]，以計(jì)算tk的狀態(tài)。此時(shí)可用最小作用量原理。時(shí)間區(qū)段k#:[tk-1,tk]的作用量是

作用量可從例如[16]找到。動(dòng)力系統(tǒng)通常是沒有約束條件的。

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

然后，按一般理論，有

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

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

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

一般的動(dòng)力學(xué)積分是沒有約束的，相對(duì)比較簡(jiǎn)單些。而祖沖之類算法并不要求區(qū)段k#:[tk-1,tk]內(nèi)滿足約束，而只要求在離散時(shí)間點(diǎn)處滿足約束條件。

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

這方面有許多數(shù)值例題，不再重復(fù)。本節(jié)是初步的分析動(dòng)力學(xué)的內(nèi)容，時(shí)間有限元么。對(duì)于祖沖之類算法，不管約束條件而取“動(dòng)力學(xué)意義下的短程線”來說，是不可缺少的基礎(chǔ)。有了基礎(chǔ)，就可講動(dòng)力學(xué)的祖沖之類算法了。

三，祖沖之類算法

現(xiàn)代化推進(jìn)使機(jī)器人成為發(fā)展的熱點(diǎn)。機(jī)器人要運(yùn)動(dòng)，當(dāng)然需要?jiǎng)恿W(xué)分析。機(jī)器人起碼是機(jī)構(gòu)，機(jī)械工程離不開機(jī)構(gòu)。機(jī)器各個(gè)部件皆可運(yùn)動(dòng)，但相互間是有位移幾何約束的，稱為完整約束。因此我們面臨受約束的動(dòng)力學(xué)系統(tǒng)分析問題。國(guó)外發(fā)展了稱為微分-代數(shù)方程(DAE,Differential-Algebraic Equation)的求解方法。DAE的求解成為一個(gè)重要問題，祖沖之算法可推廣用于此類重要問題。著作[20]的副標(biāo)題就直接寫上了stiff and differential-algebraic problems。論述了DAE的求解，但效果不理想。因?yàn)檫@些著作的求解方法，是先進(jìn)行微商，將約束方程歸化到微分方程，其微商次數(shù)稱為Index，可稱為Index法。看起來約束條件處處滿足，而實(shí)際上數(shù)值結(jié)果的約束條件滿足不行。

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

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

采用祖沖之類算法的思路，時(shí)間積分的格點(diǎn)處，要求位移約束嚴(yán)格滿足；而在時(shí)間區(qū)段內(nèi)則不管位移約束了，要求“動(dòng)力學(xué)意義下的短程線”即可?，F(xiàn)在是動(dòng)力學(xué)，要求的短程線應(yīng)是2n維狀態(tài)空間的短程線，因?yàn)樵跁r(shí)間積分的區(qū)段內(nèi)不管位移約束條件了。

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

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

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

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

導(dǎo)出的方程組是

對(duì)應(yīng)的DAE為(這些是國(guó)外的方法論)

H(q,p)=E是保守的。此即在狀態(tài)空間q,p下的微分-代數(shù)方程[5,20]。國(guó)外求解方法論的基礎(chǔ)是微分方程，要推出聯(lián)立微分方程引起指標(biāo)(Index)問題等，給理論與計(jì)算上帶來了復(fù)雜的附加因素。

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

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

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

同樣的思路也可運(yùn)用于DAE的方程求解。先進(jìn)行離散，約束條件只要求在節(jié)點(diǎn)處嚴(yán)格滿足，而軌道則可按無約束動(dòng)力系統(tǒng)積分，(動(dòng)力學(xué)意義下的短程線)，例如可用時(shí)間有限元近似積分。只要節(jié)點(diǎn)劃分足夠密，則軌道定會(huì)逼近真實(shí)解的。這樣的思路下導(dǎo)出的算法，可稱之為祖沖之方法論、祖沖之類算法。我們理應(yīng)接住1500多年前祖師爺傳來的球，挖掘繼承，發(fā)揚(yáng)光大，與近代數(shù)學(xué)融合?！肮艦榻裼?，洋為中用”，向前沖之。

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

基于分析結(jié)構(gòu)力學(xué)理論的方法，可從最簡(jiǎn)單的例題著手闡明。

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

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

其中g(shù)r是重力加速度。

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

擴(kuò)展Lagrange函數(shù)的作用量，也要相應(yīng)地離散。區(qū)段k#:(tk-1,tk)內(nèi)按祖沖之思路，不用考慮約束，因此位移仍是n維。但到k號(hào)時(shí)間節(jié)點(diǎn)tk時(shí)，要將nc維向量λk考慮進(jìn)去。

分別對(duì)區(qū)段位移分量進(jìn)行簡(jiǎn)單的線性插值，單元內(nèi)部的約束條件則放松掉。運(yùn)用有限元法近似計(jì)算作用量。線性插值有。但還要加入時(shí)間節(jié)點(diǎn)tk處的約束對(duì)偶向量λk，于是與完全無約束的(2.2)對(duì)比，應(yīng)修改為

它相當(dāng)于結(jié)構(gòu)力學(xué)的區(qū)段變形能。整體變形能無非是全部單元變形能之和

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

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

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

根據(jù)

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

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

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

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

可歸納分析結(jié)構(gòu)力學(xué)(帶參變量λk)的算法如下：

1)形成Lagrange函數(shù)(動(dòng)能—?jiǎng)菽?，并引入約束的Lagrange參數(shù)，形成擴(kuò)展Lagrange函數(shù)。

2)將時(shí)間坐標(biāo)離散，得一系列的時(shí)間點(diǎn)to=0,t1,…,tk…，以各點(diǎn)的n維位移qk當(dāng)作未知數(shù)。

3)按分析結(jié)構(gòu)力學(xué)，計(jì)算區(qū)段(tk-1,tk)的作用量sk(qk-1,qk,λk)，可用有限元線性插值得作用量。

4)生成對(duì)偶向量

組成狀態(tài)向量。

5)根據(jù)

以及初始條件，對(duì)1號(hào)單元，用插值公式計(jì)算初始p0并組成初始狀態(tài)。

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

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

分析結(jié)構(gòu)力學(xué)的算法運(yùn)用時(shí)間有限元，不具體考慮約束對(duì)偶λ的微分方程，而是將λk逐點(diǎn)當(dāng)作待定的常參數(shù)。因不涉及其微商，故不必考慮其微分方程，也不會(huì)產(chǎn)生DAE微分方程理論的Index問題。通過數(shù)值例題可看到祖沖之方法論，祖沖之類算法特色思路的效果。本課題的數(shù)值結(jié)果展示了優(yōu)越性。

可提供更多的數(shù)值例題以供比較。本文展示約束保守體系的分析結(jié)構(gòu)力學(xué)DAE有限元保辛算法的解。

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

解：按[5]取重力加速度取gr=1，擺的長(zhǎng)度為1，質(zhì)點(diǎn)質(zhì)量為1，z坐標(biāo)以向下為正。Hamilton函數(shù)為

初始位置為

初始動(dòng)量為

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

例題3：空間雙擺問題：質(zhì)點(diǎn)1的初始位置是

初始速度是{0.1 0 0}T，質(zhì)點(diǎn)2的初始位置是

圖2.球面擺質(zhì)點(diǎn)的軌跡圖(步長(zhǎng)分別為0.03秒和0.1秒)

圖3，球面擺Hamilton函數(shù)隨時(shí)間變化，真實(shí)H=0.99680

圖4.質(zhì)點(diǎn)1的軌跡，約束條件滿足非常好

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

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

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

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

圖5.質(zhì)點(diǎn)2相對(duì)于質(zhì)點(diǎn)1的軌跡，約束條件滿足依然非常好

圖6.質(zhì)點(diǎn)2的軌跡，混沌出現(xiàn)了

圖7.Hamilton函數(shù)隨時(shí)間變化，真實(shí)H=28.85251

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

本文涉及了辛矩陣對(duì)稱群，讀者可能對(duì)群論有些怕。其實(shí)本文只用到群的最簡(jiǎn)單原理。英國(guó)M.Atiyah(1990年當(dāng)選的皇家學(xué)會(huì)會(huì)長(zhǎng)，著名純數(shù)學(xué)家)說[26]：“群在自然中產(chǎn)生，它們是使事物運(yùn)動(dòng)的東西，它們是變換或置換…。理解這些東西的本性，并且使用它們才是目的”。又說：“重要的東西常常不是技術(shù)上最困難的即最難證明的東西，而常常是較為初等的部分。因?yàn)檫@些部分與其它領(lǐng)域、分支的相互作用最廣泛，即影響面最大”?！霸谌赫撝杏性S多極端重要的，并且在數(shù)學(xué)的各個(gè)角落到處都出現(xiàn)的東西。這些是較為初等的東西：群及其同態(tài)、表示的基本觀點(diǎn)，一般的性質(zhì)，一般的方法—這些才是真正重要的”。作者也由此得益，這些觀點(diǎn)對(duì)讀者也許有益。

作為機(jī)器人基礎(chǔ)的DAE，基于祖沖之方法論得到的解，比國(guó)外著名算法的解好多了。中國(guó)祖師爺?shù)膬?yōu)秀成果應(yīng)努力予以挖掘繼承，融合近代數(shù)學(xué)，發(fā)揚(yáng)光大。中華文化博大精深，不只是懷念懷念，而應(yīng)以實(shí)實(shí)在在的東西體現(xiàn)出來。當(dāng)然還應(yīng)經(jīng)歷挖掘、品味、提煉、繼承、融合近代數(shù)學(xué)、然后才能發(fā)揚(yáng)光大。學(xué)習(xí)大學(xué)微積分，讀完后沒見到中國(guó)人實(shí)實(shí)在在的貢獻(xiàn)，遺憾；到現(xiàn)在信息時(shí)代，挖掘出祖沖之方法論，理應(yīng)占有一席之地。我們應(yīng)為此而站起來大喊大叫?！澳虾ＶT島一千多年前就在中國(guó)的管轄之下”，我們應(yīng)理直氣壯地講，因?yàn)橛惺聦?shí)依據(jù)。有人不承認(rèn)，難道我們就不敢講了嗎！今天看到祖沖之方法論的優(yōu)越性，可達(dá)到貫通古今，融合中西的境界。中國(guó)人就應(yīng)占有這一席之地，有什么不敢講的。難道這也要等著，讓洋人來發(fā)揚(yáng)光大嗎?! 他們會(huì)嗎？

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

對(duì)自己沒有信心，弱者的心態(tài)么，要不得。

當(dāng)年紅軍英勇奮斗，請(qǐng)了李德，就打敗仗。道路決定命運(yùn)么。評(píng)價(jià)體系是學(xué)術(shù)研究的指揮棒。提出一個(gè)什么SCI等的洋指標(biāo)，完全放棄了自主，又請(qǐng)來了一個(gè)“李德”，特別不爽。

中國(guó)數(shù)學(xué)界泰斗、中科院院士吳文俊認(rèn)為：“數(shù)學(xué)應(yīng)該是有助于解決實(shí)際問題的，這應(yīng)該是研究數(shù)學(xué)的主要目的，數(shù)學(xué)不是什么抽象的理論?！袊?guó)的應(yīng)用數(shù)學(xué)應(yīng)該以解決實(shí)際問題為出發(fā)點(diǎn)的，絕不能只做抽象的理論研究。”

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

[1]馮?諾依曼：數(shù)學(xué)在科學(xué)和社會(huì)中的應(yīng)用。大連理工大學(xué)出版社，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體系的辛計(jì)算格式。浙江科技出版社，2004.

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

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

[7]鐘萬(wàn)勰：彈性力學(xué)求解新體系。大連理工大學(xué)出版社，1995.

[8]姚偉岸，鐘萬(wàn)勰：辛彈性力學(xué)。高等教育出版社，2003.

[9]鐘萬(wàn)勰：力、功、能量與辛數(shù)學(xué)，3rd ed.大連理工大學(xué)出版社，2012.

[10]鐘萬(wàn)勰：應(yīng)用力學(xué)對(duì)偶體系。科學(xué)出版社，2002.

[11]鐘萬(wàn)勰：應(yīng)用力學(xué)的辛數(shù)學(xué)方法。高等教育出版社，2006.

[12]鐘萬(wàn)勰，高強(qiáng)，彭海軍：經(jīng)典力學(xué)-辛講，大連理工大學(xué)出版社，2013.

[13]鐘萬(wàn)勰，歐陽(yáng)華江，鄧子辰：計(jì)算結(jié)構(gòu)力學(xué)與最優(yōu)控制。大連理工大學(xué)出版社，1993.

[14]鐘萬(wàn)勰，吳志剛，譚述君：狀態(tài)空間控制理論與計(jì)算?？茖W(xué)出版社，2007.

[15]希爾伯特：數(shù)學(xué)問題。大連理工大學(xué)出版社，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]高強(qiáng)，鐘萬(wàn)勰：Hamilton系統(tǒng)的保辛-守恒積分算法。動(dòng)力學(xué)與控制學(xué)報(bào), 2009.

[19]鐘萬(wàn)勰，高強(qiáng)：辛破繭。大連理工大學(xué)出版社，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]鐘萬(wàn)勰，姚證：時(shí)間有限元與保辛，機(jī)械強(qiáng)度，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]高強(qiáng)，鐘萬(wàn)勰：Hamilton系統(tǒng)的保辛-守恒積分算法。動(dòng)力學(xué)與控制學(xué)報(bào), 2009.

[24]鐘萬(wàn)勰，高強(qiáng)∶ 約束動(dòng)力系統(tǒng)的分析結(jié)構(gòu)力學(xué)積分，動(dòng)力學(xué)與控制學(xué)報(bào)，4(3), 2006.

[25]杜瑞芝：數(shù)學(xué)史詞典。山東教育出版社，2000.

[26]阿蒂亞：數(shù)學(xué)的統(tǒng)一性。大連理工大學(xué)出版.

相關(guān)鏈接：

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

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

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

祖沖之計(jì)算出π的真值在3.1415926和3.1415927之間，如果用割圓術(shù)計(jì)算，需要計(jì)算圓內(nèi)接正49152邊形才可以達(dá)到3.14159265，而以當(dāng)年的計(jì)算工具顯然做不到。祖沖之到底用什么方法計(jì)算到這么精確的值現(xiàn)在已經(jīng)是個(gè)謎。