pi的计算 数值分析论文
古人计算圆周率，一般是用割圆法。即用圆的内接或外切正多边形来逼近圆的周长。Archimedes用正96边形得到圆周率小数点后3位的精度；刘徽用正3072边形得到5位精度；Ludolph Van Ceulen用正262边形得到了35位精度。这种基于几何的算法计算量大，速度慢，吃力不讨好。随着数学的发展，数学家们在进行数学研究时有意无意地发现了许多计算圆周率的公式。下面挑选一些经典的常用公式加以介绍。除了这些经典公式外，还有很多其它公式和由这些经典公式衍生出来的公式，就不一一列举了。
1、 Machin公式
[这个公式由英国天文学教授John Machin于1706年发现。他利用这个公式计算到了100位的圆周率。Machin公式每计算一项可以得到1.4位的十进制精度。因为它的计算过程中被乘数和被除数都不大于长整数，所以可以很容易地在计算机上编程实现。
Machin.c 源程序
还有很多类似于Machin公式的反正切公式。在所有这些公式中，Machin公式似乎是最快的了。虽然如此，如果要计算更多的位数，比如几千万位，Machin公式就力不从心了。下面介绍的算法，在PC机上计算大约一天时间，就可以得到圆周率的过亿位的精度。这些算法用程序实现起来比较复杂。因为计算过程中涉及两个大数的乘除运算，要用FFT(Fast Fourier Transform)算法。FFT可以将两个大数的乘除运算时间由O(n2)缩短为O(nlog(n))。
2、 Ramanujan公式
1914年，印度数学家Srinivasa Ramanujan在他的论文里发表了一系列共14条圆周率的计算公式，这是其中之一。这个公式每计算一项可以得到8位的十进制精度。1985年Gosper用这个公式计算到了圆周率的17,500,000位。
1989年，David & Gregory Chudnovsky兄弟将Ramanujan公式改良成为：
这个公式被称为Chudnovsky公式，每计算一项可以得到15位的十进制精度。1994年Chudnovsky兄弟利用这个公式计算到了4,044,000,000位。Chudnovsky公式的另一个更方便于计算机编程的形式是：
3、AGM(Arithmetic-Geometric Mean)算法
Gauss-Legendre公式：
这个公式每迭代一次将得到双倍的十进制精度，比如要计算100万位，迭代20次就够了。1999年9月Takahashi和Kanada用这个算法计算到了圆周率的206,158,430,000位，创出新的世界纪录。
4、Borwein四次迭代式：
这个公式由Jonathan Borwein和Peter Borwein于1985年发表，它四次收敛于圆周率。
这个公式简称BBP公式，由David Bailey, Peter Borwein和Simon Plouffe
于1995年共同发表。它打破了传统的圆周率的算法，可以计算圆周率的任意第n位，而不用计算前面的n-1位。这为圆周率的分布式计算提供了可行性。1997年，Fabrice Bellard找到了一个比BBP快40％的公式：
【与π有关的等式】
(π^2)/6 = 1/1^2 + 1/2^2 + 1/3^2 + ······ + 1/n^2 + ······e^(πi) + 1 = 0
e^(-x^2) 在-∞到+∞上的积分是√π
sinx/x 在0到∞上的积分是π/2
[瓦里斯公式] π/2 = lim (n→∞) [ (2n)!! / (2n-1)!! ]^2 / (2n+1)
古人计算圆周率π，一般是用割圆法。即用圆的内接或外切正多边形来逼近圆的周长。阿基米德(Archimedes)用正96边形得到圆周率小数点后3位的精度；刘徽用正3072边形得到5位精度；Ludolph Van Ceulen用正262边形得到了35位精度。这种基于几何的算法计算量大，速度慢。随着数学的发展，数学家们在进行数学研究时有意无意地发现了许多计算圆周率的公式。下面介绍一些经典的公式。
Machin公式
这个公式由英国天文学教授John Machin于1706年发现。他利用这个公式计算到了100位的圆周率。Machin公式每计算一项可以得到1.4位的十进制精度。因为它的计算过程中被乘数和被除数都不大于长整数，所以很容易在计算机上编程实现。
还有很多类似于Machin公式的反正切公式。在所有这些公式中，Machin公式似乎是最快的了。虽然如此，如果要计算更多的位数，比如几千万位，Machin 公式就力不从心了。下面介绍的算法，在PC机上计算大约一天时间，就可以得到圆周率的过亿位的精度。这些算法用程序实现起来比较复杂。因为计算过程中涉及两个大数的乘除运算，要用快速傅利叶变换（FFT，Fast Fourier Transform）算法。FFT可以将两个大数的乘除运算时间由O(n2)缩短为O(nlog(n))。
拉马努安(Ramanujan)公式
1914年，印度数学家Srinivasa Ramanujan在他的论文里发表了一系列共14条圆周率的计算公式，这是其中之一。这个公式每计算一项可以得到8位的十进制精度。1985年Gosper用这个公式计算到了圆周率的17,500,000位。
1989年，David & Gregory Chudnovsky兄弟将Ramanujan公式改良成为：
这个公式被称为Chudnovsky公式，每计算一项可以得到15位的十进制精度。1994年Chudnovsky兄弟利用这个公式计算到了4,044,000,000位。Chudnovsky 公式的另一个更方便于计算机编程的形式是：
AGM(Arithmetic-Geometric Mean)算法
高斯-勒让德（Gauss-Legendre）公式：
初值：
重复计算：
最后计算：
这个公式每迭代一次将得到双倍的十进制精度，比如要计算100万位，迭代20次就够了。1999年9月Takahashi和Kanada用这个算法计算到了圆周率的206,158,430,000位，创出新的世界纪录。
Borwein四次迭代式：
初值：
重复计算：
最后计算：
这个公式由Jonathan Borwein和Peter Borwein于1985年发表，它四次收敛于圆周率。
Bailey-Borwein-Plouffe算法
这个公式简称BBP公式，由David Bailey, Peter Borwein和Simon Plouffe 于1995年共同发表。它打破了传统的圆周率的算法，可以计算圆周率的任意第n位，而不用计算前面的n-1位。这为圆周率的分布式计算提供了可行性。1997年，Fabrice Bellard找到了一个比BBP快40％的公式：
圆周率π小数点后1000位：π=3.
圆周率，一般以π来表示，是一个在数学及物理学普遍存在的数学常数。它定义为圆形之周长与直径之比。它也等于圆形之面积与半径平方之比。是精确计算圆周长、圆面积、球体积等几何形状的关键值。
在分析学上，π可以严格地定义为满足sin(x) = 0的最小正实数x，这里的sin是正弦函数（采用分析学的定义）。
[编辑本段]
【圆周率的历史】
π=Pài(π=Pi)
古希腊欧几里德《几何原本》（约公元前3世纪初）中提到圆周率是常数，中国古算书《周髀算经》（约公元前2世纪）中有“径一而周三”的记载，也认为圆周率是常数。历史上曾采用过圆周率的多种近似值，早期大都是通过实验而得到的结果，如古埃及纸草书（约公元前1700）中取pi=（4/3）^4≒3.1604 。第一个用科学方法寻求圆周率数值的人是阿基米德，他在《圆的度量》（公元前3世纪）中用圆内接和外切正多边形的周长确定圆周长的上下界，从正六边形开始，逐次加倍计算到正96边形，得到(3+(10/71))<π<(3+(1/7)) ，开创了圆周率计算的几何方法（亦称古典方法，或阿基米德方法），得出精确到小数点后两位的π值。
中国数学家刘徽在注释《九章算术》（263年）时只用圆内接正多边形就求得π的近似值，也得出精确到两位小数的π值，他的方法被后人称为割圆术。他用割圆术一直算到圆内接正192边形。
南北朝时代数学家祖冲之进一步得出精确到小数点后7位的π值（约5世纪下半叶），给出不足近似值3.1415926和过剩近似值3.1415927，还得到两个近似分数值，密率355/113和约率22／7。其中的密率在西方直到1573才由德国人奥托得到，1625年发表于荷兰工程师安托尼斯的著作中，欧洲称之为安托尼斯率。
阿拉伯数学家卡西在15世纪初求得圆周率17位精确小数值，打破祖冲之保持近千年的纪录。
德国数学家柯伦于1596年将π值算到20位小数值，后投入毕生精力，于1610年算到小数后35位数，该数值被用他的名字称为鲁道夫数。
无穷乘积式、无穷连分数、无穷级数等各种π值表达式纷纷出现，π值计算精度也迅速增加。1706年英国数学家梅钦计算π值突破100位小数大关。1873 年另一位英国数学家尚可斯将π值计算到小数点后707位，可惜他的结果从528位起是错的。到1948年英国的弗格森和美国的伦奇共同发表了π的808位小数值，成为人工计算圆周率值的最高纪录。
电子计算机的出现使π值计算有了突飞猛进的发展。1949年美国马里兰州阿伯丁的军队弹道研究实验室首次用计算机（ENIAC）计算π值，一下子就算到2037位小数，突破了千位数。1989年美国哥伦比亚大学研究人员用克雷－2型和IBM－VF型巨型电子计算机计算出π值小数点后4.8亿位数，后又继续算到小数点后10.1亿位数，创下最新的纪录。至今，最新纪录是小数点后25769.8037亿位。
[编辑本段]
【圆周率的计算】
余古今中外，许多人致力于圆周率的研究与计算。为了计算出圆周率的越来越好的近似值，一代代的数学家为这个神秘的数贡献了无数的时间与心血。
十九世纪前，圆周率的计算进展相当缓慢，十九世纪后，计算圆周率的世界纪录频频创新。整个十九世纪，可以说是圆周率的手工计算量最大的世纪。
进入二十世纪，随着计算机的发明，圆周率的计算有了突飞猛进。借助于超级计算机，人们已经得到了圆周率的2061亿位精度。
历史上最马拉松式的计算，其一是德国的Ludolph Van Ceulen，他几乎耗尽了一生的时间，计算到圆的内接正262边形，于1609年得到了圆周率的35
位精度值，以至于圆周率在德国被称为Ludolph数；其二是英国的威廉·山克斯，他耗费了15年的光阴，在1874年算出了圆周率的小数点后707位，并将其刻在了墓碑上作为一生的荣誉。可惜，后人发现，他从第528位开始就算错了。
把圆周率的数值算得这么精确，实际意义并不大。现代科技领域使用的圆周率值，有十几位已经足够了。如果用鲁道夫算出的35位精度的圆周率值，来计算一个能把太阳系包起来的一个圆的周长，误差还不到质子直径的百万分之一。以前的人计算圆周率，是要探究圆周率是否循环小数。自从1761年兰伯特证明了圆周率是无理数，1882年林德曼证明了圆周率是超越数后，圆周率的神秘面纱就被揭开了。
现在的人计算圆周率, 多数是为了验证计算机的计算能力，还有，就是为了兴趣。
[编辑本段]
【圆周率的计算方法】
古人计算圆周率，一般是用割圆法。即用圆的内接或外切正多边形来逼近圆的周长。阿基米德用正96边形得到圆周率小数点后3位的精度；刘徽用正3072边形得到5位精度；鲁道夫用正262边形得到了35位精度。这种基于几何的算法计算量大，速度慢，吃力不讨好。随着数学的发展，数学家们在进行数学研究时有意无意地发现了许多计算圆周率的公式。下面挑选一些经典的常用公式加以介绍。除了这些经典公式外，还有很多其它公式和由这些经典公式衍生出来的公式，就不一一列举了。
1、马青公式
π=16arctan1/5-4arctan1/239
这个公式由英国天文学教授约翰·马青于1706年发现。他利用这个公式计算到了100位的圆周率。马青公式每计算一项可以得到1.4位的十进制精度。因为它的计算过程中被乘数和被除数都不大于长整数，所以可以很容易地在计算机上编程实现。
还有很多类似于马青公式的反正切公式。在所有这些公式中，马青公式似乎是最快的了。虽然如此，如果要计算更多的位数，比如几千万位，马青公式就力不从心了。
2、拉马努金公式
1914年，印度天才数学家拉马努金在他的论文里发表了一系列共14条圆周率的计算公式。这个公式每计算一项可以得到8位的十进制精度。1985年Gosper 用这个公式计算到了圆周率的17,500,000位。
1989年，大卫·丘德诺夫斯基和格雷高里·丘德诺夫斯基兄弟将拉马努金公式改良，这个公式被称为丘德诺夫斯基公式，每计算一项可以得到15位的十进制精度。1994年丘德诺夫斯基兄弟利用这个公式计算到了4,044,000,000位。丘德诺夫斯基公式的另一个更方便于计算机编程的形式是：
3、AGM（Arithmetic-Geometric Mean）算法
高斯-勒让德公式：
这个公式每迭代一次将得到双倍的十进制精度，比如要计算100万位，迭代20次就够了。1999年9月，日本的高桥大介和金田康正用这个算法计算到了圆周率的206,158,430,000位，创出新的世界纪录。
4、波尔文四次迭代式：
这个公式由乔纳森·波尔文和彼得·波尔文于1985年发表的。
5、bailey-borwein-plouffe算法
这个公式简称BBP公式，由David Bailey, Peter Borwein和Simon Plouffe 于1995年共同发表。它打破了传统的圆周率的算法，可以计算圆周率的任意第n位，而不用计算前面的n-1位。这为圆周率的分布式计算提供了可行性。
6、丘德诺夫斯基公式[1]
这是由丘德诺夫斯基兄弟发现的，十分适合计算机编程，是目前计算机使用较快的一个公式。以下是这个公式的一个简化版本：
丘德诺夫斯基公式
[编辑本段]
【圆周率的计算历史】
时间纪录创造者小数点后位数所用方法
前2000 古埃及人 0
前1200 中国 0
前500 《圣经》0（周三径一）
前250 阿基米德 3
263 刘徽 5 古典割圆术
480 祖冲之 7
1429 Al-Kashi 14
1593 Romanus 15
1596 鲁道夫 20 古典割圆术
1609 鲁道夫 35
1699 夏普 71 夏普无穷级数
1706 马青 100 马青公式
1719 （法）德·拉尼 127（112位正确）夏普无穷级数
1794（奥地利）乔治·威加 140 欧拉公式
1824 （英）威廉·卢瑟福 208（152位正确）勒让德公式
1844 Strassnitzky & Dase 200
1847 Clausen 248
1853 Lehmann 261
1853 Rutherford 440
1874 威廉·山克斯 707(527位正确)
20世纪后
年月纪录创造者所用机器小数点后位数
1946 （英）弗格森 620
1947 1 （英）弗格森 710
1947 9 Ferguson & Wrench 808
1949 Smith & Wrench 1,120
1949 Reitwiesner et al ENIAC 2,037
1954 Nicholson & Jeenel NORC 3,092
1957 Felton Pegasus 7,480
1958 1 Genuys IBM704 10,000
1958 5 Felton Pegasus 10,021
1959 Guilloud IBM 704 16,167
1961 Shanks & Wrench IBM 7090 100,265
1966 Guilloud & Filliatre IBM 7030 250,000
1967 Guilloud & Dichampt CDC 6600 500,000
1973 Guilloud & Bouyer CDC 7600 1,001,250
1981 Miyoshi & Kanada FACOM M-200 2,000,036
1982 Guilloud 2,000,050
1982 Tamura MELCOM 900II 2,097,144
1982 Tamura & Kanada HITACHI M-280H 4,194,288
1982 Tamura & Kanada HITACHI M-280H 8,388,576
1983 Kanada, Yoshino & Tamura HITACHI M-280H 16,777,206 1985 10 Gosper Symbolics 3670 17,526,200
1986 1 Bailey CRAY-2 29,360,111
1986 9 Kanada & Tamura HITACHI S-810/20 33,554,414
1986 10 Kanada & Tamura HITACHI S-810/20 67,108,839
1987 1 Kanada, Tamura & Kubo et al NEC SX-2 134,217,700 1988 1 Kanada & Tamura HITACHI S-820/80 201,326,551
1989 5 Chudnovskys CRAY-2 & IBM-3090/VF 480,000,000
1989 6 （美）丘德诺夫斯基兄弟 IBM 3090 525,229,270
1989 7 Kanada & Tamura HITACHI S-820/80 536,870,898
1989 8 （美）丘德诺夫斯基兄弟 IBM 3090 1,011,196,691
1989 11 Kanada & Tamura HITACHI S-820/80 1,073,741,799 1991 8 （美）丘德诺夫斯基兄弟 2,260,000,000
1994 5 （美）丘德诺夫斯基兄弟 4,044,000,000
1995 8 Takahashi & Kanada HITACHI S-3800/480 4,294,967,286 1995 10 Takahashi & Kanada 6,442,450,938
1997 7 Takahashi & Kanada 51,539,600,000
1999 4 Takahashi & Kanada 68,719,470,000
1999 9 Takahashi & Kanada HITACHI SR8000 206,158,430,000 2002 金田康正团队 1,241,100,000,000
2009 日本筑波大学 2,576,980,370,000
[编辑本段]
【圆周率的最新计算纪录】
1、新世界纪录
圆周率的最新计算纪录由日本筑波大学所创造。他们于2009年算出π值2,576,980,370,000 位小数，这一结果打破了由日本人金田康正的队伍于2002年创造的1,241,100,000,000位小数的世界纪录。
2、个人计算圆周率的世界纪录
１１月２０日，在位于陕西杨凌的西北农林科技大学，生命科学学院研究生吕超结束背诵圆周率之后，戴上了象征成功的花环。当日，吕超同学不间断、无差错背诵圆周率至小数点后６７８９０位，此前，背诵圆周率的吉尼斯世界纪录为无差错背诵小数点后４２１９５位。整个过程用时２４小时０４分。（新华社报道）
[编辑本段]
【一些数字序列在π小数点后出现的位置】
数字序列出现的位置
01234567891 26,852,899,245 41,952,536,161 99,972,955,571
102,081,851,717 171,257,652,369
01234567890 53,217,681,704 148,425,641,592
432109876543 149,589,314,822
543210987654 197,954,994,289
98765432109 123,040,860,473 133,601,569,485 150,339,161,883 183,859,550,237
09876543210 42,321,758,803 57,402,068,394 83,358,197,954
10987654321 89,634,825,550 137,803,268,208 152,752,201,245
27182818284 45,111,908,393
[编辑本段]
【PC机上的计算】
1、PiFast
目前PC机上流行的最快的圆周率计算程序是PiFast。它除了计算圆周率，还可以计算e和sqrt(2)。PiFast可以利用磁盘缓存，突破物理内存的限制进行
超高精度的计算，最高计算位数可达240亿位，并提供基于Fabrice Bellard 公式的验算功能。
2、PC机上的最高计算记录
最高记录：12,884,901,372位
时间：2000年10月10日
记录创造者：Shigeru Kondo
所用程序：PiFast ver3.3
机器配置：Pentium III 1G, 1792M RAM，WindowsNT4.0，
40GBx2(IDE,FastTrak66)
计算时间：1,884,375秒 (21.809895833333333333333333333333天) 验算时间：29小时
【C++编译器中的运算程序】
微机ＷindowsXP中Dev-cpp中的运算程序（30000位）（C++）
#include
#include
#include
#define N 30015
using namespace std;
void mult (int *a,int b,int *s)
{
for (int i=N,c=0;i>=0;i--)
{
int y=(*(a+i))*b+c;
c=y/10;
*(s+i)=y%10;
}
}
void divi (int *a,int b,int *s)
{
for (int i=0,c=0;i<=N;i++)
{
int y=(*(a+i))+c*10;
c=y%b;
*(s+i)=y/b;
}
}
void incr(int *a,int *b,int *s)
{
for (int i=N,c=0;i>=0;i--)
{
int y=(*(a+i))+(*(b+i))+c;
c=y/10;
*(s+i)=y%10;
}
}
bool eqs(int *a,int *b)
{
int i=0;
while (((*(a+i))==(*(b+i)))&&(i<=N)) i++;
return i>N;
}
int _tmain(int argc, char *argv[])
{
cout << "正在计算 . . . (0%)";
int lpi[N+1],lls[N+1],lsl[N+1],lp[N+1];
int *pi=lpi,*ls=lls,*sl=lsl,*p=lp;
for (int i=0;i<=N;i++)*(pi+i)=*(ls+i)=*(sl+i)=*(p+i)=0; memset(pi,0,sizeof(pi));
memset(ls,0,sizeof(ls));
memset(sl,0,sizeof(sl));
memset(p,0,sizeof(p));
*pi=*ls=*sl=1;
for (int i=1;true;i++)
{
mult(ls,i,sl);
divi(sl,2*i+1,ls);
incr(pi,ls,p);
if (eqs(pi,p))
{
cout << "\b\b\b\b100%)\n";
break;
}
int *t;
t=p;
p=pi;
pi=t;
//if (i%1000==0) cout << i << " ";
if(i%1000 == 0)
{
/*cout << i/1000 << "% ";
if(i%5000 == 0)
cout << endl;*/
if(i/1000 < 11)
{
cout << "\b\b\b";
} else {
cout << "\b\b\b\b";
}
cout << i/1000 << "%)";
}
}
cout << endl;
cout << "计算完成\n正在保存 . . .\n";
mult(p,2,pi);
ofstream fout("pi.txt");
fout << *pi << ".";
for (int i=1;i <= N - 15;i++)
{
fout << *(pi+i);
if (i%10==0) fout << " ";
if (i%80==0) fout << endl;
}
cout << "保存完成\n";
cout<< "按回车键退出";
cin.peek();
return EXIT_SUCCESS;
}
注：①运行时会有数据弹出，这无关紧要，只为了加快了感觉速度；
注：程序中有语法错误。请高人改正。
[编辑本段]
【背圆周率的口诀】
3 . 1
4 1
5 9 2
6 5 3 5 8 9
7 9 3 2 3
8 4 6 2 6
山巅一寺一壶酒，尔乐苦煞吾，把酒吃，酒杀尔，杀不死，乐尔乐。
4 3 3 8 3 2 7 9
5 0 2 8 8 4 1 9 7 1
6 9 3 9 9 3 7
死珊珊，霸占二妻。救我灵儿吧！不只要救妻，一路救三舅，救三妻。
5 1 0 5 8 2 0 9 7 4 9 4 4 5 9 2 3 0 7
我一拎我爸，二拎舅（其实就是撕我舅耳）三拎妻。
8 1 6 4 0 6 2 8 6 2 0 8 9 9 8 6
不要溜！司令溜，儿不溜！儿拎爸，久久不溜！
（作者华罗庚）
来历：有个教书先生，喜欢喝酒，每次总是给学生留道题，就到私塾的后山上找山上的老和尚喝酒。这天，他给学生留了道题,就是背这个圆周率，然后自
己提壶酒就到山上的庙里去了。圆周率位数这么多，不好背啊，其中有个聪明的学生就想出了一个办法，把圆周率编了个打油诗：山巅一寺一壶酒，尔乐苦煞吾，把酒吃；酒杀尔杀不死，乐尔乐。其实就是3.1415926535897932384626的谐音。先生一回来，学生居然都把这个给背了下来，很是奇怪，一想，就什么都明白了，原来是在讽刺他呀……
中国人用的是谐音记忆法，外国人(母语为英语的)一般用字长记忆法。例：
3. 1 4 1 5 9
Now I, even I, would celebrate
2 6 5
3 5
In rhymes inapt, the great
8 9 7 9
Immortal Syracusan, rivaled nevermore,
3 2 3 8 4
Who in his wondrous lore,
6 2 6
Passed on before,
4 3 3 8
Left men his guidance
3 2 7 9
How to circles mensurate.
【背圆周率小数点后位数多的人】
背诵圆周率最多的人：日本人原口证（于2006年10月3日至4日背诵圆周率小数後第100,000位数，总计背诵时间为16个小时半）
一学生背圆周率至小数点后6万位
截至２０日１４时５６分，西北农林科技大学硕士研究生吕超用２４小时零４分钟，不间断无差错地背诵圆周率至小数点后６７８９０位，从而刷新由一名日本学生于１９９５年创造的无差错背诵圆周率至小数点后４２１９５位的吉尼斯世界纪录。
生于１９８２年１１月的吕超，２００１年由湖北省枣阳市考入西北农林科技大学生命科学2００５年被推荐免试攻读本校的应用化学硕士学位。他有较强的记忆能力，特别擅长背诵和默写数字，通常记忆１００位数字只需１０分钟。吕超从４年前开始背诵圆周率，近１年来加紧准备，目前能够记住的圆周率位数超过９万位。在２０日的背诵中，吕超背诵至小数点后６７８９０位时将“０”背为“５”发生错误，挑战结束。
圆周率是一个无穷小数，到目前为止，专家利用超级电脑已计算圆周率到小数点后约１００万兆位。据介绍，挑战背诵圆周率吉尼斯世界纪录的规则是：必须大声地背出；背诵过程中不能给予帮助或（视觉与听觉方面的）提示，也不能有任何形式的协助；背诵必须连续，两个数字之间的间隔不得超过１５秒；背诵
出错时可以更正，但更正必须是在说出下一个数字之前；任何错误（除非错误被立刻更正）都将使挑战失败。因此，吕超在背诵前进行了全面体检，并由家长签字同意，背诵过程中还使用了尿不湿和葡萄糖、咖啡、巧克力来解决上厕所和进食等生理问题。
东方网11月25日消息：昨日，记者从西北农林科技大学获悉，该校学生吕超于去年11月成功创造的“背诵圆周率”吉尼斯世界新纪录，最近被英国吉尼斯总部正式认可，并于今年10月26日向吕超颁发了吉尼斯世界纪录证书。在背诵圆周率的吉尼斯纪录历史上，第一次留下了中国人的名字。
现年24岁的吕超是西北农林科技大学理学院应用化学专业在读硕士生。2005年11月20日，吕超经过连续24小时04分的艰苦努力，无差错背诵圆周率达到小数点后第67890位，打破了“背诵圆周率”吉尼斯世界纪录。此前，背诵圆周率的吉尼斯世界纪录，为无差错背诵小数点后第42195位，是日本人友寄英哲于1995年创造的。
据了解，吕超于2004年利用各种记忆方法开始准备背诵圆周率。2005年暑假，他每天花费10多个小时对圆周率反复记忆、复习，经过两个多月的准备，能够准确背诵小数点9万位以上，遂决定向“背诵圆周率”世界纪录发起挑战。
2006年1月初，吕超向英国吉尼斯总部寄送了全部申报材料。经过详细审核，2006年10月，吉尼斯总部正式认可吕超的挑战纪录，并向吕超颁发了吉尼斯世界纪录证书。
昨日面对鲜花和来自老师、同学们的掌声，吕超格外激动地说：“这是我们集体的荣誉，收获最大的不是这个成绩，而是创造这个纪录的过程。”
吕超透露，在练习背诵圆周率过程中，他多次想到了放弃，背到第二周的时候开始失眠，背到一个月的时候掉头发。但为了实现目标，最终还是坚持下来。
当问及下一步是否还打算刷新自己保持的纪录时，吕超说：“没必要把这个纪录一次次刷新。我希望有更多人具备这个能力，这是对人类记忆能力的一种挑战。”
3月14日，在英国牛津大学科学历史博物馆礼堂内众多专家和观众面前，为了替英国“癫痫症治疗协会”募集资金，英国肯特郡亨里湾的丹尼尔·塔曼特在5小时之内成功地将圆周率背诵到了小数点后面22514位！据悉，塔曼特是世界上25位拥有这项“惊人绝技”的记忆专家之一！
据报道，现年25岁的塔曼特是在小时候患了癫痫症后，才突然发现自己拥有“记忆数字”的惊人能力的。长大并战胜自己的疾病后，塔曼特成了一名记忆专家，他不仅精通多种语言，还成立了一间“记忆技巧公司”。
塔曼特是欧洲背诵圆周率小数点后数字最多的人，但却并不是世界第一。据称，最厉害的人是一名马来西亚大学生，他曾在15小时内将圆周率背诵到小数点后67053位.
[编辑本段]
【与π有关的等式】
(π^2)/6 = 1/1^2 + 1/2^2 + 1/3^2 + ······ + 1/n^2 + ······e^(πi) + 1 = 0
e^(-x^2) 在-∞到+∞上的积分是√π
sinx/x 在0到∞上的积分是π/2
[瓦里斯公式] π/2 = lim (n→∞) [ (2n)!! / (2n-1)!! ]^2 / (2n+1)
祖冲之和圆周率的计算
所谓“圆周率”是指一个圆的周长与其直径的比值。古今中外，许多
人致力于圆周率的研究与计算。为了计算出圆周率的越来越好的近似值，
一代代的数学家为这个神秘的数贡献了无数的时间与心血。
一、计算圆周率的各种方法
早在我国的三国时代，数学家刘徽就用“割圆术”求出了比较精确的
圆周率。他发现：当圆内接正多边形的边数不断增加后，多边形的周长会
越来越逼近圆周长，而多边形的面积也会越来越逼近圆面积。于是，刘徽
利用正多边形面积和圆面积之间的关系，从正六边形开始，逐步把边数加
倍：正十二边形、正二十四边形，正四十八边形……，一直到正三○七二
边形，算出圆周率等于三点一四一六，将圆周率的精度提高到小数点后第
四位。
祖冲之（公元429－500年），是中国南北朝时期著名的数
学家、天文学家。他在刘徽研究的基础上，进一步地发展，经
过既漫长又烦琐的计算，一直算到圆内接正二四五七六边形，
而得到一个结论：圆周率的值介于三点一四一五九二六和三点
一四一五九二七之间，成为世界上最早把圆周率推算出七位数
字的科学家，直到一千年以后，才有西方的数学家达到和超过
祖冲之的成就。同时，他还找到了圆周率的约率：22∕7、密率：
355∕113。
以前人们计算圆周率，是要探究圆周率是否循环小数。自从1761年
Lambert证明了圆周率是无理数，1882年Lindemann证明了圆周率是超越
数后，圆周率的神秘面纱就被揭开了。现在人们计算圆周率，多是为了验
证计算机的计算能力。
古人计算圆周率，一般是用割圆法。但这种基于几何的算法计算量大，速度慢，吃力不讨好。随着数学的发展，数学家们在进行数学研究时有意无意地发现了许多计算圆周率的公式。我们选取其中的一个公式，用VB 编程来实现这个公式。
英国天文学教授John Machin于1706年发现了一个计算圆周率的公式，称为Machin公式，他利用这个公式计算到了100位的圆周率。
还有很多类似于Machin公式的反正切公式。
以下即为Machin公式：
分析其中的arctgx公式可以知道，这是一个级数公式，而在程序设Array计中则可以用一个“累加器”算法来实现。
用流程图来表现，则在流程图中，必定有判别框，并根据判别条件成立与否分别设置了重复部分操作内容的分支流程。
二、算法的程序实现
为了实现这个算法，则需要编制相应的程序，在程序中除了需要用到赋值语句、输入输出语句、其它计算语句外，还必须用到循环语句。
范例：我使用VB来编写程序实现这个算法。
算法中用到了一条输入语句、两个循环语句、一个输出语句以
及多个赋值语句。
（1）建立窗体和输入、输出、命令按钮组件对象。
（2）编写“Command1”触发的程序代码。
在“Private Sub command1_click()”和“End Sub”之间输入以下
第一行，定义了两个整数类型的数值变量I和n，一个单精度浮点数变量x，以及三个双精度浮点数变量pi、arc1和arc2。其中pi用于表示圆周率的值。
第二行，将text1文本框中的数据转换为整型数值并赋值给整型变量i.。
第三行，将x赋值为“1/5”。
第四行，将arc1赋值为“0”。因为arc1是一个乘加器，所以其初值应该是0。
第五行，表示开始一个循环，循环变量n从1开始，步长为1，依次
取值到I，一共循环I次。
第六行，arc1 = arc1 + (-1) ^ (n - 1) * x ^ (2 * n - 1) / (2 * n - 1)，
这是一个累加器的算法，它将变量arc1的原值加上表达式的值，然后将加法运算的结果重新赋值给变量arc1作为arc1的新值。
Machin公式中的级数代数式，转换成表达式则为：
(-1) ^ (n - 1) * x ^ (2 * n - 1) / (2 * n - 1)
第七行。NEXT，表示循环变量n增加一个步长的值1，然后判断“n<=I”是否成立，如果成立则继续循环，否则不再循环直接执行下一个语句。
第八行，将x赋值为“1/239”。“Loop While I<=n”表示当I<=n成立时继续循环，从第五行“Do”的下面一行继续执行。如果I<=n不成立，即I比n大时，则不再循环，直接执行下一行即第九行的语句，从而结束循环。
第九行，将arc2赋值为“0”。Arc2也是一个乘加器。
第十~十二行，通过循环计算arc2的值。
第十三行，pi = 16 * arc1 - 4 * arc2。通过Machin公式计算圆周率的值。
第十四行，将表示圆周率的变量pi的值赋值给“Label1”组件对象的“Caption”属性，输出圆周率。
（3）运行程序。
将第一个文本框中的“Text1”删除，重新输入“9”；然后单击“Command1”，就能在原来“Label1”的位置上输出pi的值“3．1415926998682”。输入“10”时，输出pi的值“
3.14159269986815”。输入“1000”时，输出pi的值仍为“3.14159269986815”。说明，应该这个程序，当n大于10时，在双精度浮点数的数值范围内，圆周率的精度至少可以达到小数点后的14位。
数值分析小论文 董安
数值分析作业 课题名称代数插值法-拉格朗日插值法班级Y110201 研究生姓名董安 学号S2******* 学科、专业机械制造及其自动化 所在院、系机械工程及自动化学院2011 年12 月26日
代数插值法---拉格朗日插值法 数值分析中的插值法是一种古老的数学方法，它来自生产实践。利用计算机解决工程问题与常规手工计算的差异就在于它特别的计算方法.电机设计中常常需要通过查曲线、表格或通过作图来确定某一参量,如查磁化曲线、查异步电动机饱和系数曲线等.手工设计时,设计者是通过寻找坐标的方法来实现.用计算机来完成上述工作时,采用数值插值法来完成。因此学好数值分析的插值法很重要。 插值法是函数逼近的重要方法之一，有着广泛的应用 。在生产和实验中，函数f(x)或者其表达式不便于计算复杂或者无表达式而只有函数在给定点的函数值(或其导数值) ，此时我们希望建立一个简单的而便于计算的函数 (x)，使其近似的代替f(x)，有很多种插值法,其中以拉格朗日(Lagrange)插值和牛顿(Newton)插值为代表的多项式插值最有特点,常用的插值还有Hermit 插值,分段插值和样条插值.本文着重介绍拉格朗日(Lagrange)插值法。 1.一元函数插值概念 定义 设有m+1个互异的实数1x ，2x ，···，m x 和n+1 个实值函数()0 x j ， ()1 x j ， ···()n x j ，其中n ￡m 。若向量组 k f =（()0k x j ，()1k x j ，···，() k m x j ）T （k=0,1，，n ） 线性无关，则称函数组{()k x j （k=0,1， ，n ）}在点集{i x (i=0,1, ,m)}上线性无关；否 则称为线性相关。 例如，函数组{2+x ，1-x ，x+2 x }在点集{1,2,3,4}上线性无关。 又如，函数组{sin x ，n2x ，sin 3x }在点集{0， 3p ，2 3 p ，p }上线性相关。 给点n+1个互异的实数0x ，1x ，···，n x ，实值函数() f x 在包含0x ，1x ，···，n x 的某个区间[] ,a b 内有定义。设函数组 {()k x j （k=0,1， ，n ）} 是次数不高于n 的多项式组，且在点集{0x ，1x ，···，n x }上线性无关。
数值分析_数值计算小论文
Runge-Kutta 法的历史发展与应用 摘要Runge-Kutta 法是极其重要的常微分方程数值解法，本文仅就其起源及发展脉络加以简要研究。对Runge 、Heun 以及Kutta 等人的贡献做出适当评述，指出Runge-Kutta 方法起源于Euler 折线法。同时对Runge-Kutta 法的应用做简要研究。 关键词 Euler 折线法 标准四阶Runge-Kutta 法 应用 一、发展历史[1] 1.1 Euler 折线法 在微分方程研究之初，瑞士数学家L.Euler(1707.4—1783.9)做出了开创性的工作。他和其他一些数学家在解决力学、物理学问题的过程中创立了微分方程这门学科。在常微分方程方面，Euler 在1743年发表的论文中，用代换kx y e =给出了任意阶常系数线性微分方程的古典解法，最早引入了“通解”和“特解”的概念。 1768年，Euler 在其有关月球运行理论的著作中，创立了广泛用于求初值问题 00 (,), (1.1)() (1.2)y f x y x x X y x a '=<≤??=? 的数值解的方法，次年又把它推广到二阶方程。欧拉的想法如下：我们选择0h >，然后在00x x x h ≤≤+情况下用解函数的切线 0000()()(,)l x y x x f x y =+- 代替解函数。这样对于点 10x x h =+ 就得到 1000(,)y y hf x y =+。 在11(,)x y 重复如上的程序再次计算新的方向就会得到所谓的递推公式： 11, (,),m m m m m m x x h y y hf x y ++=+=+
数值分析之幂法及反幂法C语言程序实例
数值分析之幂法及反幂法C 语言程序实例 1、算法设计方案： ①求1λ、501λ和s λ的值： s λ：s λ表示矩阵的按模最小特征值，为求得s λ直接对待求矩阵A 应用反幂法即可。 1λ、501λ：已知矩阵A 的特征值满足关系 1n λλ<< ，要求1λ、及501λ时，可 按如下方法求解： a ． 对矩阵A 用幂法，求得按模最大的特征值1m λ。 b ． 按平移量1m λ对矩阵A 进行原点平移得矩阵1m B A I λ=+，对矩阵B 用反幂法 求得B 的按模最小特征值2m λ。 c ． 321m m m λλλ=- 则：113min(,)m m λλλ=，13max(,)n m m λλλ=即为所求。 ②求和A 的与数5011 140 k k λλμλ-=+最接近的特征值 ik λ（k=0，1，…39）： 求矩阵A 的特征值中与k μ最接近的特征值的大小，采用原点平移的方法： 先求矩阵 B=A-k μI 对应的按模最小特征值k β，则k β+k μ即为矩阵A 与k μ最接近的特征值。 重复以上过程39次即可求得ik λ（k=0，1，…39）的值。 ③求A 的（谱范数）条件数2cond()A 和行列式det A ： 在（1）中用反幂法求矩阵A 的按模最小特征值时，要用到Doolittle 分解方法，在Doolittle 分解完成后得到的两个矩阵分别为L 和U ，则A 的行列式可由U 阵求出，即：det(A)=det(U)。 求得det(A)不为0，因此A 为非奇异的实对称矩阵，则： max 2()s cond A λλ= ，max λ和s λ分别为模最大特征值与模最小特征值。
数值分析小论文
“数值分析”课程 第一次小论文 郑维珍2015210459 制研15班（精密仪器系）内容：数值分析在你所在研究领域的应用。 要求：1）字数2500以上；2）要有摘要和参考文献；3）截至10.17，网络学堂提交，过期不能提交！ 数值分析在微流控芯片研究领域的应用 摘要： 作者在硕士期间即将参与的课题是微流控芯片的研制。当前，微流控芯片发展十分迅猛，而其中涉及到诸多材料学、电子学、光学、流体力学等领域的问题，加上微纳尺度上的尺寸效应，理论研究和数值计算都显得困难重重。发展该领域的数值计算，成为重中之重。本文从微流体力学、微传热学、微电磁学、微结构力学等分支入手，简要分析一下数值分析方法在该领域的应用。 微流控芯片(Microfluidic Chip)通常又称芯片实验室(Lab-On-a-Chip )，它是20世纪90年代初由瑞士的Manz和Widmer提出的[1-2]，它通过微细加工技术，将微管道、微泵、微阀、微电极、微检测元件等功能元件集成在芯片材料(基片)上，完成整个生化实验室的分析功能，具有减少样品的消耗量、节省反应和分析的时间、高通量和便携性等优点。 通常一个微流控芯片系统都会执行一个到多个微流体功能，如泵、混合、热循环、扩散和分离等，精确地操纵这些流体过程是微流控芯片的关键。因此它的研究不仅需要生命科学、MEMS、材料学、电子学、光学、流体力学等多学科领域的基础理论的支持，还需要很多数学计算。
1）微流体力学计算[3]： 对微管里的流体动力的研究主要包含了以下几个方面：(1)微管内流体的粘滞力的研究；(2)微管内气流液流的传热活动；(3)在绝热或传热的微管内两相流的流动和能量转换。这三方面的研究涵盖了在绝热、传热和多相转换条件下，可压缩和不可压缩流体在规则或不规则的微管内的流动特性研究。 由此，再结合不同的初值条件和边界条件，我们可以得到各种常微分方程或偏微分方程，而求解这些方程，就是需要很多数值分析的知识。例如，文献[4]里就针对特定的初值和边界条件，由软件求解了Navier-Stodes方程： 文献[4]专门有一章节讨论了该方程的离散化和数值求解。 微流体力学主要向两个方面发展：一方面是研究流动非定常稳定特性、分叉解及微尺寸效应下的湍流流动的机理，更为复杂的非定常、多尺度的流动特征，高精度、高分辨率的计算方法和并行算法；另一方面是将宏观流体力学的基本模型，结合微纳效应，直接用于模拟各种实际流动，解决微纳芯片生产制造中提出来的各种问题。 2）微传热方程计算： 常微分、偏微分方程的数值求解应用较为广泛的另一问题就是微流体传热问题。由传热学的相关知识，我们可以达到如下的传热学基本方程： 该方程在二维情况下经过简化和离散，可以得到如教材第三章所讲的“五点差分格式”的方程组，从而采取数值方法求解[5]。 除此之外，微结构芯片在加工和制造过程中也会有很多热学方面的问题，例如文献[6]所反映的注塑成型工艺中，就有大量的类似问题的解决。 3）微电磁学计算： 由于外加电场的作用，电渗流道中会产生焦耳热效应。许多研究者对电渗流道中的焦耳热效应进行了数值模拟研究。新加坡南洋理工大学的G. Y. Tang等在电渗流模型的基础上，考虑了与温度有关的物理系数，在固一液祸合区域内利用
演讲稿数值分析应用实例.doc
非线性方程求根 问题：在相距100m的两座建筑物（高度相等的点）之间悬挂一根电缆，仅允许电缆在中间最多下垂1m，试计算所需电缆的长度。 设空中电缆的曲线（悬链线）方程为 ] , [ , ) ( 50 50 2 - ∈ + = - x e e a y a x a x （1） 由题设知曲线的最低点)) ( , (0 0y与最高点)) ( , (50 50y之间的高度差为1m，所以有 1 2 50 50 + = +- a e e a a a) ( （2） 由上述方程解出a后，电缆长度可用下式计算： ) ( ) (a a a x a x L e e a dx e e dx x y ds L 50 50 50 50 50 2 1- - - - = ? ? ? ? ? ? + = ' + = =? ? ?（3） 相关Matlab命令： 1、描绘函数] , [ , ) ( ) (1500 500 1 2 50 50 ∈ - - + = - a a e e a a y a a 的图形；
2、用fzero 命令求方程在1250=a 附近的根的近似值x ，并计算)(x y 的函数值； 3、编写二分法程序，用二分法求0=)(a y 在],[13001200内的根，误差不超过310-，并给出对分次数； 4、编写Newton 迭代法程序，并求0=)(a y 在],[13001200内的根，误差不超过310-，并给出迭代次数。 5、编写Newton 割线法程序，并求0=)(a y 在],[13001200内的根，误差不超过310-，并给出迭代次数。
线性方程组求解应用实例 问题：投入产出分析 国民经济各个部门之间存在相互依存的关系，每个部门在运转中将其他部门的产品或半成品（称为投入）经过加工变为自己的产品（称为产出），如何根据各部门间的投入产出关系，确定各部门的产出水平，以满足社会需求，是投入产出分析中研究的课题。考虑下面的例子： 设国民经济由农业、制造业和服务业三个部门构成，已知某年它们之间的投入产出关系、外部需求、初始投入等如表1所示（数字表示产值）。 表1 国民经济三个部门间的关系单位：亿元 假定总投入等于总产出，并且每个部门的产出与它的投入成正比，由上表可以确定三个部门的投入产出表：如表2所示。 表2 三个部门的投入产出表
研究生《数值分析》教学大纲
研究生《数值分析》教学大纲 课程名称：数值分析 课程编号：S061005 课程学时：64 学时 课程学分： 4 适用专业：工科硕士生 课程性质：学位课 先修课程：高等数学，线性代数，计算方法，Matlab语言及程序设计 一、课程目的与要求 “数值分析”课是理工科各专业硕士研究生的学位课程。主要介绍用计算机解决数学问题的数值计算方法及其理论。内容新颖，起点较高，并加强了数值试验和程序设计环节。通过本课程的学习，使学生熟练掌握各种常用的数值算法的构造原理和过程分析，提高算法设计和理论分析能力，并且能够根据数学模型，提出相应的数值计算方法编制程序在计算机上算出结果。力求使学生掌握应用数值计算方法解决实际问题的常用技巧。 二、教学内容、重点和难点及学时安排： 第一章? 数值计算与误差分析( 4学时) 介绍数值分析的研究对象与特点，算法分析与误差分析的主要内容。 第一节数值问题与数值方法 第二节数值计算的误差分析 第三节数学软件工具----MATLAB 语言简介 重点：误差分析 第二章? 矩阵分析基础( 10学时) 建立线性空间、赋范线性空间、内积空间的概念，为学习以后各章打好基础。矩阵分解是解决数值代数问题的常用方法，掌握矩阵的三角分解、正交分解、奇异值分解，并能够编写算法程序。 第一节? 矩阵代数基础
第二节? 线性空间 第三节? 赋范线性空间 第四节? 内积空间和内积空间中的正交系 第五节矩阵的三角分解 第六节矩阵的正交分解 第七节矩阵的奇异值分解 难点：内积空间中的正交系。矩阵的正交分解。 重点：范数，施密特(Schmidt) 正交化过程，正交多项式，矩阵的三角分解, 矩阵的正交分解。 第三章? 线性代数方程组的数值方法( 12学时) 了解研究求解线性代数方程组的数值方法分类及直接法的应用范围。高斯消元法是解线性代数方程组的最常用的直接法，也是其它类型直接法的基础。在此方法基础上加以改进，可得选主元的高斯消元法、按比例增减的高斯消元法，其数值稳定性更高。掌握用列主元高斯消元法解线性方程组及计算矩阵的行列式及逆，并且能编写算法程序。掌握矩阵的直接三角分解法：列主元LU 分解，Cholesky分解。了解三对角方程组的追赶法的分解形式及数值稳定性的充分条件。掌握矩阵条件数的定义，并能利用条件数判别方程组是否病态以及对方程组的直接方法的误差进行估计。 迭代解法是求解大型稀疏方程组的常用解法。熟练掌握雅可比迭代法、高斯- 塞德尔迭代法及SOR 方法的计算分量形式、矩阵形式，并能在计算机上编出三种方法的程序用于解决实际问题。了解极小化方法：最速下降法、共轭斜量法。迭代法的收敛性分析是研究解线性代数方程组的迭代法时必须考虑的问题。对于上述常用的迭代法，须掌握其收敛的条件。而对一般的迭代法，掌握其收敛性分析的基本方法和主要结果有助于进一步探究新的迭代法。 第一节求解线性代数方程组的基本定理 第二节高斯消元法及其计算机实现 第三节矩阵分解法求解线性代数方程组 第三节? 误差分析和解的精度改进 第四节? 大型稀疏方程组的迭代法 第五节? 极小化方法 难点：列主元高斯消元法，直接矩阵三角分解。迭代法的收敛性，雅可比迭代法，高斯-塞德尔迭代法，SOR 迭代法。
《数值分析》课程设计报告
《数值分析》课程设计实验报告 龙格—库塔法分析Lorenz 方程 200820302033 胡涛 一、问题叙述 考虑著名的Lorenz 方程 () dx s y x dt dy rx y xz dt dz xy bz dt ?=-???=--???=-?? 其中s ，r ，b 为变化区域内有一定限制的实参数，该方程形式简单，表面上看并无惊人之处，但由该方程揭示出的许多现象，促使“混沌”成为数学研究的崭新领域，在实际应用中也产生了巨大的影响。 二、问题分析 Lorenz 方程实际上是一个四元一阶常微分方程，用解析法精确求解是不可能的，只能用数值计算，最主要的有欧拉法、亚当法和龙格- 库塔法等。为了得到较高精度的，我们采用经典四阶龙格—库塔方法求解该问题。 三、实验程序及注释 （1）算法程序 function [T]=Runge_Kutta(f,x0,y0,h,n) ％定义算法，其中f 为待解方程组， x0是初始自变量，y0是初始函数 值，h 是步长，n 为步数 if nargin<5 n=100; ％如果输入参数个数小于5，则步数 n=100 end r=size(y0);r=r(1); ％返回初始输出矩阵的行列数，并将 值赋给r(1) s=size(x0);s=s(1); ％返回初始输入矩阵的行列数，并 将值赋给s(1) r=r+s; T=zeros(r,n+1); T(:,1)=[y0;x0]; for t=2:n+1 ％以下是具体的求解过程 k1=feval(f,T(1:r-1,t-1)); k2=feval(f,[k1*(h/2)+T(1:r-1,t-1);x0+h/2]); k3=feval(f,[k2*(h/2)+T(1:r-1,t-1);x0+h/2]); k4=feval(f,[k3*h+T(1:r-1,t-1);x0+h]); x0=x0+h; T(:,t)=[T(1:r-1,t-1)+(k1+k2*2+k3*2+k4)*(h/6);x0]; end
数值分析在生活中的应用举例及Matlab实现
Matlab 实验报告 学院：数学与信息科学学院班级：信息班 学号：20135034027 姓名：马永杉
最小二乘法，用MATLAB实现 1.数值实例 下面给定的是郑州最近1个月早晨7：00左右的天气预报所得到的温度，按照数据找出任意次曲线拟合方程和它的图像。下面用MATLAB编程对上述数据进行最小二乘拟合。 2、程序代码 x=[1:1:30]; y=[9,10,11,12,13,14,13,12,11,9,10,11,12,13,14,12,11,10,9,8,7,8,9,11,9 ,7,6,5,3,1]; a1=polyfit(x,y,3) %三次多项式拟合% a2= polyfit(x,y,9) %九次多项式拟合% a3= polyfit(x,y,15) %十五次多项式拟合% b1=polyval(a1,x) b2=polyval(a2,x) b3=polyval(a3,x) r1= sum((y-b1).^2) %三次多项式误差平方和% r2= sum((y-b2).^2) %九次次多项式误差平方和% r3= sum((y-b3).^2) %十五次多项式误差平方和% plot(x,y,'*') %用*画出x,y图像% hold on plot(x,b1, 'r') %用红色线画出x,b1图像% hold on plot(x,b2, 'g') %用绿色线画出x,b2图像% hold on plot(x,b3, 'b:o') %用蓝色o线画出x,b3图像% 2.流程图
4．数值结果分析 不同次数多项式拟合误差平方和为： r1=67.6659 r2=20.1060 r3=3.7952 r1、r2、r3分别表示三次、九次、十五次多项式误差平方和。 5、拟合曲线如下图
数值分析论文 (8)
牛顿迭代法及其应用 [摘要]本文研究应用泰勒展开式构造出牛顿迭代法，论证了它的局部收敛性和收敛阶。分别讨论了单根情形和重根情形，给出了实例应用。最后给出了离散牛顿法的具体做法。 [关键词] 关键词：泰勒展开式，牛顿迭代法及其收敛性，重根，离散牛顿法。 1.牛顿法及其收敛性 求方程f（x）=0的根，如果已知它的一个近似,可利用Taylor展开式求出f(x)在附近的线性近似，即 ，ξ在x与之间 忽略余项，则得方程的近似 右端为x的线性方程，若，则解,记作,它可作为的解的新近似，即 (2.4.1) 称为解方程的牛顿法.在几何上求方程的解，即求曲线y=f(x)与x轴交点.若已知的一个近似,通过点(，f())作曲线y=f(x)的切线，它与 x轴交点为,作为的新近似，如图1所示
图1 关于牛顿法收敛性有以下的局部收敛定理. 定理1设是f(x)=0的一个根，f(x)在附近二阶导数连续，且，则牛顿法(2.4.1)具有二阶收敛，且 (2.4.2) 证明由式(2.4.1)知迭代函数,， ,而，由定理可知，牛顿迭代(2.4.1)具有二阶收敛，由式可得到式(2.4.2).证毕. 定理表明牛顿法收敛很快，但在附近时才能保证迭代序列收敛.有关牛顿法半局部收敛性与全局收敛定理.此处不再讨论. 例1用牛顿法求方程的根. ,牛顿迭代为 取即为根的近似，它表明牛顿法收敛很快.
例2设＞0，求平方根的过程可化为解方程.若用牛顿法求解，由式(2.4.1)得 (2.4.3) 这是在计算机上作开方运算的一个实际有效的方法，它每步迭代只做一次除法和一次加法再做一次移位即可，计算量少，又收敛很快，对牛顿法我们已证明了它 的局部收敛性，对式(2.4.3)可证明对任何迭代法都是收敛的，因为当 时有 即,而对任意，也可验证,即从k=1开始，且 所以｛｝从k=1起是一个单调递减有下界的序列，｛｝有极限.在式(2.4.3) 中令k→∞可得，这就说明了只要，迭代(2.4.3)总收敛到,且是二阶收敛. 在例2.4的迭代法(3)中，用式(2.4.3)求只迭代3次就得到 =1.732 051,具有7位有效数字. 求非线性方程f(x)=0的根x*,几何上就是求曲线y=f(x)与x轴交点x*，若已知曲线上一点过此点作它的切线。方程为 此切线与x轴交点记作，它就是(2,4,1)给出的牛顿迭代法，由图2-3 看到牛顿法求根就是用切线近似曲线，切线与x轴交点xk+1作为方程f(x)=0 根x*的新近似。 根据定理2.3可以证明牛顿法是二阶收敛的，这就是定理4.1给出的结果，牛顿法由于收敛快，它是方程求根最常用和最重要的方法，在计算机上用牛顿法解方程的计算步骤： 算法如下：（牛顿法） 步0： 给初始近似,计算精度最大迭代步数N,0→k.
数值分析学习心得体会.doc
数值分析学习感想 一个学期的数值分析，在老师的带领下，让我对这门课程有了深刻的理解和感悟。这门 课程是一个十分重视算法和原理的学科，同时它能够将人的思维引入数学思考的模式，在处 理问题的时候，可以合理适当的提出方案和假设。他的内容贴近实际，像数值分析，数值微 分，求解线性方程组的解等，使数学理论更加有实际意义。 数值分析在给我们的知识上，有很大一部分都对我有很大的帮助，让我的生活和学习有 了更加方便以及科学的方法。像第一章就讲的误差，在现实生活中，也许没有太过于注意误 差，所以对误差的看法有些轻视，但在学习了这一章之后，在老师的讲解下，了解到这些误 差看似小，实则影响很大，更如后面所讲的余项，那些差别总是让人很容易就出错，也许在 别的地方没有什么，但是在数学领域，一个小的误差，就很容易有不好的后果，而学习了数 值分析的内容，很容易就可以将误差锁定在一个很小的范围内，在这一范围内再逼近，得出 的近似值要准确的多，而在最开始的计算中，误差越小，对后面的影响越小，这无疑是好的。 数值分析不只在知识上传授了我很多，在思想上也对我有很大的影响，他给了我很多数 学思想，很多思考的角度，在看待问题的方面上，多方位的去思考，并从别的例子上举一反三。像其中所讲的插值法，在先学习了拉格朗日插值法后，对其理解透彻，了解了其中 的原理和思想，再学习之后的牛顿插值以及三次样条插值等等，都很容易的融会贯通，很容 易的就理解了其中所想，他们的中心思想并没有多大的变化，但是使用的方式却是不同的， 这不仅可以学习到其中心内容，还可以去学习他们的思考方式，每个不同的思考方式带来的 都是不同的算法。而在看待问题上，不同的思考方式总是可以快速的全方位的去看透彻问题， 从而知道如何去解决。 在不断的学习中，知识在不断的获取，能力在不断的提升，同时在老师的不懈讲解下， 我逐渐的发现数值分析所涵盖的知识面特别的广泛，而我所需要学习的地方也更加的多，自 己的不足也在不断的体现，我知道这只是我刚刚接触到了数学的那一角，在以后我还会接触 到更多，而这求知的欲望也在不停的驱赶我，学习的越多，对今后的生活才会有更大的帮助。 计算132 2013014923 张霖篇二：数值分析学习报告 数值分析学习心得报告 班级：11级软工一班 姓名： * * * 学号: 20117610*** 指导老师：* * * 学习数值分析的心得体会 无意中的一次选择，让我接触了数值分析。 作为这学期的选修课，我从内心深处来讲，数值分析真的有点难。感觉它是在高等数学 和线性代数的基础上，又加深了探讨。虽然这节课很难，我学的不是很好，但我依然对它比 较感兴趣。下面就具体说说我的学习体会，让那些感兴趣的同学有个参考。 学习数值分析，我们首先得知道一个软件——matlab。matrix laboratory，即矩阵实验 室，是math work公司推出的一套高效率的数值计算和可视化软件。它是当今科学界最具影 响力、也是最具活力的软件，它起源于矩阵运算，并高速发展成计算机语言。它的优点是强 大的科学运算、灵活的程序设计流程、高质量的图形可视化与界面、便捷的与其他程序和语 言接口。 根据上网搜集到的资料，你就会发现matlab有许多优点： 首先，编程简单使用方便。到目前为止，我已经学过c语言，机器语言，java语言，这
数值分析课程设计学生题目
《数值分析》课程设计
本课程设计的内容为：每个小组的同学均应完成以下五个案例； 目标：能将数值分析课程中所学的算法知识熟练应用于实际问题中。 案例1 土木工程和环境工程师在设计一条排水渠道时必须考虑渠道的各种参数（如宽度，深度，渠道内壁光滑度）及水流速度、流量、水深等物理量之间的关系。 假设修一条横断面为矩形的水渠，其宽度为B ，假定水流是定常的，也就是说水流速度不随时间而变化。 根据质量守恒定律可以得到 Q=UBH （1.1） 其中Q 是水的流量（s m /3 ），U 是流速（s m /），H 是水的深度（m ）。 在水工学中应用的有关流速的公式是 3 /23 /22/1)2()(1H B BH S n U += （1.2） 这里n 是Manning 粗糙系数，它是一个与水渠内壁材料的光滑性有关的无量纲量；S 是水渠 的斜度系数，也是一个无量纲量，它代表水渠底每米内的落差。 把（1.2）代入（1.1）就得到 3 /23 /52/1)2()(1H B BH S n U += （1.3） 为了不同的工业目的（比如说要把污染物稀释到一定的浓度以下，或者为某工厂输入一定量 的水），需要指定流量Q 和B ，求出水的深度。这样，就需要求解 0) 2()(1)(3 /23 /52/1=-+=Q H B BH S n H f （1.4） 一个具体的案例是 s m Q S n m B /5 ,0002.0 ,03.0 ,203==== 求出渠道中水的深度H 。 所涉及的知识——非线性方程解法。 案例2 在化学工程中常常研究在一个封闭系统中同时进行的两种可逆反应 C D A C B A ?+?+2 其中A ，B ，C 和D 代表不同的物质。反应达到平衡是有如下的平衡关系： d a c b a c C C C k C C C k == 22 1 , 其中2 24 1107.3 ,104--?=?=k k 称为平衡常数，),,,(d c b a n C n =代表平衡状态时该物质的浓度。假定反应开始时各种物质的浓度为：
数值分析论文
插值方法总结 摘 要：本文是对学过的插值方法进行了总结使我们更清楚的知道那一种方法适合那一种型。 关键词：插值；函数；多项式；余项 （一）Lagrange 插值 1．Lagrange 插值基函数 n+1个n 次多项式 ∏≠=--= n k j j j k j k x x x x x l 0)( n k ,,1,0 = 称为Lagrange 插值基函数 2．Lagrange 插值多项式 设给定n+1个互异点))(,(k k x f x ，n k ,,1,0 =，j i x x ≠，j i ≠，满足插值条件 )()(k k n x f x L =，n k ,,1,0 = 的n 次多项式 ∏∏ ∏=≠==--==n k n k j j j k j k k n k k n x x x x x f x l x f x L 0 00 ))(()()()( 为Lagrange 插值多项式，称 ∏=+-+=-=n j j x n n x x n f x L x f x E 0)1()()!1()()()()(ξ 为插值余项，其中),()(b a x x ∈=ξξ （二）Newton 插值 1．差商的定义 )(x f 关于i x 的零阶差商 )(][i i x f x f = )(x f 关于i x ，j x 的一阶差商 i j i j j i x x x f x f x x f --= ][][],[ 依次类推，)(x f 关于i x ，1+i x ，……，k i x +的k 阶差商
i k i k i i k i i k i i i x x x x f x x f x x x f --= +-+++++] ,,[],,[],,,[111 2．Newton 插值多项式 设给定的n+1个互异点))(,(k k x f x ，n k ,,1,0 =，j i x x ≠，j i ≠， 称满足条件 )()(k k n x f x N =，n k ,,1,0 = 的n 次多项式 )()](,,,[)](,[][)(10100100---++-+=n n n x x x x x x x f x x x x f x f x N 为Newton 插值多项式，称 ],[,)(],,,[)()()(0 10b a x x x x x x f x N x f x E n j j n n ∈-=-=∏= 为插值余项。 （三）Hermite 插值 设],[)(1b a C x f ∈，已知互异点0x ，1x ，…，],[b a x n ∈及所对应的函数值为 0f ，1f ，…，n f ，导数值为'0f ，' 1f ，…，' n f ，则满足条件 n i f x H f x H i i n i i n ,,1,0,)(,)(' '1212 ===++ 的12+n 次Hermite 插值多项式为 )()()(0 '12x f x f x H j n j j j n j i n βα∏∏=++= 其中 )())((,)]()(21[)(2 2'x l x x x l x l x x x j j j j j j j j ---=βα 称为Hermite 插值基函数，)(x l j 是Lagrange 插值基函数，若],[22b a C f n +∈，插值误差为 220) 22(12)()()! 22() ()()(n x n n x x x x n f x H x f --+= -++ ξ，),()(b a x x ∈=ξξ （四）分段插值 设在区间],[b a 上给定n+1个插值节点 b x x x a n =<<<= 10 和相应的函数值0y ，1y ，…，n y ，求作一个插值函数)(x ?，具有性质
数值分析课程课程设计汇总
课 程 设 计 我再也回不到大二了， 大学是那么短暂 设计题目 数值分析 学生姓名 李飞吾 学 号 x x x x x x x x 专业班级 信息计x x x x x 班 指导教师 设 计 题 目 共15题如下 成绩
数值分析课程设计 1．1 水手、猴子和椰子问题：五个水手带了一只猴子来到南太平洋的一个荒岛上，发现那里有一大堆椰子。由于旅途的颠簸，大家都很疲惫，很快就入睡了。第一个水手醒来后，把椰子平分成五堆，将多余的一只给了猴子，他私藏了一堆后便又去睡了。第二、第三、第四、第五个水手也陆续起来，和第一个水手一样，把椰子分成五堆，恰多一只猴子，私藏一堆，再去入睡，天亮以后，大家把余下的椰子重新等分成五堆，每人分一堆，正好余一只再给猴子，试问原先共有几只椰子？（15621） 试分析椰子数目的变化规律，利用逆向递推的方法求解这一问题 解：算法分析：解该问题主要使用递推算法，关于椰子数目的变化规律可以设起初的椰子数为0p ，第一至五次猴子在夜里藏椰子后，椰子的数目分别为01234,,,,p p p p p 再设最后每个人分得x 个椰子,由题: 14 (1)5 k k p p +=- （k=0,1,2,3,4）51(1)5 x p =- 所以551p x =+,11k k p p +=+利用逆向递推方法求解 15 1,4 k k p p +=+ （k=0,1,2,3,4） MATLAB 代码： n=input('n= '); n= 15621 for x=1:n p=5*x+1; for k=1:5 p=5*p/4+1; end if p==fix(p), break end end disp([x,p]) 1．2 设，1 5n n x I dx x =+? （1）从0I 尽可能精确的近似值出发，利用递推公式： 11 5(1,2,20)n n I I n n -=-+= 计算机从1I 到20I 的近似值； （2）从30I 较粗糙的估计值出发，用递推公式：
数值分析小论文论文
对于牛顿型方法的改进 对于函数f(x),假定已给出极小点* x 的一个较好的近似点0x ，则在0x 处将f(x)泰勒展开到二次项，得二次函数()x φ。按极值条件'()0x φ=得()x φ的极小点，用它作为*x 的第一个近似点。然后再在1x 处进行泰勒展开，并求得第二个近似点2x 。如此迭代下去，得到一维情况下的牛顿迭代公式'k 1''k ()() k k f x x x f x +=- (k=0,1,2,…) 对于多元函数f(x),设k x 为f(x)极小点*x 的一个近似值，在k x 处将f(x)进行泰勒展开，保留到二次项得21()()()()()()()()2T T k k k k k k f x x f x f x x x x x f x x x ?≈=+?-+ -?-， 式中 2()k f x ?—f(x)在k x 处的海赛矩阵。 设1k x +为()x ?的极小点，它作为f(x)极小点*x 的下一个近似点，根据极值必要条件 1()0k x ?+?=即21()()()k k k k f x f x x x +?+?-得1 21()()k k k k x x f x f x -+??=-???? （k=0,1,2，…） 上式为多元函数求极值的牛顿法迭代公式。 对于二次函数，f(x)的上述泰勒展开式不是近似的，而是精确地。海赛矩阵是一个常矩阵，其中各元素均为常数。因此，无论从任何点出发，只需一步就可以找到极小点。因为若某一迭代法能使二次型函数在有限次迭代内达到极小点，则称此迭代方法是二次收敛的，因此牛顿方法是二次收敛的。 从牛顿法迭代公式的推演中可以看到，迭代点的位置是按照极值条件确定的，其中并未含有沿下降方向搜寻的概念。因此对于非二次函数，如果采用上述牛顿法公式，有时会使函数值上升，即出现1>k k f f +（x ）（x ） 现象。为此对上述牛顿方法进行改进，引入数学规划法的概念。 如果把1 2()()k k k d f x f x -??=-????看作是一个搜索方向，则采取如下的迭代公式 121()()k k k k k k k k x x a d x a f x f x -+??=-=-???? （k=0,1,2,…） 式中 k a —沿牛顿方向进行以为搜索的最佳步长k a 可通过如下极小化过程求得1()()()min k k k k k k k a f x f x a d f x a d +=+=+。由于此种方法每次迭代都在牛顿方向上进 行一维搜索，这就避免了迭代后函数值上升的现象，从而保持了牛顿法二次收敛的特性，而对初始点的选取并没有苛刻的要求。其计算步骤如下：
中北大学数值分析小论文
中北大学 《数值分析》 常微分方程初值问题的数值解法 专业： 班级： 学号： 姓名: 日期: 2012.12.26
常微分方程初值问题的数值解法 摘 要 微分方程的数值解法在科学技术及生产实践等多方面应用广泛. 文章分析了构造常微分方程初值问题数值解法的三种常用基本方法，差商代替导数法，数值积分法及待定系数法，推导出了Euler 系列公式及三阶龙格-库塔公式，指出了各公式的优劣性及适用条件，并对Euler 公式的收敛性、稳定性进行了分析。 Abstract The numerical solution of differential equations is widely used in science, technology, production practices and many other fields. This paper analyzed three kinds of basic methods for constructing numerical solutions for initial value problem of ordinary differential equations ：difference quotient instead of derivative method, numerical integral method and undetermined coefficients method. At the same time, the paper deduces the Euler series formula and the classical third order Runge-Kutta formula. In addition, the paper pointed out the advantages and disadvantages of each formula and application condition, it also analyzed the convergence and stability of the Euler formula. 1.引言 科学技术及实际生产实践中的许多问题都可归结为微分方程的求解问题，使用较多的是常微分方程初值问题的求解。对于一阶常微分方程的初值问题 000dy /dx f (x,y),y(x )y ,x x b ==<<，其中f 为已知函数，0y 是初始值。如 果函数f 关于变量y 满足Lipschitz 条件，则初值问题有唯一解。只有当f 是一些特殊类型的函数时，才能求出问题的解析解，但一般情况下都满足不了生产实践与科学技术发展的需要，因此通常求其数值解法。 2.主要算法 数值解法是一种离散化的方法，可以求出函数的精确解在自变量一系列离散点处的近似值。基本思想是离散化，首先要将连续区间离散化，对连续区域[]0x ,b 进行剖分01n 1n x x x x b -<<Λ<<=，n n 1n h x x +=-为步长；其次将其函离散
数值分析课程报告
插值法和多项式拟合的研究 摘要 在科研和生产实践中，常常需要通过一组测量数据来寻找变量x与y的函数关系近似表达式。解决这类问题的方法有两种：一种是插值法,另一种是拟合法。插值法的原理是用一个简单函数逼近被计算函数，然后用该简单函数的函数值近似替代被计算函数的函数值。拟合法能够是从给定的一组实验数据出发，寻找函数的一个近似表达式,该近似表达式能反映数据的基本趋势而又不一定过全部的点，即曲线拟合。本文主要介绍拉格朗日插值法、埃尔米特插值法、三次样条插值法以及基于最小二乘法的多项式拟合。 关键词：拉格朗日插值，埃尔米特插值，样条插值，多项式拟合
1方法的意义 在许多实际问题及科学研究中，因素之间往往存在着函数关系，然而，这种关系经常很难有明显的解析表达，通常只是由观察与测试得到一些离散数值。有时，即使给出了解析表达式，却由于表达式过于复杂，不仅使用不便，而且不易于进行计算与理论分析。解决这类问题的方法有两种：一种是插值法,另一种是拟合法。插值法的原理是用一个简单函数逼近被计算函数，然后用该简单函数的函数值近似替代被计算函数的函数值。它要求给出函数的一个函数表，然后选定一种简单的函数形式，比如多项式、分段线性函数及三角多项式等，通过已知的函数表来确定一个简单的函数()x ?作为()f x 的近似，概括地说，就是用简单函数为离散数组建立连续模型。插值法在实际应用中非常广泛，但是它也有明显的缺陷，一是测量数据常常带有测试误差，而插值多项式又通过所有给出的点，这样就是插值多项式保留了这些误差；二是如果实际得到的数据过多，则必然得到次数较高的插值多项式，这样近似的效果并不理想。拟合法能够很好的解决这些问题，它从给定的一组实验数据出发，寻找函数的一个近似表达式y=()x ?,该近似表达式能反映数据的基本趋势而又不一定过全部的点，即曲线拟合的问题，函数的近似表达式y=()x ?称为拟合曲线。常用最小而二乘法来确定拟合曲线。 2插值法的介绍 2.1 插值法定义 设 f (x )为[a ，b ]上的函数，在互异点n x x x ,...,,10处的函数值分别为 )(),...,(),(10n x f x f x f ,构造一个简单函数 ?(x ) 作为函数 f (x ) 的近似表达式y = f (x ) ≈ ?(x ),使 )()(i i x f x =? , i ＝0, 1, 2, …,n (1.0) 则称?(x ) 为关于节点n x x x ,...,,10的插值函数；称n x x x ,...,,10 为插值节点；称 ))((i i x f x , i =1,2,… , n 为插值点；f (x ) 称为被插值函数。式（1.0）称为插值条 件。这类问题称为插值问题。插值的任务就是由已知的观测点,为物理量(未知量)建立一个简单的、连续的解析模型，以便能根据该模型推测该物理量在非观测点
数值分析小论文
基于MATLAB曲线拟合对离散数据的处理和研究 摘要：曲线拟合是数值分析中的一种普遍且重要的方法，求解拟合曲线的方法也有很多，这里主要介绍利用MATLAB曲线拟合工具箱对离散数据点做你和处理，并与利用最小二乘法求相应的拟合曲线的方法做对比，突出MATLAB曲线拟合工具箱的优点，并阐述了其适用的范围，最后通过利用MATLAB曲线拟合工具箱对实例中离散数据点的拟合来具体说明它的使用方法和优点。 关键字：数值分析；MATLAB；曲线拟合；最小二乘法 一问题探究 在很多的实际情况中，两个变量之间的关系往往很难用具体的表达式把它表示出来，通常只能通过实际测量得到一些互不相同的离散数据点，需需要利用这些已知的数据点估计出两个变量的关系或工件的具体轮廓，并要得到任意未知数据点的具体数据，这个过程就需要用到拟合或差值方法来实现，这里主要讨论拟合的方法。 曲线拟合可以通过MATLAB编程来完成，通常为了达到更好的讷河效果需要做多次重复修改，对于非线性曲线拟合还需要编写复杂的M-文件，运用MATLAB曲线拟合工具箱来实现离散数据点的曲线拟合是一种直观并且简洁的方法。 二曲线拟合的最小二乘法理论 假设给定了一些数据点（Xi,Yi），人们总希望找到这样的近似的函数，它既能反映所给数据的一般趋势，又不会出现较大的偏差，并且要使构造的函数与被逼近函数在一个给定区间上的偏差满足某种要求。这种思想就是所谓的“曲线拟合”的思想。 曲线拟合和差值不同，若要求通过所有给定的数据点是差值问题，若不要求曲线通过所有给定的数据点，而只要求反映对象整体的变化趋势，拟合问题，曲线拟合问题最常用的解决方法是线性最小二乘法[1],步骤如下： 第一步：先选定一组函数r1(x),r2(x),…,rm(x),m