作為解方程的基本通用算法,牛頓法是應(yīng)用數(shù)學(xué)和計算數(shù)學(xué)最重要的算法之一。這個簡單神速的算法被冠以牛頓大名,那它真是牛頓發(fā)現(xiàn)的嗎?
傳說三千五百年前,聰明的巴比倫人想出了一個計算平方根簡單漂亮的竅門:如果一個正數(shù)比 小一點,則A除以它就比 大一點,這樣它們的平均值就有希望更靠近 。于是,先取 的一個大致估計記為a, 只需求 a 和的平均,寫成公式就是
算出的 b 就是 更精確的近似值 。
神奇的是,如果你的初始估計不是太差,比如猜到作為 a 有兩位數(shù)字精確,那么 b≈1.414 就有四位數(shù)字精確。用同樣辦法改進(jìn) b 得到 c≈1.4142136 的精確數(shù)字就有八位。每次計算精確數(shù)位翻一番!
猜不出初始近似怎么辦呢?沒關(guān)系,從任何非零正數(shù) a 出發(fā)都可以。隨意亂猜的結(jié)果無非是多算幾步,很快就會進(jìn)入精度加倍狀態(tài)。這個方法流傳至今成為經(jīng)典傳奇,叫做巴比倫方法。為什么說是個傳奇?因為沒有文字記載。直到公元60年,古希臘數(shù)學(xué)家Hero of Alexandria (c. 10-c.70) 對這個方法才給出了有案可查的第一個明確的表述,所以巴比倫法也稱為赫倫方法 (Heron’s method)?,F(xiàn)在我們知道, 本文主題牛頓迭代法最古老的源頭來自四大文明之一的巴比倫文明。
什么是牛頓迭代法?
牛頓 (Isaac Newton,1643-1727) 的大名就無需介紹了。他以發(fā)現(xiàn)力學(xué)三大定律、萬有引力定律和微積分躋身從古到今最偉大的科學(xué)家行列。牛頓無疑是有史以來最杰出的數(shù)學(xué)家之一。作為一個重大貢獻(xiàn),他率先發(fā)現(xiàn)的微積分可以上升到人類文明瑰寶的高度。所有的理工科大學(xué)生都應(yīng)該知道赫赫有名的“牛頓法”,也稱“牛頓迭代”或“牛頓近似”。自然科學(xué)家、應(yīng)用和計算數(shù)學(xué)家及工程學(xué)家們一旦需要求解非線性方程和方程組,腦子里首先應(yīng)該想到的就會是牛頓法。
什么是牛頓法呢?設(shè)想我們要求出一元非線性方程f(x) = 0的解,比如說x – cos x = 0,這里f(x) =x – cos x 。數(shù)學(xué)史上有個著名的阿貝爾不可能定理,說的是非線性方程一般來說是不可能保證找到精確解的,門都沒有。所以我們需要所謂“數(shù)值方法”一步一步地逼近解,算到精度夠了就行。假如f(x)有導(dǎo)函數(shù)f'(x) ,牛頓法就是這樣的迭代程式:先取一個初始點x0作為解的近似,然后按下面的簡單公式依次迭代:
(1)就得到一個序列x0, x1, x2, … 。只要滿足三個并不苛刻的條件:(i)函數(shù) f(x) 二次連續(xù)可導(dǎo),(ii)在所求解x*處導(dǎo)數(shù)非零,加上(iii)初始近似x0足夠接近x*,則這個序列將快速趨近 ,基本上是每一步精確數(shù)位加倍。因此,實際計算中大多三五步迭代就可以獲得足夠精確的近似解。讀者可以很容易驗證,巴比倫方法其實就是求解平方根方程x2-A=0的牛頓迭代。當(dāng)然啦,巴比倫人不大可能知道“迭代”這個概念。
你要作科學(xué)計算和工程計算嗎?幾乎所有問題要么本身就是個方程,要么一定會在某個步驟需要解方程。作為解方程的基本通用算法,牛頓法是應(yīng)用數(shù)學(xué)和計算數(shù)學(xué)最重要的算法之一。這個簡單神速的算法被冠以牛頓大名,那它真是牛頓發(fā)現(xiàn)的嗎?這是個歷史悠久又頗有爭議的傳奇故事,其間包含數(shù)學(xué)史上一個個如雷貫耳的名字。
經(jīng)久不衰的“謬誤”
現(xiàn)在大家用的代數(shù)符號和表達(dá)式體系,創(chuàng)始人是法國人韋達(dá) (Fran?ois Viète,1541-1603) 。他的謀生職業(yè)是律師,他還做過亨利三世和四世的王室智囊,掙足了錢給自己提供經(jīng)費(fèi)研究數(shù)學(xué)并將結(jié)果出版。除了代數(shù)上的造詣,他還是方程理論的大師。他計算能力超強(qiáng),在歐洲首次算出十位精確圓周率值。韋達(dá)在十七世紀(jì)初提出了一個多項式方程求根的算法,每一次計算精度數(shù)位增加一位。用現(xiàn)在的話說叫線性收斂或一階收斂。韋達(dá)的算法解釋起來很費(fèi)功夫,有興趣的讀者可以參考引文[1] 。
據(jù)文獻(xiàn)[2]考證,牛頓于1664年左右讀到韋達(dá)的技巧,約1669年寫入《分析論》(De analysi)。但這部書直至1711年才由威爾士數(shù)學(xué)家瓊斯(William Jones,1675-1749)為他編輯出版。牛頓在書中改進(jìn)了韋達(dá)的思路,提出一個近似求解多項式方程的新方法。以三次方程x3 – 2x – 5 = 0為例。他首先注意到在2與3之間有個解(讀者可以用介值定理驗證),于是他把這個解寫成x = 2 + p, 代入原方程化簡后得到p的三次方程p3 + 6 P2 + 10p – 1 = 0。當(dāng)然,解這個新方程看起來跟老方程一樣困難。但p的方程可以用上微積分的思路求解:因為p很小,它的平方和立方就更小,于是三次函數(shù)p3 + 6 p2 + 10p – 1 可以用線性部分 10p – 1近似。解10p-1=0 得到 p ≈ 0.1。也就是 x 的近似從 2 到 2.1, 精確數(shù)位翻倍。
然后,牛頓依法炮制,即令p = 0.1 + q,代入0 = p3 + 6 p2 + 10p – 1:
0= (0.1+q)3+6(0.1+q)2+10(0.1+q)-1=q3+6.3q2+11.23q+0.061≈11.23q+0.061
得到q ≈ -0.0054,也就是x=2.1+q≈2.0956 。再令q = -0.0054 + r,同法得r ≈-0.00004852。這樣經(jīng)過三步以后,牛頓找到原方程的一個8位精確的近似解:
x = 2 + p + q + r ≈ 2 + 0.1 – 0.0054 - 0.00004852 = 2.09455147
每算一步精確數(shù)位加倍!
在其1687年首版的輝煌巨著《自然哲學(xué)的數(shù)學(xué)原理》中,牛頓求解了用于天文學(xué)的x – e sin x = M這個超越方程。他將其中的正弦函數(shù)用級數(shù)展開,得到一個近似的多項式方程。然后他就可以用上述數(shù)值逼近多項式方程解的法子得到原方程的近似解。正如他在二十年前所做的那樣,他沒有明確地用到函數(shù)的導(dǎo)數(shù)概念來推導(dǎo)這個數(shù)值方法,也沒有明確提出迭代概念和公式。
這就是后人所知牛頓參與創(chuàng)立“牛頓法”的過程。牛頓的貢獻(xiàn)是用微積分思路,在韋達(dá)方法的基礎(chǔ)上把巴比倫方法從平方根方程 x2-A=0 推廣到一般的多項式求根。
英國科學(xué)史專家Nicholas Kollerstrom于1992年發(fā)表了一篇關(guān)于牛頓法的考證文章[3]。文章的標(biāo)題很有意思:《托馬斯·辛普森和“牛頓近似法”:一個經(jīng)久不衰的迷思》。意思是說把公式(1)稱為“牛頓法”是個迷思(myth),也就是一個廣泛流傳的謬誤,而且這個謬誤“經(jīng)久不衰”。他指出,牛頓法(1)有兩個重要的特征:1. 它是一個迭代過程;2. 它采用了微分表達(dá)式。而這兩個特征中的哪一個,都沒有在牛頓的《分析論》里出現(xiàn)。迭代法在理論上是一個無窮極限過程,牛頓只給出了三步計算演示。其實還可以加上一條:由于沒有使用微分,牛頓提出的方法只能用于多項式,不是一個通用算法。
第一個真正的迭代法
數(shù)值計算求解方程的第一個真正意義上的迭代法是跟牛頓同在英國的約瑟夫·拉夫森(Joseph Raphson,1648-1715)在他于1690年發(fā)表的文章《方程分析通論》中給出的近似方法。但它同樣沒有求導(dǎo)運(yùn)算,因此不符合“牛頓法”的第二個特征。然而,一些慷慨的后人,包括部分現(xiàn)代數(shù)學(xué)家,把“牛頓法”的勛章切成一半分給了拉夫森——稱之為“牛頓-拉夫森法”。
那么,拉夫森是怎么獲得他的數(shù)值方法的呢?我們用求解三次方程x3 – ax + b = 0來描述他的求解方案。在每次迭代中,他分兩步走。設(shè)目前的近似解為u,則將下一個的近似解寫成u + d。然后用x=u + d代入方程并按二項式公式展開,這是第一步。在第二步,合并同類項得到d的一次項的系數(shù)3u2 – a,然后令 ,這樣得到下一個近似解。
拉夫森強(qiáng)調(diào)用他的上述辦法周而復(fù)始地迭代下去,就可以計算出滿足任意精確的方程解。然而我們依然看不到求導(dǎo)數(shù)運(yùn)算的影子。此外,他僅僅對多項式方程提出了這個迭代法,用到的二項式公式無法直接推廣到像求解超越方程這樣的情形。在他最初由倫敦皇家學(xué)會發(fā)表的那篇文章的前言里,他提到他的方法與牛頓之前的做法有類似之處。然而,七年后的1697年,當(dāng)他把這個方法著書時沒提牛頓的名字,而說韋達(dá)是他的方法的始祖。
如果我們比較牛頓和拉夫森的做法,不難發(fā)現(xiàn),牛頓用到一個經(jīng)過代入步驟而導(dǎo)出的一個似乎更加復(fù)雜的多項式,再丟掉高階項求得近似;拉夫森從頭到尾都是用給定的原多項式,運(yùn)算要簡單得多。拉夫森感覺自己的方法跟牛頓是完全不一樣的推導(dǎo),無需歸功于牛頓。類似的比較也陸續(xù)出現(xiàn)。比如,在1796年發(fā)表的文章《關(guān)于拉夫森先生的方法的觀察》中,作者費(fèi)倫德(W. Frend)比較了兩法的各自優(yōu)點:“考慮到兩種方法的簡單性和概念性......我認(rèn)為總的來說,拉夫森先生求解方程的方法比艾薩克·牛頓爵士的更為方便。”
在1798年,法國大數(shù)學(xué)家拉格朗日(Joseph-Louis Lagrange,1736-1813)發(fā)表了頗具影響力的論文《數(shù)值求解方程》。他精細(xì)化并推廣了牛頓著作《分析論》中的方法,但依然沒有用到導(dǎo)數(shù)或微分術(shù)語。
眼尖的讀者很快就發(fā)現(xiàn),在上述拉夫森得到的d的分?jǐn)?shù)表達(dá)式中,分子就是函數(shù)x3 – ax + b在當(dāng)前迭代點u的值,而分母恰恰就是這個函數(shù)的導(dǎo)數(shù)在此迭代點的負(fù)值,因而這就給出了“牛頓法”在多項式函數(shù)的全部內(nèi)容。然而,不能因此就說拉夫森發(fā)明了今日所稱的牛頓法!原因就在于他和牛頓一樣,都沒有使用導(dǎo)數(shù)的記號和運(yùn)算而得出一般的牛頓法格式,仍然無法直接應(yīng)用到一般的非線性方程。
被忘卻的發(fā)明人
那么,誰才是“牛頓法”當(dāng)之無愧的發(fā)明人呢?
此人的全名是托馬斯·辛普森(Thomas Simpson,1710-1761),他是比牛頓和拉夫森遲了幾十年的英國數(shù)學(xué)家。他就是近似數(shù)值積分著名的“辛普森法則”的那個辛普森。有意思的是,辛普森在牛頓法貢獻(xiàn)恐怕最大,卻被后人差不多忘得一干二凈。他反而在數(shù)值積分法獲得并非實至名歸的榮譽(yù)。該得的沒得到,不該得的反而拿著了。早他一百年,德國天文學(xué)家開普勒(Johannes Kepler,1571-1630)就已經(jīng)發(fā)現(xiàn)了近似計算“曲邊矩形”面積的該項法則。因此,德國人把我們叫慣了的辛普森法則自豪地稱作為“開普勒的桶法則”,就像我們常常把關(guān)于二項展開式各項系數(shù)的“帕斯卡三角形”稱為“楊輝三角形”那樣異曲同工。
辛普森構(gòu)造出現(xiàn)代意義下的牛頓法是在1740年,此時牛頓已經(jīng)去世了十三年。那年他出了一本關(guān)于數(shù)學(xué)的論文集,其中一篇描述了“求解方程的一個新方法”,卻沒有列出任何先驅(qū)者的姓名。在前言部分,他斷言:“因為它比以往的任何方法都更普遍,它不能不具有相當(dāng)大的用途。”這聽上去口氣很大。他是自信而非吹牛:“取給定方程的流數(shù)……”此處,流數(shù)的英文是fluxion,正是牛頓當(dāng)年用來表示今天我們所稱的“導(dǎo)數(shù)”的那個東西——函數(shù)的瞬時變化率。接著,他給出了和上面公式(1)式實際一模一樣的迭代程序,除了沒有采用當(dāng)今標(biāo)準(zhǔn)的、微積分另一發(fā)明人萊布尼茨(Gottfried Leibniz,1646-1716)所引進(jìn)的導(dǎo)數(shù)記號。
辛普森用這個普遍方法做了五個例子,包括求解三次方程、平方根計算、指數(shù)方程等。更進(jìn)一步,他第一個將他的方法用于求解含有兩個未知數(shù)和兩個方程的方程組!既然他是有史以來第一個完整地提出和今日所指的牛頓法有完全相同格式的迭代法,數(shù)學(xué)史專家Kollerstrom得出結(jié)論:辛普森才是牛頓法的發(fā)現(xiàn)者。
辛普森版的牛頓法跟現(xiàn)代教科書的差別僅僅是所用的符號。他應(yīng)當(dāng)之無愧地被授予創(chuàng)造該法的榮譽(yù)。然而,到底是誰寫出了現(xiàn)代形式的牛頓法呢?
他就是在近代數(shù)學(xué)向前邁步的崎嶇道路上留下巨大腳印的傅里葉(Joseph Fourier,1768-1830)。這位法國數(shù)學(xué)家在辛普森提出標(biāo)準(zhǔn)的牛頓法后的第二十八個年頭才出生。在十九世紀(jì)初,他首次用當(dāng)今世界通用的導(dǎo)數(shù)記號f’重新敘述了迭代法(1),同時把它說成是“牛頓法”。由于拉格朗日的那篇雄文,后來,有些英國數(shù)學(xué)家將此法稱為“牛頓和拉格朗日的方法”,而對拉夫森只字不提。十九世紀(jì)的數(shù)學(xué)史名家、德國人康托爾(Moritz Cantor,1829-1920)考察了牛頓、拉夫森等人的方程近似求解法,把拉夫森描繪為“牛頓的絕對仰慕者和模仿者”,認(rèn)為他的近似法“與牛頓的方法極其類似”。
瑞士出生的美國數(shù)學(xué)史家弗洛里安·卡喬里(Florian Cajori,1859-1930)在1911年的《美國數(shù)學(xué)月刊》上發(fā)文提出這個方法理應(yīng)被稱為“牛頓-拉夫森法”。但是,他的命名論據(jù)受到Kollerstrom的質(zhì)疑,依據(jù)正是那個“兩個特征”,后者認(rèn)為榮譽(yù)只能歸于辛普森。然而,著名的數(shù)學(xué)史家博耶(Carl Boyer,1906-1976)在他1968年出版的大作《數(shù)學(xué)史》中這樣斷言:“方程近似求解的牛頓法可在《分析論》中發(fā)現(xiàn)。”由于牛頓在科學(xué)史包括數(shù)學(xué)史上的巨大名望,擁有“牛頓法”的真正主人辛普森在數(shù)學(xué)史上幾乎失去了立足之地,拉夫森也只能偶然出現(xiàn)在牛頓名字的后面。
這種張冠李戴的命名在科學(xué)和數(shù)學(xué)史中比比皆是。例如,學(xué)過初等微積分的人都知道求不定式極限的“洛必達(dá)法則”,實際上是瑞士數(shù)學(xué)家伯努利(Johann Bernoulli,1667-1748) 的杰作。伯努利的數(shù)學(xué)功力可不是法國數(shù)學(xué)家洛必達(dá)(Guillaume de l’H?pital,1661-1704) 所能望其項背的。洛必達(dá)于1694年3月17日在給伯努利的信中,提出每年給他三百法郎換取他的最新數(shù)學(xué)發(fā)現(xiàn),并且不能透露給第三者。當(dāng)年伯努利告訴了洛必達(dá)這個求極限定理,兩年后洛必達(dá)將它寫進(jìn)了自己的著作《曲線的無窮小分析》,據(jù)說這是全世界的第一本微積分教材。盡管洛必達(dá)在書中感謝了萊布尼茨和伯努利,尤其感謝伯努利,或許作者有意無意地沒有明確承認(rèn)“洛必達(dá)法則”是誕生于別人家的嬰兒,不明就里的后人就把這條極其有用的求極限法則冠上了他的名字。當(dāng)然伯努利在數(shù)學(xué)史上大名鼎鼎,少了個洛必達(dá)法則也不至于淪落成籍籍無名的歷史過客。
牛頓法的基本思想和深層含義是什么?哪些激情探索者又續(xù)寫傳奇?且待下篇
后記:作者以此文紀(jì)念我們的導(dǎo)師李天巖教授(1945年6月28日-2020年6月25日)逝世一周年。
參考文獻(xiàn)
[1] Tjalling J. Ypma, Historical Development of the Newton-Raphson Method, SIAM Review, 37, 531-551, 1995
[2] P. Deuflhard,“A short history of Newton’s method,”Documenta Math.,Extra Volume ISMP,25,25-30,2012.
[3] N. Kollerstrom,“Thomas Simpson and 'Newton’s method of approximation’:an enduring myth,”BJHS,25, 347-354,1992.
聯(lián)系客服