程序设计解决数学问题的研究.doc

用程序设计解决数学问题课题组成员：张昊誉、聂俊奇、姚望一、研究背景作为曾经的程序设计竞赛的一分子，程序设计社团的一分子，我们无可非议地选择了本课题的研究。二、研究目的及意义增强自己对理论知识的实践运用能力，要知道，在如今的课堂上，我们学习的都是理论知识，而不是动手实践能力，只知道，但不会用，那么学了又有什么作用呢？因此，为了锻炼这种能力，我们选择了本课题的研究。在这里，我们的专业知识能够对研究起到极佳的辅助作用，本研究活动也能加深我们对程序设计这一知识的更深入的掌握和运用。三、研究方法1作为一门比较冷门的知识，程序设计比较不为大众所了解，因此，我们取消了常规的的问卷调查方法，而是集中在少数掌握了本门知识的人群中，进行单个的调查。2同时，最重要的则是自己动手编写程序，解决学习过程中的一些用常规方法无法有效解决的数学问题。3在本过程中，我们遇到了许多困难，这要求我们依靠向老师求助与组员交流解决。4网络也是一个非常重要的知识来源，在这里，我们可以将自己的疑问与收获与众多有相同志趣的网友分享，加深我们对程序设计的理解。5、 研究过程：时间地点内容学校图书馆图书馆中查找相关文献学校机房学些编程相关知识学校机房上机实践，比较各算法差距寒假网上分头找资料，汇总讨论寒假网上向大学学长提问、讨论学校小结研究情况，准备中期评估一用程序计算圆周率数学上最迷人之一的东西就是神秘的圆周率，于是，我们。（你懂的）首先，圆周率是什么！？有谁能回答？呵呵，恐怕真的没人能准确说出来。请问数学老师知道么。以下是我们在网上搜索的结果：圆周率，一般以来表示，是一个在数学及物理学普遍存在的数学常数。它定义为圆形之周长与直径之比。它也等于圆形之面积与半径平方之比。是精确计算圆周长、圆面积、球体积等几何形状的关键值。 在分析学上，可以严格地定义为满足sin(x) = 0的最小正实数x。那么古人又是如何计算圆周率的的呢？带着这一疑问，我们首先对周围的同学们进行了调查。（调查中。）A：同学，你好，能回答我们一个问题吗？B：。好吧。问吧A：你知道古人是如何计算圆周率的吗？B：。A：好吧，但还是非常感谢你的参与，谢谢。A：同学，你好，能回答我们一个问题吗？C：恩？A：你知道古人是如何计算圆周率的吗？C：。听说过，好象是什么用圆的内切和外切多边形什么的。A：恩，好的，非常感谢你的参与！A：同学，你好，能回答我们一个问题吗？D：可以，但别问我隐私！A：。- -!不是不是，你想多了。你知道古人是如何计算圆周率的吗？D：我就知道祖之冲。（众人：祖冲之！ D：寒。）A：谢谢。看来同学们都对此不是很了解，那么数学老师呢？（于是，我们采访了高二数学组某数学老师）数学老师不愧是数学老师，TA果断给出了回答：用圆内接和外切正多边形的周长确定圆周长的上下界。（以下省略100字）好了，让我们看看网上的权威资料是怎么说的：(pai)是第十六个希腊字母，本来它是和圆周率没有关系的，但大数学家欧拉在一七三六年开始，在书信和论文中都用来代表圆周率。既然他是大数学家，所以人们也有样学样地用来表示圆周率了。但除了表示圆周率外，也可以用来表示其他事物，在统计学中也能看到它的出现。=Pai（=Pi） 古希腊欧几里德几何原本（约公元前3世纪初）中提到圆周率是常数，中国古算书周髀算经（ 约公元前2世纪）中有“径一而周三”的记载，也认为圆周率是常数。历史上曾采用过圆周率的多种近似值，早期大都是通过实验而得到的结果，如古埃及纸草书（约公元前1700）中取pi=（4/3）43.1604 。第一个用科学方法寻求圆周率数值的人是阿基米德，他在圆的度量（公元前3世纪）中用圆内接和外切正多边形的周长确定圆周长的上下界，从正六边形开始，逐次加倍计算到正96边形，得到（3+（10/71）（3+（1/7） ，开创了圆周率计算的几何方法（亦称古典方法，或阿基米德方法），得出精确到小数点后两位的值。 中国数学家刘徽在注释九章算术（263年）时只用圆内接正多边形就求得的近似值，也得出精确到两位小数的值，他的方法被后人称为割圆术。他用割圆术一直算到圆内接正192边形，得出根号10 （约为3.16）。 南北朝时代著名数学家祖冲之进一步得出精确到小数点后7位的值（约5世纪下半叶），给出不足近似值3.1415926和过剩近似值3.1415927，还得到两个近似分数值，密率355/113和约率22/7。他的辉煌成就比欧洲至少早了1000年。其中的密率在西方直到1573才由德国人奥托得到，1625年发表于荷兰工程师安托尼斯的著作中，欧洲不知道是祖冲之先知道密率的，将密率错误的称之为安托尼斯率。 阿拉伯数学家卡西在15世纪初求得圆周率17位精确小数值，打破祖冲之保持近千年的纪录。 德国数学家柯伦于1596年将值算到20位小数值，后投入毕生精力，于1610年算到小数后35位数，该数值被用他的名字称为鲁道夫数。 无穷乘积式、无穷连分数、无穷级数等各种值表达式纷纷出现，值计算精度也迅速增加。1706年英国数学家梅钦计算值突破100位小数大关。1873 年另一位英国数学家尚可斯将值计算到小数点后707位，可惜他的结果从528位起是错的。到1948年英国的弗格森和美国的伦奇共同发表了的808位小数值，成为人工计算圆周率值的最高纪录。（/view/3287.htm）在查找中，我们发现，目前计算圆周率有很多公式，如马青公式 拉马努金公式 AGM（Arithmetic-Geometric Mean）算法 波尔文四次迭代式 bailey-borwein-plouffe算法 丘德诺夫斯基公式 莱布尼茨公式 ，而研究再三，我们决定采用算法是用泰勒公式计算反正切值。于是，说干就干，我们很快写出了程序如下：#include #include main(int argc, char * argv) long * pi, * t, m, n, r, s; int t03 = 48, 32, 20, 24, 8, 4, k03 = 1, 1, 0, 1, 1, 1; int n03 = 18, 57, 239, 8, 57, 239, d, i, j, k, p, q; d = (argc 1) ? (i = atoi(argv1) 2) ? 1 : 0; printf(%snn, Nature (R) Pi value compute Program (C) Tue 1999.11.30); printf(pi= %s%d * arctg(1/%d) %s %d * arctg(1/%d) %s %d * arctg(1/%d) %sn, k0q0 ? : -, t0q0, n0q0, k0q1 ? + : -, t0q1, n0q1, k0q2 ? + : -, t0q2, n0q2, q ? Stomer : Gauss); if (t = (long *)calloc(d += 5) + 1, sizeof(long) = NULL) return 1; if (pi = (long *)calloc(d + 1, sizeof(long) = NULL) return 2; for (i = d; i = 0; i-) pii = 0; for (p = 0; p = 0; i-) ti = 0 for (r = 0, i = j; i = 0; i-) r = (m = 10 * r + ti) % n; ti = m / n; k ? (pii += ti) : (pii -= ti); while (j 0 & tj = 0) j-; for (k = !k, s = 3, n *= n; j 0; k = !k, s += 2) for (r = 0, i = j; i = 0; i-) r = (m = 10 * r + ti) % n; ti = m / n; while (j 0 & tj = 0) j-; for (r = 0, i = j; i = 0; i-) r = (m = 10 * r + ti) % s; m /= s; k ? (pii += m) : (pii -= m); for (n = i = 0; i = d; pii+ = r) n = (m = pii + n) / 10; if (r = m % 10) = 5; i-) printf(%ld%s, pii, (m = d - i + 5) % 65) ? (m % 5) ? : ) : n); printf(%sDIGITS: %dn, (m % 65) ? n : , d - 5); return 0;而指导老师姜老师也对这个程序提出了自己的看法：在命令行不跟参数执行该程序则使用Gauss公式计算前1000位圆周率的值，如果带一个命令行参数，则该值为要计算的位数。如果还有第二个命令行参数，则使用Stomer公式计算，可作为验算。因为该程序只涉及到纯数学计算，可以在Linux、Unix、Windows等操作系统下编译并运行。当时写这个程序时，int是2个字节的，现在大多数的C编译器int都是4个字节，不过这不影响程序的正确性。相信大家都看得很迷茫吧。其实我们自己也是！于是，我们不管三七二十一，便运行了程序，但是却得到了程序出错的提示！Error!众人顿时倒吸了一口冷气，倒不是因为它的破坏有多大，而是因为要检查出此错误要花相当多的时间，根据专业的程序员的说法：他们88%的时间都在校验程序的正确性，只有12%的时间用于录入程序！但命运之神还是垂青于我们的，呵呵，我们很快发现了这个bug，那么到底是什么呢？请先让我们卖个关子，以下是修改后的程序：#include #include main(int argc, char * argv) long * pi, * t, m, n, r, s; int t03 = 48, 32, 20, 24, 8, 4, k03 = 1, 1, 0, 1, 1, 1; int n03 = 18, 57, 239, 8, 57, 239, d, i, j, k, p, q; d = (argc 1) ? (i = atoi(argv1) 2) ? 1 : 0; printf(%snn, Nature (R) Pi value compute Program (C) Tue 1999.11.30); printf(pi= %s%d * arctg(1/%d) %s %d * arctg(1/%d) %s %d * arctg(1/%d) %sn, k0q0 ? : -, t0q0, n0q0, k0q1 ? + : -, t0q1, n0q1, k0q2 ? + : -, t0q2, n0q2, q ? Stomer : Gauss); if (t = (long *)calloc(d += 5) + 1, sizeof(long) = NULL) return 1; if (pi = (long *)calloc(d + 1, sizeof(long) = NULL) return 2; for (i = d; i = 0; i-) pii = 0; for (p = 0; p = 0; i-) ti = 0; for (r = 0, i = j; i = 0; i-) r = (m = 10 * r + ti) % n; ti = m / n; k ? (pii += ti) : (pii -= ti); while (j 0 & tj = 0) j-; for (k = !k, s = 3, n *= n; j 0; k = !k, s += 2) for (r = 0, i = j; i = 0; i-) r = (m = 10 * r + ti) % n; ti = m / n; while (j 0 & tj = 0) j-; for (r = 0, i = j; i = 0; i-) r = (m = 10 * r + ti) % s; m /= s; k ? (pii += m) : (pii -= m); for (n = i = 0; i = d; pii+ = r) n = (m = pii + n) / 10; if (r = m % 10) = 5; i-) printf(%ld%s, pii, (m = d - i + 5) % 65) ? (m % 5) ? : ) : n); printf(%sDIGITS: %dn, (m % 65) ? n : , d - 5); return 0;相信你即使看了好几遍也不一定找得到到底哪里发生了改动，那就让我们来告诉你吧！在 for (k=k0qp, n=n0qp, ti=j=d=t0qp, i-; i = 0; i-) ti = 0;一句中，我们漏了个;！是的，就是这么简单！一个字符的缺失，竟导致了整个程序的崩溃！经过此次教训，我们在之后的程序录入中，变得更加小心谨慎，避免犯下类似的错误。耕耘之后，总会有收获的，当我们运行了程序之后，也如愿以偿地获得了我们想要的结果：Nature (R) Pi value compute Program (C) Tue 1999.11.30pi= 48 * arctg(1/18) + 32 * arctg(1/57) - 20 * arctg(1/239) Gausspi= 3.14159 26535 89793 23846 26433 83279 50288 41971 69399 37510 58209 7494459230 78164 06286 20899 86280 34825 34211 70679 82148 08651 32823 06647 0938446095 50582 23172 53594 08128 48111 74502 84102 70193 85211 05559 64462 2948954930 38196 44288 10975 66593 34461 28475 64823 37867 83165 27120 19091 4564856692 34603 48610 45432 66482 13393 60726 02491 41273 72458 70066 06315 5881748815 20920 96282 92540 91715 36436 78925 90360 01133 05305 48820 46652 1384146951 94151 16094 33057 27036 57595 91953 09218 61173 81932 61179 31051 1854807446 23799 62749 56735 18857 52724 89122 79381 83011 94912 98336 73362 4406566430 86021 39494 63952 24737 19070 21798 60943 70277 05392 17176 29317 6752384674 81846 76694 05132 00056 81271 45263 56082 77857 71342 75778 96091 7363717872 14684 40901 22495 34301 46549 58537 10507 92279 68925 89235 42019 9561121290 21960 86403 44181 59813 62977 47713 09960 51870 72113 49999 99837 2978049951 05973 17328 16096 31859 50244 59455 34690 83026 42522 30825 33446 8503526193 11881 71010 00313 78387 52886 58753 32083 81420 61717 76691 47303 5982534904 28755 46873 11595 62863 88235 37875 93751 95778 18577 80532 17122 6806613001 92787 66111 95909 21642 01989DIGITS: 1000（这些数字可不是我们一个一个照着结果打出来的哦！只需要对程序稍加修改，就可以让它直接将结果输出在txt文件里，这样就省力多了，哈哈）但是，直接照搬前人的公式的难度并不大，在姜老师的要求下，我们决定采用自己的方法进行又一次编程。这是一个直角坐标系，我们可以看到如果满足X方+Y方小于等于r方的那就可以判断点在圆内（这里为了方便编程，所以用了四分之一圆）图上s1和s2的比值正好是0.25我们可以用模拟试验的方法，在取很多个点中记下在圆内的点的个数，通过在圆内与总的点的个数比值就可以估算s1和s2的比值自然也就可以算出来了当然，这种算法效率很低，精度也低但原理十分简单，十分适合我们上手编程程序如下：#include #include #include /使用当前时钟做种子 void main( void ) int n=1;double x,y;double sum1=0,sum2=0;long i,j; srand( (unsigned)time( NULL ) ); /初始化随机数 printf(你想试几次？);scanf(%d,&n);for(j=1;j=n;j+)for( i=1; i=1000000;i+ ) /打印出n个随机数x=rand();y=rand();if(x*x+y*y=32767*32767)/距离公式再熟悉不过了，这边32767是随机数的最大范围sum1+;/计数范围内的点sum1=4*sum1/n;/printf(%lfn,sum1);sum2=sum2+sum1;sum1=0;printf(n%lfn,sum2/1000000);/最后取个平均值具体见我们的视频。就这样，我们顺利地完成了第一个目标，大家都十分兴奋，而姜老师也对我们的成功表示肯定，并鼓励我们向更高难度冲刺二二分法大家应该不会对二分法感到陌生吧！那恶心的计算量一定会让每个同学感到无可奈何！但是，我们程序设计的最大优势就是计算量巨大而且又快又对！因此，二分法无疑是程序设计的一个比较好的运用。而姜老师也相当支持我们的第二个目标。可能有些人对二分法不熟悉，或者已经忘了！（我们就碰到了这个小尴尬）因此，在这里，我们从网上查找了一些相关资料帮助大家了解一般地，对于函数f(x),如果存在实数c,当x=c时，若f(c)=0,那么把x=c叫做函数f(x)的零点。 解方程即要求f(x)的所有零点。 假定f(x)在区间（x，y）上连续 先找到a、b属于区间（x，y），使f(a)，f(b)异号，说明在区间(a,b)内一定有零点，然后求f(a+b)/2, 现在假设f(a)0,ab 如果f(a+b)/2=0，该点就是零点， 如果f(a+b)/2a，从开始继续使用 中点函数值判断。 如果f(a+b)/20，则在区间(a,(a+b)/2)内有零点，(a+b)/2=b，从开始继续使用 中点函数值判断。 这样就可以不断接近零点。 通过每次把f(x)的零点所在小区间收缩一半的方法，使区间的两个端点逐步迫近函数的零点，以求得零点的近似值，这种方法叫做二分法。 从以上可以看出，每次运算后，区间长度减少一半，是线形收敛。另外，二分法不能计算复根和重根。 （/view/75441.htm）首先还是常规采访：A：同学！你在干什么呢？E：没看到我在做数学呢！A：哦？用二分法求解方程的近似解？E：是的，好恶心啊A：我告诉你，用程序设计可以很快的做出这道题目哦E：程序设计是什么？可以吃吗？A：。E：开玩笑，呵呵。不过我可不会那什么程序，只能自己死算喽（前面部分同上）F：哦？那你快点去编一下把答案告诉我！A：。（前面部分同上）G：这是锻炼我计算能力的大好机会，让电脑做不是浪费了么！？（众人：好学生啊！ 老师：T-T感动啊）A：。好吧，那您慢做然后是老师：X老师：二分法不考的，你放心好了！A：。（赶紧溜）我们首先先写好编程思路（这真的很重要）给定精确度,用二分法求函数f(x)零点近似值的步骤如下: 1 确定区间a,b,验证f(a)f(b)0,给定精确度. 2 求区间(a,b)的中点c. 3 计算f(c). (1) 若f(c)=0,则c就是函数的零点; (2) 若f(a)f(c)0,则令b=c; (3) 若f(c)f(b)0,则令a=c. (4) 判断是否达到精确度:即若f(a)或者f(b),则得到零点近似值a(或b),否则重复2-4.恩，闲话不多说了，直接动手编程，以下是我们编写的一个比较简单的程序#include#define f(x) (x*x*x-2*x*x+3*x-4)void main() float a=-10,b=10,c,eps=1e-5; while (b-a)eps) c=(a+b)/2; if(f(c)=0) break;else if(f(a)*f(c)0) b=c; else a=c; printf(root=%fn,c);这段代码是求解方程f(x)=0在区间-10,10上的根的数值解。方法的思想就是：一直选取区间中间的数值，如果发现中间的函数值与一侧函数值，异号，那么说明解在这个更小的区间中，采用eps=1e-5作为区间的极限大小，通过迭代的方法求解这个方程的数值解。是的，就这么简单，一小段程序，在不到1S的时间内，就能把我们笔算要花很久的问题瞬间求解出来。而其实，其他复杂的二分法程序，与我们上面写的那段小程序原理都是相同的。由于我们的目的是对学习中可能遇到的二分法用编程进行解决，因此，就不深入研究高次方程二分法求解的过程了。但是，二分法其实还有其他用处！比如，在程序设计中，用二分法进行排序，可以比常规的排序方法快很多，因此也被称为快排其基本思路也用到了二分法的原理以下是我们小组成员写的快排的伪代码（写起来比较简单，而且大家也比较好看懂）QUICKSORT(A，p，r) 1ifpr 2thenq PARTITION(A，p，r) 3QUICKSORT(A，p，q-1) 4QUICKSORT(A，q+1，r) 为排序一个完整的数组A，最初的调用是QUICKSORT(A，1，lengthA)。 快速排序算法的关键是PARTITION过程，它对子数组Ap.r进行就地重排： PARTITION(A，p，r) 1x Ar 2i p-1 3forj ptor-1 4do ifAjx 5theni i+1 6exchange AiAj 7exchange Ai+1Ar 8returni+1当然啦，这些都是在专业编程中所使用的，在我们的研究性学习中，作用并不是很大，因此不再多做叙述了三其他一些数学问题（杂题）递归的运用：递归做为一种算法在程序设计语言中广泛应用.是指函数/过程/子程序在运行过程中直接或间接调用自身而产生的重入现象.递归是计算机科学的一个重要概念,递归的方法是程序设计中有效的方法,采用递归编写程序能使程序变得简洁和清晰.。恩，这个相信给大家说了概念也不能很好的解释，给大家举个实例把已知f(n)=n-1 求f(f(.f(f(10000).)的值要求整个式子的值，我们要先求出f(10000)的值，再将结果代回到f(n)，既f(9999)，再求f(9998)。直到求出最后一个f(n)的值，这样就得出了最后的结果。当然，这个问题其实相当简单，不需要我们用程序设计来解决，那么，来看以下这个问题 总共有100级楼梯，从第0级开始走，每次可以走1步或者2步，请问走完这些楼梯一共有多少种方法？采访同学中。H：不会。I：你啊无聊的啊！J：（做了半天）是不是这个答案？！ （A：其实我们自己也还没编程算呢。） 那么用程序设计该如何解决这个问题呢？根据上面的思路和题意，我们只要求出到第99级楼梯的走法数和第98级楼梯的走法数之和，那么到99级楼梯的走法有多少呢？我们又可以划分为到第98级楼梯和第97级楼梯的走法数之和。依次类推，知道最后划分为N个走到第1级楼梯有多少种走法，很明显，是1种。最后把结果相加即可。接下来是另一道题目海滩上有一堆桃子，五只猴子来分，第一只猴子把这堆桃子凭据分成五分，多了一个，这只猴子就把多的一个扔到了海里！，拿走了一份，第二只猴子把剩下的桃子又平均分成五分，又多了一个，它同样把多的一个扔到了海里，拿走了一份，第三只，第四只，第五只都是这样做的，问海滩上原来有多少桃子直接给程序#include int test(int total, int count) if(count = 6) /如果分五次之后仍满足要求，则该数为所求 return 1; if(total % 5 != 1) /如果不满足分五分后剩一个则该数不满足要求 return 0; else if( test(total-1) * 4 / 5, count+1) /将该数的5分之1继续尝试分 printf(第%d次分时还剩%d个桃子n, count, total); return 1; else return 0; int main() int i = 4; while(1) if(test(+i, 1) break; /看看此数是否能够符合题意 printf(the num is %dn, i); return 1;四不知不觉我们的研究也进入尾声了，除了自己编程解决数学问题外，我们对一些数学软件也十分感兴趣在种类繁多的数学软件中，我们了解到：数学软件基本分为三类： 1 数值计算的软件，如matlab（商业软件），scilab(开源自由软件）等等； 2 统计软件，如SAS（商业软件）、minitab（商业软件）、SPSS（商业软件），R（开源自由软件）等； 3 符号运算软件，这种是最绝妙的，不像前两种那样只能计算出数值，而是可以把符号表达成的公式、方程进行推导和化简，可以求出微分积分的表达式，代表有maple（商业软件）、mathematica（商业软件），maxima（开源自由软件），mathcad（商业软件）等等。其中，我们对mathematica十分感兴趣，而且在上学期的一次展示中，我们也曾经用这个软件进行了一些展示Mathematica是一款科学计算软件，很好地结合了数值和符号计算引擎、图形系统、编程语言、文本系统、和与其他应用程序的高级连接。很多功能在相应领