知其然,更知其所以然.中國(guó)先哲

1、 樂(lè)經(jīng)良 知其然知其然, ,更知其所更知其所以然以然. . - -中國(guó)先哲中國(guó)先哲 哪里有數(shù)哪里有數(shù), ,哪里就哪里就有美有美. . - Proclus 樂(lè)經(jīng)良數(shù)學(xué)實(shí)驗(yàn)數(shù)學(xué)實(shí)驗(yàn)上海交大數(shù)學(xué)科學(xué)學(xué)院上海交大數(shù)學(xué)科學(xué)學(xué)院 的的 計(jì)計(jì) 算算 樂(lè)經(jīng)良你也許能寫(xiě)出你也許能寫(xiě)出 = 3.1415926535實(shí)際問(wèn)題實(shí)際問(wèn)題圓周率圓周率, , 我們十分熟悉的常數(shù)我們十分熟悉的常數(shù). .用用Matlab 可以求出可以求出 到到幾百位幾百位 digits(100) vpa(pi)但你會(huì)計(jì)算但你會(huì)計(jì)算 的值嗎？你又能用幾種方法計(jì)算的值嗎？你又能用幾種方法計(jì)算？ans = 3.1415926535897932384

2、62643383279502884197169399375105820974944592307816406286208998628034825342117068 樂(lè)經(jīng)良劉徽割圓法劉徽割圓法從正六邊形開(kāi)始，逐步求邊長(zhǎng)與從正六邊形開(kāi)始，逐步求邊長(zhǎng)與面積面積oABCD22222)(ODOCADDCADAC如圖2222142)2(11)2(nnnnaaaa相應(yīng)相應(yīng)OAC的的面積面積4221211nnnaaADOCS設(shè)邊數(shù)為設(shè)邊數(shù)為62n的正多邊形邊長(zhǎng)為的正多邊形邊長(zhǎng)為 an遞推法遞推法 樂(lè)經(jīng)良 于是于是 的值的值 ( (劉徽計(jì)算到劉徽計(jì)算到1 192邊形面積，得到邊形面積，得到 3.141) )nnn

3、naS232611用用MatMatlablab計(jì)算計(jì)算 m m文件文件function calpi(n)a(1)=1;for i=1:n-1 a(i+1)=sqrt(2-sqrt(4-a(i)2);end S=3*2(n-1)*a(n)命令窗口輸入命令窗口輸入 format long gcalpi(5) 樂(lè)經(jīng)良如何提高精度如何提高精度提高多邊形的邊數(shù)提高多邊形的邊數(shù) 不能完全達(dá)到目的不能完全達(dá)到目的在在Matlab文件中解決文件中解決 符號(hào)運(yùn)算符號(hào)運(yùn)算function calpi1(n)a(1)=sym(1);for i=1:n-1 a(i+1)=sym(sqrt(2-sqrt(4-a(i)2

4、); endS=3*2(n-1)*a(n);vpa(S,60) %最后進(jìn)行數(shù)值計(jì)算，60為數(shù)值計(jì)算過(guò)程中保留的有效數(shù)字 樂(lè)經(jīng)良任務(wù)任務(wù)1 1 德國(guó)人魯?shù)婪蛴靡簧?jì)算圓周率。德國(guó)人魯?shù)婪蛴靡簧?jì)算圓周率。他同樣是用圓的內(nèi)接多邊形逼近圓周，他同樣是用圓的內(nèi)接多邊形逼近圓周，不過(guò)他是從正方形開(kāi)始成倍增加邊數(shù)。不過(guò)他是從正方形開(kāi)始成倍增加邊數(shù)。試推導(dǎo)出他計(jì)算所采用的遞推公式，然試推導(dǎo)出他計(jì)算所采用的遞推公式，然后求后求的近似值到的近似值到1010位和位和2020位位. . 樂(lè)經(jīng)良121)1(5131141nn12)1(53arctan12153nxxxxxnn221422) 1(111nnxxxx利用

5、冪級(jí)數(shù)計(jì)算利用冪級(jí)數(shù)計(jì)算積分導(dǎo)出積分導(dǎo)出取取 x=1 樂(lè)經(jīng)良121)1(51311 41nnnS4偶數(shù)奇數(shù)nnSnnSSnnn1211211111S ( Sn的迭代格式 ) 樂(lè)經(jīng)良用用Matlab計(jì)算計(jì)算創(chuàng)建創(chuàng)建m文件文件 calpi 2.m，內(nèi)容如下內(nèi)容如下：function calpi2(n)S=0;for i=1:n if mod(i,2)=0 S=S-1/(2*i-1); else S=S+1/(2*i-1); endendS=4*S 樂(lè)經(jīng)良 calpi2(1000)ans = 3.14059265383979 calpi2(10000)ans = 3.14149265359003結(jié)果

6、如何結(jié)果如何 ? calpi2(20000)ans = 3.14154265358982精度提高很慢精度提高很慢! ! 樂(lè)經(jīng)良 能不能算得更快一點(diǎn)、更精確一點(diǎn)？能不能算得更快一點(diǎn)、更精確一點(diǎn)？431arctan21arctan12153)21(12)1()21(51)21(31214nnnMachin Machin 公式公式42391arctan51arctan4簡(jiǎn)單公式簡(jiǎn)單公式)31(12)1()31(51)31(313112153nnn 樂(lè)經(jīng)良用用MatMatlablab創(chuàng)建創(chuàng)建m m文件文件 function calpi2_1(n)S=0;for i=1:n if mod(i,2)=0

7、S=S-1/(2*i-1)*(1/(2(2*i-1)+1/(3(2*i-1); else S=S+1/(2*i-1)*(1/(2(2*i-1)+1/(3(2*i-1); endendS=vpa(4*S,30)%觀察30位有效數(shù)字 樂(lè)經(jīng)良 calpi2_1(10)ans = 3.14159257960635063255949717131計(jì)算結(jié)果計(jì)算結(jié)果 calpi2_1(20)ans = 3.14159265358975625659354591335 calpi2_1(50)ans = 3.14159265358979323846264338328 樂(lè)經(jīng)良一個(gè)結(jié)論一個(gè)結(jié)論算法很重要算法很重要計(jì)算

8、機(jī)速度計(jì)算機(jī)速度300300次次/ /秒秒 33.8633.86 10104040 兆兆/ /秒秒從從1950 1950 2000 2000年年10104 4次次 / /秒秒 10 101212次次 / /秒，提高秒，提高1 1億倍億倍 算法算法 ( (解線性方程組解線性方程組 高斯消去法高斯消去法 多重網(wǎng)格法多重網(wǎng)格法) ) 計(jì)算機(jī)速度計(jì)算機(jī)速度運(yùn)算次數(shù)：運(yùn)算次數(shù)：10101818次次 10 106 6次，提高次，提高1 1萬(wàn)億倍萬(wàn)億倍 樂(lè)經(jīng)良任務(wù)任務(wù)2 281arctan51arctan21arctan42) 2) 驗(yàn)證公式驗(yàn)證公式1) 1) 用反正切函數(shù)的冪級(jí)數(shù)展開(kāi)式結(jié)合有關(guān)公式用反正切

9、函數(shù)的冪級(jí)數(shù)展開(kāi)式結(jié)合有關(guān)公式 簡(jiǎn)單公式和簡(jiǎn)單公式和Machin公式所用的項(xiàng)數(shù)公式所用的項(xiàng)數(shù). .求求 ，若要精確到，若要精確到4040位、位、5050位數(shù)字，試比較位數(shù)字，試比較試試用此公式右端作冪級(jí)數(shù)展開(kāi)完成任務(wù)試試用此公式右端作冪級(jí)數(shù)展開(kāi)完成任務(wù)1)1)所需的項(xiàng)數(shù)所需的項(xiàng)數(shù) 樂(lè)經(jīng)良3 3）回憶在微積分中學(xué)習(xí)到的其它級(jí)數(shù)形）回憶在微積分中學(xué)習(xí)到的其它級(jí)數(shù)形式是否可用來(lái)求式是否可用來(lái)求 的值到的值到1010位、位、2020位、位、3030位，相應(yīng)需要級(jí)數(shù)的多少項(xiàng)？位，相應(yīng)需要級(jí)數(shù)的多少項(xiàng)？ 樂(lè)經(jīng)良dxxA102114211)(xxy設(shè)將0,1區(qū)間n等分，取xk=k/n, )(220121nn

10、yyyyynA梯形法)(2)(61224220mmyyyyymASimpson法 利用數(shù)值積分方法利用數(shù)值積分方法yk= 1/ (1+xk2)(41231myyy 樂(lè)經(jīng)良Matlab 計(jì)算計(jì)算創(chuàng)建創(chuàng)建 m 文件文件梯形法梯形法function calpi3(n)x=0:1/n:1;y=1./(1+x.2);S=2*sum(y)-1-0.5;2*S/n 樂(lè)經(jīng)良 calpi3(100)ans = 3.14157598692313 calpi3(500)ans = 3.14159198692313 calpi3(10000)ans = 3.14159265192314 樂(lè)經(jīng)良 用數(shù)值積分計(jì)算用數(shù)值積

11、分計(jì)算 ，分別用梯形法和，分別用梯形法和Simpson法精確到法精確到1010位數(shù)字，用位數(shù)字，用Simpson法精確到法精確到1515位數(shù)字位數(shù)字. .任務(wù)任務(wù)3 3 樂(lè)經(jīng)良針與平行線相交的次數(shù)為針與平行線相交的次數(shù)為n Monte Carlo 法法從從Buffon落針實(shí)驗(yàn)談起：落針實(shí)驗(yàn)談起： 紙上一組平行線距離為紙上一組平行線距離為1 1，將長(zhǎng)度為將長(zhǎng)度為1 1的針多次地扔到的針多次地扔到紙上。若扔針次數(shù)為紙上。若扔針次數(shù)為m, ,而其中而其中Buffon指出：指出： 的數(shù)值與的數(shù)值與 m/n 有關(guān)，他由此有關(guān)，他由此求出求出 的近似值為的近似值為3.1423.142 樂(lè)經(jīng)良設(shè)計(jì)方案設(shè)計(jì)方

12、案 4 4 m/n 計(jì)算機(jī)模擬：產(chǎn)生區(qū)間計(jì)算機(jī)模擬：產(chǎn)生區(qū)間0,10,1上數(shù)目為上數(shù)目為n 的一組的一組 在正方形在正方形 0 x 1, 0 y1上隨機(jī)的投大量的點(diǎn)，那么上隨機(jī)的投大量的點(diǎn)，那么落在四分之一園內(nèi)的點(diǎn)數(shù)落在四分之一園內(nèi)的點(diǎn)數(shù)數(shù)數(shù)m與在正方形內(nèi)的點(diǎn)數(shù)與在正方形內(nèi)的點(diǎn)數(shù)n之比之比 m/n 應(yīng)為這兩部分圖形應(yīng)為這兩部分圖形面積之比面積之比 /4, /4, 故故隨機(jī)數(shù)隨機(jī)數(shù)(x, y)，計(jì)算滿足計(jì)算滿足 x2 + y2 1 的點(diǎn)數(shù)的點(diǎn)數(shù)m 樂(lè)經(jīng)良Matlab計(jì)算計(jì)算創(chuàng)建創(chuàng)建m文件文件function y=calpi4(k)m=0;for n=1:k if rand(1)2+rand(1)

13、2=1 m=m+1; end;end;y=4*m/k;觀察結(jié)果，思考為什么？ 樂(lè)經(jīng)良2) 2) 設(shè)計(jì)方案用計(jì)算機(jī)模擬設(shè)計(jì)方案用計(jì)算機(jī)模擬Buffon實(shí)驗(yàn)實(shí)驗(yàn)任務(wù)任務(wù)4 4看能否求得看能否求得5 5位精確數(shù)字？位精確數(shù)字？1) 1) 用用Monte Carlo 法計(jì)算法計(jì)算 ，除了加大隨機(jī)數(shù)，除了加大隨機(jī)數(shù)，在隨機(jī)數(shù)一定時(shí)可重復(fù)算若干次后求平均值，在隨機(jī)數(shù)一定時(shí)可重復(fù)算若干次后求平均值， 樂(lè)經(jīng)良其他方法其他方法 1/ 的展開(kāi)式的展開(kāi)式 Ramanujan 公式公式 044396) !()263901103()!4(9801221nnnnn 算術(shù)幾何平均值迭代法算術(shù)幾何平均值迭代法nnnnnnnn

14、nnbaMbabbaabalimlim,2,21, 1110020222)(211Mbannnn 樂(lè)經(jīng)良利用積分利用積分dxxn20sin為奇數(shù)nnn,2! !)!1(推導(dǎo)公式推導(dǎo)公式122122543432122nnnn任務(wù)任務(wù)5 5用此公式計(jì)算用此公式計(jì)算 的近似值，效果如何？的近似值，效果如何？ 樂(lè)經(jīng)良有必要計(jì)算那么精確嗎有必要計(jì)算那么精確嗎十位小數(shù)就足以使地球周界準(zhǔn)確到一英十位小數(shù)就足以使地球周界準(zhǔn)確到一英寸以內(nèi)，三十位小數(shù)便能使整個(gè)可見(jiàn)宇宙周寸以內(nèi)，三十位小數(shù)便能使整個(gè)可見(jiàn)宇宙周邊準(zhǔn)確到連最強(qiáng)大的顯微鏡都不能分辨的一個(gè)邊準(zhǔn)確到連最強(qiáng)大的顯微鏡都不能分辨的一個(gè)度量，度量， 三十五位小數(shù)

15、的三十五位小數(shù)的值計(jì)算能把太陽(yáng)系值計(jì)算能把太陽(yáng)系包圍起來(lái)的圓的周長(zhǎng)，誤差還不到質(zhì)子直徑的包圍起來(lái)的圓的周長(zhǎng)，誤差還不到質(zhì)子直徑的百萬(wàn)分之一百萬(wàn)分之一 樂(lè)經(jīng)良計(jì)算計(jì)算 的意義的意義 反映數(shù)學(xué)和計(jì)算技術(shù)發(fā)展的一個(gè)側(cè)面 “歷史上一個(gè)國(guó)家所算得的圓周率的準(zhǔn)確程度，可以作為衡量這個(gè)國(guó)家當(dāng)時(shí)數(shù)學(xué)發(fā)展水平的指標(biāo)?！?3.1415926 3.1415927 113355722(領(lǐng)先世界900余年) 樂(lè)經(jīng)良位數(shù)位數(shù) 100100萬(wàn)萬(wàn) 10.110.1億億 20612061億億 1241112411億億 2.72.7萬(wàn)億萬(wàn)億 1010萬(wàn)億萬(wàn)億年代年代 1973 1989 1999 2002 2010.1 2011

16、.101973 1989 1999 2002 2010.1 2011.10國(guó)家國(guó)家 法國(guó)法國(guó) 美國(guó)美國(guó) 日本日本 日本日本 法國(guó)法國(guó) 日本日本人工計(jì)算：實(shí)驗(yàn)法人工計(jì)算：實(shí)驗(yàn)法 幾何法幾何法 分析法分析法用計(jì)算機(jī)用計(jì)算機(jī): : 始于始于19491949 2035 2035 位位 （Von Neumann）最高記錄：最高記錄：808808位位(1948 )(1948 )家用電腦家用電腦運(yùn)用超級(jí)計(jì)算機(jī)運(yùn)用超級(jí)計(jì)算機(jī) 已算到小數(shù)位已算到小數(shù)位6060萬(wàn)億位萬(wàn)億位 樂(lè)經(jīng)良 引發(fā)新的概念、方法和思想引發(fā)新的概念、方法和思想 , ,產(chǎn)生新的問(wèn)題產(chǎn)生新的問(wèn)題 測(cè)試或檢驗(yàn)超級(jí)計(jì)算機(jī)的各項(xiàng)性能測(cè)試或檢驗(yàn)超級(jí)計(jì)算機(jī)的各項(xiàng)性能( (Super PI) ) IntelIntel公司推出奔騰（公司推出奔騰（Pentium）時(shí)發(fā)現(xiàn)問(wèn)題）時(shí)發(fā)現(xiàn)問(wèn)題 雅虎科技公司的研究員尼古拉斯雅虎科技公司的研究員尼古拉斯斯則斯則(Nicholas Sze)(Nicholas Sze)，采用，采用“云計(jì)算云計(jì)算”技術(shù)，利用技術(shù)，利用10001000臺(tái)電腦同時(shí)計(jì)算，歷時(shí)臺(tái)電腦同時(shí)計(jì)算，歷時(shí)2323天，將圓周率精確天，將圓周率精確到小數(shù)點(diǎn)后到小數(shù)點(diǎn)后2 2千萬(wàn)億位千萬(wàn)億位 樂(lè)經(jīng)良1) 1) 利用學(xué)習(xí)過(guò)的知識(shí)利用學(xué)習(xí)過(guò)的知