等无理数近似计算与分析

浙江科技学院本科毕业设计 论文 I 等无理数近似计算与分析等无理数近似计算与分析 众所周知 圆周率是平面上圆的周长与直径之比 其值为 3 1415926 在生产实践中 经常遇到求圆周长 圆面积 球体积等 这类问题都接触到圆周率 的计算和理论研究反映 了一个民族的数学水平 对古代人来说 的计算是一件繁重的工作 古代人们把 3 作为它的 近似值 古希腊阿基米德曾得到 我国宋代的祖冲之得到的近似值 约率 101 33 717 22 7 和 密率 后者化为小数后等于 3 141592 与的准确值的误差在以下 迄今为止 355 113 6 10 仍被认为是理想的近似值 355 113 后来 数学家们不满足于此 他们认为不可能是一个分式表示的有理数 并且在上个世 纪 丹麦数学家康托证明了属于无理数中的 超越数 它的真值是一个无限不循环小数 像 这样神妙莫测的无理数比比皆是 在漫长几千年的历史中许多数学家都对这类数的近似计算 作不断的研究 出现了很多计算无理数近似值的方法 随着角的弧度制的采用 微积分的发现 特别是计算机的发明 计算等无理数的位数大大地增加了 的近似计算已提高到了千亿位 以上 在这样一个计算机科学迅速发展的时代 将数学思想 方法 理论知识与计算机技术更好 地链接并使数值计算方法得到更好的见证与应用无疑是一个值得讨论的课题 数值分析的主要 内容是研究适合用计算机求解各种数学问题的数值计算方法 包括构造算法及数值计算的误差 分析问题 随着计算机科学的迅速发展 数值计算方法已广泛地应用于自然科学 工程技术的 各个领域之中 它是求解各种数学问题进行科学计算必不可少的主要方法 运用所学知识并借 鉴前人的方法将近似计算的方法应用到数学实验中 通过 Mathematica 数学软件来实现其结果 将数学思想渗透到高科技领域 必将有新的收获 浙江科技学院本科毕业设计 论文 2 第一章第一章 现状与分析现状与分析 随着计算机技术的迅速发展 数值方法在各个技术领域中的应用越来越广泛 并已成 为数学与计算机之间的桥梁 求解一个科学技术的实际问题都是先把实际问题抽象近似地归纳为数学问题 建立一 个数学模型 因直接求解数学问题很难得到原问题能用数学解析式表达的解析解 故需要 通过离散化将数学问题转换成一个数值问题 这才是真正实用而有效的方法 所谓数值问 题是指可以在原问题自变量的原始数据和用输出数据表示的结果之间建立的一个确定的函 数关系 对该数值问题进行数值运算即可得到数值解 即计算解或近似解 在计算机成为数值计算主要工具的今天 求解数值问题必须要构造适合于计算机的数 值方法 这种数值方法是指求解数值问题的 按确定顺序可以在计算机上依次执行的计算 公式 即通常所说的数值计算方法 求解同一个数学问题可以构造多种算法 选取不同的 数值方法 得到的结果会有很大差别 其好坏与原数学问题的性态有关 与数值方法本身 的稳定性 计算解的精度和可靠性 以及计算量的大小有关 研究如何较好地用计算机求 解数学问题的数值方法和理论是数值分析的主要内容 它是数学的一个重要分支 其内容 不像纯数学那样只研究理论 而是着重研究求解的数值方法和相关的理论 这些理论包括 方法的收敛性 稳定性和误差分析 具有高度抽象性 严密科学性 广泛应用性 与计算 机科学的密切性以及实用性 数值计算方法的重要性已越来越显著 浙江科技学院本科毕业设计 论文 3 第二章第二章 计算历程与研究意义计算历程与研究意义 2 1 的历史 1 圆周率是一个极其驰名的数 从有文字记载的历史开始 这个数就引进了外行人和学 者们的兴趣 作为一个非常重要的常数 圆周率最早是出于解决有关圆的计算问题 仅凭 这一点 求出它的尽量准确的近似值 就是一个极其迫切的问题了 事实也是如此 几千 年来作为数学家们的奋斗目标 古今中外一代代的数学家为此献出了自己的智慧和劳动 回顾历史 人类对的认识过程 反映了数学和计算技术发展情形的一个侧面 的研究 在一定程度上反映这个地区或时代的数学水平 德国数学史家康托说 历史上一个国家所 算得的圆周率的准确程度 可以作为衡量这个国家当时数学发展水平的指标 下面就计算 的历史过程讨论一下值的计算方法 几何学几何学 若圆的半径为 r 圆周为 2Cr 若圆的半径为 r 其面积为 2 Ar 若球的半径为 r 其体积为 3 43 Vr 若球的半径为 r 其表面积为 r 2 4Ar 分析数学 10 i e 数论数论 任意两个自然数 互质的概率是 6 概率论 概率论 取一枚长为 l 的针 再取一张白纸在上面画上一些距离为 2 的平行线 把针从一定高 度释放 让其自由落体到纸面上 针与平行线相交的概率是圆周率的倒数 蒲丰投针 物理学 物理学 海森堡测不准原理 4 h x p A A 爱因斯坦相对论场方程 4 8 2 ik ikikik g RG RgT c 实验时期计算实验时期计算 通过实验对值进行估算 这是计算的的第一阶段 这种对值的估算基本上都是 以观察或实验为根据 是基于对一个圆的周长和直径的实际测量而得出的 在古代世界 实际上长期使用这个数值 最早见于文字记载的有基督教 圣经 中的章节 其上3 取圆周率为 3 这一段描述的事大约发生在公元前 950 年前后 其他如巴比伦 印度 中 浙江科技学院本科毕业设计 论文 4 国等也长期使用 3 这个粗略而简单实用的数值 在我国刘徽之前 圆径一而周三 曾广泛 流传 我国第一部 周髀算经 中 就记载有圆 周三径一 这一结论 在我国 木工师 傅有两句从古流传下来的口诀 叫做 周三径一 方五斜七 意思是说 直径为 1 的圆 周长大约是 3 边长为 5 的正方形 对角线之长约为 7 这正反映了早期人们对圆周率 和这两个无理数的粗略估计 东汉时期官方还明文规定圆周率取 3 为计算面积的标 2 准 后人称之为 古率 早期的人们还使用了其它的粗糙方法 如古埃及 古希腊人曾用谷粒摆在圆形上 以 数粒数与方形对比的方法取得数值 或用匀重木板锯成圆形和方形以秤量对比取值 由 此 得到圆周率的稍好些的值 如古埃及人应用了约四千年的 4 8 9 2 3 1605 在印度 公元前六世纪 曾取 在我国东 西汉之交 新朝王莽令刘歆制造量的103 162 容器 律嘉量斛 刘歆在制造标准容器的过程中就需要用到圆周率的值 为此 他大约 也是通过做实验 得到一些关于圆周率的并不划一的近似值 现在根据铭文推算 其计算 值分别取为 3 1547 3 1992 3 1498 3 2031 比径一周三的古率已有所进步 人类的这种探 索的结果 当主要估计圆田面积时 对生产没有太大影响 但以此来制造器皿或其它计算 就不合适了 几何法计算几何法计算 凭直观推测或实物度量 来计算值的实验方法所得到的结果是相当粗略的 真正使 圆周率计算建立在科学的基础上 首先应归功于阿基米德 他是科学地研究这一常数的第 一个人 是他首先提出了一种能够借助数学过程而不是通过测量的 能够把的值精确到 任意精度的方法 由此 开创了圆周率计算的第二阶段 图图 2 1 的值域 2 24 圆周长大于内接正四边形而小于外切正四边形 因此 当然 这是一个差2 24 劲透顶的例子 据说阿基米德用到了正 96 边形才算出他的值域 阿基米德求圆周率更精确的近似值的方法 体现在他的一篇论文 圆的测定 之中 在这一书中 阿基米德第一次创用上 下界来确定的近似值 他用几何方法证明了 圆 周长与圆直径之比小于 3 1 7 而大于 3 10 71 他还提供了误差的估计 重要的是 浙江科技学院本科毕业设计 论文 5 这种方法从理论而言 能够求得圆周率更准确的值 到公元 150 年左右 希腊天文学家托 勒密得出 取得了自阿基米德以来的巨大进步 3 1416 图图 2 2 割圆术 割圆术 多次利用勾股定理 来计算正边形的边长 N 在我国 首先是由数学家刘徽得出较精确的圆周率 公元 263 年前后 刘徽提出著名 的割圆术 得出 通常称为 徽率 他指出这是不足近似值 虽然他提出割圆术3 14 的时间比阿基米德晚一些 但其方法确有着较阿基米德方法更美妙之处 割圆术仅用内接 正多边形就确定出了圆周率的上 下界 比阿基米德用内接同时又用外切正多边形简捷得 多 另外 有人认为在割圆术中刘徽提供了一种绝妙的精加工办法 以致于他将割到 192 边形的几个粗糙的近似值通过简单的加权平均 竟然获得具有 4 位有效数字的圆周率 而这一结果 正如刘徽本人指出的 如果通过割圆计算得出这个3927 12503 1416 结果 需要割到 3072 边形 这种精加工方法的效果是奇妙的 这一神奇的精加工技术是割 圆术中最为精彩的部分 令人遗憾的是 由于人们对它缺乏理解而被长期埋没了 可能大家更加熟悉的是祖冲之所做出的贡献吧 对此 隋书 律历志 有如下记载 宋末 南徐州从事祖冲之更开密法 以圆径一亿为丈 圆周盈数三丈一尺四寸一分五厘 九毫二秒七忽 朒数三丈一尺四寸一分五厘九毫二秒六忽 正数在盈朒二限之间 密率 圆径一百一十三 圆周三百五十五 约率 圆径七 周二十二 这一记录指出 祖冲之关于圆周率的两大贡献 其一是求得圆周率 3 14159263 1415927 其二是 得到的两个近似分数即 约率为 密率为 7 22355 113 他算出的的 8 位可靠数字 不但在当时是最精密的圆周率 而且保持世界记录九百 多年 以致于有数学史家提议将这一结果命名为 祖率 这一结果是如何获得的呢 追根溯源 正是基于对刘徽割圆术的继承与发展 祖冲之 才能得到这一非凡的成果 因而当我们称颂祖冲之的功绩时 不要忘记他的成就的取得是 因为他站在数学伟人刘徽的肩膀上的缘故 后人曾推算若要单纯地通过计算圆内接多边形 浙江科技学院本科毕业设计 论文 6 边长的话 得到这一结果 需要算到圆内接正 12288 边形 才能得到这样精确度的值 祖 冲之是否还使用了其它的巧妙办法来简化计算却已不得而知 因为记载其研究成果的著作 缀术 早已失传 这在中国数学发展史上是一件极令人痛惜的事 2 2 现代科学家计算值的动力与意义 现在打破记录 不管推进到多少位 也不会令人感到特别的惊奇了 实际上 把的 数值算得过分精确 应用意义并不大 现代科技领域使用的值 有十几位已经足够 如 果用鲁道夫的 35 位小数的值计算一个能把太阳系包围起来的圆的周长 误差还不到质 子直径的百万分之一 我们还可以引美国天文学家西蒙 纽克姆的话来说明这种计算的实用 价值 十位小数就足以使地球周界准确到一英寸以内 三十位小数便能使整个可见宇宙的 四周准确到连最强大的显微镜都不能分辨的一个量 那么为什么数学家们还象登山运动员那样 奋力向上攀登 一直求下去而不是停止对 的探索呢 为什么其小数值有如此的魅力呢 分析其主要原因有分析其主要原因有 它现在可以被人们用来测试或检验超级计算机的各项性能 特别是运算速度与计算 过程的稳定性 这对计算机本身的改进至关重要 就在几年前 当 Intel 公司推出奔腾 Pentium 时 发现它有一点小问题 这问题正是通过运行的计算而找到的 这正是超 高精度的计算直到今天仍然有重要意义的原因之一 计算的方法和思路可以引发新的概念和思想 虽然计算机的计算速度超出任何人的 想象 但毕竟还需要由数学家去编制程序 指导计算机正确运算 实际上 确切地说 当 我们把的计算历史划分出一个电子计算机时期时 这并非意味着计算方法上的改进 而 只是计算工具有了一个大飞跃而已 因而如何改进计算技术 研究出更好的计算公式 使 公式收敛得更快 能极快地达到较大的精确度仍是数学家们面对的一个重要课题 在这方 面 本世纪印度天才数学家拉马努扬得出了一些很好的结果 他发现了许多能够迅速而精 确地计算近似值的公式 他的见解开通了更有效地计算近似值的思路 现在计算机计 算值的公式就是由他得到的 这些极富传奇色彩的数学家与的故事讲述的是人类的胜 利 而不是机器的胜利 还有一个关于的计算的问题是 我们能否无限地继续算下去 答案是 不行 根据朱达偌夫斯基的估计 我们最多算 1077 位 虽然 现在我们离这一极限还相差很远很 远 但这毕竟是一个界限 为了不受这一界限的约束 就需要从计算理论上有新的突破 前面我们所提到的计算 不管用什么公式都必须从头算起 一旦前面的某一位出错 后面 的数值完全没有意义 那么 能否计算时不从头开始 而是从半截开始呢 这一根本性的 想法就是寻找并行算法公式 1996 年 圆周率的并行算法公式终于找到 但这是一个 16 浙江科技学院本科毕业设计 论文 7 进位的公式 这样很容易得出的 1000 亿位的数值 只不过是 16 进位的 是否有 10 进位的 并行计算公式 仍是未来数学的一大难题 作为一个无穷数列 数学家感兴趣的把展开到上亿位 能够提供充足的数据来 验证人们所提出的某些理论问题 可以发现许多迷人的性质 如 在的十进展开中 10 个数字 哪些比较稀 哪些比较密 的数字展开中某些数字出现的频率会比另一些高吗 或许它们并非完全随意 这样的想法并非是无聊之举 只有那些思想敏锐的人才会问这种 貌似简单 许多人司空见惯但却不屑发问的问题 数学家弗格森最早有过这种猜想 在的数值式中各数码出现的概率相同 正是 他的这个猜想为发现和纠正向克斯计算值的错误立下了汗马功劳 然而 猜想并不等于 现实 弗格森想验证它 却无能为力 后人也想验证它 也是苦于已知的值的位数太少 甚至当位数太少时 人们有理由对猜想的正确性做出怀疑 如 数字 0 的出现机会在开始 时就非常少 前 50 位中只有 1 个 0 第一次出现在 32 位上 可是 这种现象随着数据的 增多 很快就改变了 100 位以内有 8 个 0 200 位以内有 19 个 0 1000 万位以内有 999 440 个 0 60 亿位以内有 599 963 005 个 0 几乎占 1 10 其他数字又如何呢 结果显示 每一个都差不多是 1 10 有的多一点 有的少一点 虽然有些偏差 但都在 1 10000 之内 人们还想知道 的数字展开真的没有一定的模式吗 我们希望能够在十进制展 开式中通过研究数字的统计分布 寻找任何可能的模型 如果存在这种模型的话 迄今 为止尚未发现有这种模型 同时我们还想了解 的展开式中含有无穷的样式变化吗 或 者说 是否任何形式的数字排列都会出现呢 著名数学家希尔伯特在没有发表的笔记本中 曾提出下面的问题 的十进展开中是否有 10 个 9 连在一起 以现在算到的 60 亿位数字 来看 已经出现 连续 6 个 9 连在一起 希尔伯特的问题答案似乎应该是肯定的 看来任 何数字的排列都应该出现 只是什么时候出现而已 但这还需要更多的数位的计算才能 提供切实的证据 在这方面 还有如下的统计结果 在 60 亿数字中已出现连在一起的 8 个 8 9 个 7 10 个 6 小数点后第 710150 位与 3204765 位开始 均连续出现了七个 3 小数点 52638 位起连续出现了 14142135 这八个数字 这恰是的前八位 小数点后第 2747956 位起 出现 了有趣的数列 876543210 遗憾的是前面缺个 9 还有更有趣的数列 123456789 也出现了 如果继续算下去 看来各种类型的数列组合可能都会出现 通过几何 微积分 概率 等广泛的范围和渠道发现 这充分显示了数学方法的奇异美 竟然与这么些表面看来 风马牛不相及的试验 沟通在一起 这的确使人惊讶不已 难怪有这么多人探索他的奥妙 浙江科技学院本科毕业设计 论文 8 第三章第三章 近似值的三种数值计算方法近似值的三种数值计算方法 3 1 蒙特卡洛 Monte Carlo 方法 3 1 1 蒙特卡罗 Monte Carlo 方法简介 蒙特卡罗 Monte Carlo 方法 或称计算机随机模拟方法 是一种基于 随机数 的计算 方法 它实质上是利用服从某种分布的随机变量来模拟现实系统中可能出现的随机现象 在项目管理中 可以用来模拟计算不确定性很强的项目收益 进度和成本 以及评估不确 定因素对项目结果的影响 这一方法源于美国在第一次世界大战进研制原子弹的 曼哈顿计 划 该计划的主持人之一 数学家冯 诺伊曼用驰名世界的赌城 摩纳哥的 Monte Carlo 来命名这种方法 为它蒙上了一层神秘色彩 Monte Carlo 方法的基本思想很早以前就被人们所发现和利用 早在 17 世纪 人们就 知道用事件发生的 频率 来决定事件的 概率 19 世纪人们用投针试验的方法来决定圆周 率 本世纪 40 年代电子计算机的出现 特别是近年来高速电子计算机的出现 使得用 数学方法在计算机上大量 快速地模拟这样的试验成为可能 考虑平面上的一个边长为 1 的正方形及其内部的一个形状不规则的 图形 如何求出 这个 图形 的面积呢 Monte Carlo 方法是这样一种 随机化 的方法 向该正方形 随机地 投掷 N 个点落于 图形 内 则该 图形 的面积近似为 M N 可用民意测验来作一个不严格的比喻 民意测验的人不是征询每一个登记选民的意见 而是通过对选民进行小规模的抽样调查来确定可能的优胜者 其基本思想是一样的 科技计算中的问题比这要复杂得多 比如金融衍生产品 期权 期货 掉期等 的定 价及交易风险估算 问题的维数 即变量的个数 可能高达数百甚至数千 对这类问题 难度随维数的增加呈指数增长 这就是所谓的 维数的灾难 Course Dimensionality 传统 的数值方法难以对付 即使使用速度最快的计算机 Monte Carlo 方法能很好地用来对付 维数的灾难 因为该方法的计算复杂性不再依赖于维数 以前那些本来是无法计算的问题 现在也能够计算量 为提高方法的效率 科学家们提出了许多所谓的 方差缩减 技巧 3 1 2 蒙特卡罗 Monte Carlo 方法的基本原理 12 蒙特卡罗 Monte Carlo 方法 又称随机抽样或统计试验方法 属于计算数学的一个 分支 它是在本世纪四十年代中期为了适应当时原子能事业的发展而发展起来的 传统的 经验方法由于不能逼近真实的物理过程 很难得到满意的结果 而蒙特卡罗方法由于能够 浙江科技学院本科毕业设计 论文 9 真实地模拟实际物理过程 故解决问题与实际非常符合 可以得到很圆满的结果 这也是 我们采用该方法的原因 当所要求解的问题是某种事件出现的概率 或者是某个随机变量的期望值时 它们可 以通过某种 试验 的方法 得到这种事件出现的频率 或者这个随机变数的平均值 并用 它们作为问题的解 这就是蒙特卡罗方法的基本思想 蒙特卡罗方法通过抓住事物运动的 几何数量和几何特征 利用数学方法来加以模拟 即进行一种数字模拟实验 它是以一个 概率模型为基础 按照这个模型所描绘的过程 通过模拟实验的结果 作为问题的近似解 由概率定义知 某事件的概率可以用大量试验中该事件发生的频率来估算 当样本容 量足够大时 可以认为该事件的发生频率即为其概率 因此 可以先对影响其可靠度的随 机变量进行大量的随机抽样 然后把这些抽样值一组一组地代入功能函数式 确定结构是 否失效 最后从中求得结构的失效概率 蒙特卡罗法正是基于此思路进行分析的 设有统计独立的随机变量 其对应的概率密度函数分别为 1 2 3 i Xik 1 x f 功能函数式为 2 x f k x f 12 k Zg x xx 首先根据各随机变量的相应分布 产生 N 组随机数 值 计算功能函 1 x 2 x k x 数值 若其中有 L 组随机数对应的功能函数值 12 1 2 ik Zg x xxiN 0 i Z 则当时 根据伯努利大数定理及正态随机变量的特性有 结构失效概率 可靠指标 N 从蒙特卡罗方法的思路可看出 该方法回避了结构可靠度分析中的数学困难 不管状 态函数是否非线性 随机变量是否非正态 只要模拟的次数足够多 就可得到一个比较精 确的失效概率和可靠度指标 特别在岩土体分析中 变异系数往往较大 与 JC 法计算的 可靠指标相比 结果更为精确 并且由于思路简单易于编制程序 3 1 3 蒙特卡罗 Monte Carlo 解题主要步骤 蒙特卡罗解题归结为三个主要步骤 构造或描述概率过程 实现从已知概率分布抽样 建立各种估计量 3 1 3 1 构造或描述概率过程 对于本身就具有随机性质的问题 如粒子输运问题 主要是正确描述和模拟这个概率 过程 对于本来不是随机性质的确定性问题 比如计算定积分 就必须事先构造一个人为 的概率过程 它的某些参量正好是所要求问题的解 即要将不具有随机性质的问题转化为 随机性质的问题 浙江科技学院本科毕业设计 论文 10 3 1 3 2 实现从已知概率分布抽样 构造了概率模型以后 由于各种概率模型都可以看作是由各种各样的概率分布构成的 因此产生已知概率分布的随机变量 或随机向量 就成为实现蒙特卡罗方法模拟实验的基 本手段 这也是蒙特卡罗方法被称为随机抽样的原因 最简单 最基本 最重要的一个概 率分布是 0 1 上的均匀分布 或称矩形分布 随机数就是具有这种均匀分布的随机变量 随机数序列就是具有这种分布的总体的一个简单子样 也就是一个具有这种分布的相互独 立的随机变数序列 产生随机数的问题 就是从这个分布的抽样问题 在计算机上 可以 用物理方法产生随机数 但价格昂贵 不能重复 使用不便 另一种方法是用数学递推公 式产生 这样产生的序列 与真正的随机数序列不同 所以称为伪随机数 或伪随机数序 列 不过 经过多种统计检验表明 它与真正的随机数 或随机数序列具有相近的性质 因此可把它作为真正的随机数来使用 由已知分布随机抽样有各种方法 与从 0 1 上均匀 分布抽样不同 这些方法都是借助于随机序列来实现的 也就是说 都是以产生随机数为 前提的 由此可见 随机数是我们实现蒙特卡罗模拟的基本工具 3 1 3 3 建立各种估计量 一般说来 构造了概率模型并能从中抽样后 即实现模拟实验后 我们就要确定一个 随机变量 作为所要求的问题的解 我们称它为无偏估计 建立各种估计量 相当于对模 拟实验的结果进行考察和登记 从中得到问题的解 3 3 4 适用范围 3 1 蒙特卡罗的特点是模拟次数越多 计算结果的可靠性越大 特别适用于在计算机上 对大型项目 新产品项目和其他含有大量不确定因素的复杂决策系统进行风险模拟分析 2 蒙特卡罗模拟法不可能使计算结果发生实质性变化 但是可以给出计算结果的概率 分布 便于预测达到预期目标的可能性 3 1 5 蒙特卡洛方法的工作过程 在解决实际问题的时候应用蒙特 卡罗方法主要有两部分工作 用蒙特卡罗方法模拟某一过程时 需要产生各种概率分布的随机变量 用统计方法把模型的数字特征估计出来 从而得到实际问题的数值解 3 1 6 蒙特卡洛方法在数学中的应用 通常蒙特 卡罗方法通过构造符合一定规则的随机数来解决数学上的各种问题 对于那 浙江科技学院本科毕业设计 论文 11 些由于计算过于复杂而难以得到解析解或者根本没有解析解的问题 蒙特 卡罗方法是一种 有效的求出数值解的方法 一般蒙特 卡罗方法在数学中最常见的应用就是蒙特 卡罗积分 蒲丰投针实验 1 在用传统方法难以解决的问题中 有很大一部分可以用概率模型进行描述 由于这类 模型含有不确定的随机因素 分析起来通常比确定性的模型困难 有的模型难以作定量分 析 得不到解析的结果 或者是虽有解析结果 但计算代价太大以至不能使用 在这种情 况下 可以考虑采用 Monte Carlo 方法 下面通过例子简单介绍 Monte Carlo 方法的基本 思想 Monte Carlo 方法是计算机模拟的基础 它的名字来源于世界著名的赌城 摩纳哥 的蒙特卡洛 其历史起源于 1777 年法国科学家蒲丰提出的一种计算圆周率 p 的方法 随机投针法 即著名的蒲丰投针问题 这一方法的步骤是 1 取一张白纸 在上面画上许多条间距为 d 的平行线 2 取一根长度为的针 随机地向画有平行直线的纸上掷n次 观察针与直线 L Ld 相交的次数 记为 m 3 计算针与直线相交的概率 4 经统计实验估计出概率 m P n 1 蒲丰投针试验 在平面上画有等距离为的一些平行线 向平面上随 0d d 机投一长为的针 求针与平行线相交的概率 L Ld P A 解 若以 M 表示针的中点 以表示 M 距最近平行线的距离 表示针与平行线的交x 角 则针与平行线相交的充分必要条件是满足 x 0sin 2 0 xL 于是蒲丰投针试验就相当于向平面区域投点的几何型 0 0 2 d Gxx 随机试验 如图 浙江科技学院本科毕业设计 论文 12 图图 3 1 蒲丰投针简图 建立直角坐标系 上述条件在坐标系下将是曲线所围成的曲边梯形区域 见图 j x 8 1 2 由几何概率知 0 1 sin 22 2 d gl p d Gd 的面积 的面积 4 经统计实验估计出概率由 式即 m P n 2 ml nd Monte Carlo 方法的基本思想是首先建立一个概率模型 使所求问题的解正好是该模 型的参数或其他有关的特征量 然后通过模拟统计试验 即多次随机抽样试验 确定 m 和 n 统计出某事件发生的百分比 只要试验次数很大 该百分比便近似于事件发生的概 率 这实际上就是概率的统计定义 利用建立的概率模型 求出要估计的参数 蒙特卡洛 方法属于试验数学的一个分支 提示 设 x 是一个随机变量 它服从区间是的均匀分布 同理 j 是一个随机 02d 变量 它服从区间上的均匀分布 按照某种抽样法 产生随机变量的可能取值 例 0p 如进行 n 次抽样 得到样本值统计出满足不等式的次数 1 2 ii xjin sin 2 ii d x 从而可以计算出 p 的估计值 m mn pm n 同一原理 利用蒙特卡洛算法计算的值 基本思路 利用求单位圆的 1 4 的面积来求得从而得到 单位圆的 1 4 面积是4 一个扇形 它是边长为 1 单位正方形的一部分 只要能求出扇行面积在正方形面积中 1 SS 占的比例就立即能得到 从而得到的值 怎样求出扇形面积在正方形面积 1 KSS 1 S 浙江科技学院本科毕业设计 论文 13 中占的比例 K 呢 一个办法是在正方形中随机投入很多点 使所投的点落在正方形中每一 个位置的机会相等看其中有多少个点落在扇形内 将落在扇形内的点数 m 与所投点的总数 n 的比作为 k 的近似值 P 落在扇形内的充要条件是 m n 22 1xy 取不同的做上面的实验 将所得的的近似值记录下来 并与精确值比较 n 3 2 数值积分法 数值积分法 是从近似计算的角度 采用某种数值过程来求出定积分的近似值 我们知道 如果函数 f x 在区间 a b 上连续 且原函数为 则可用牛顿 莱 F x 布尼兹公式 aFbFdxxf b a 来求得定积分 然而 对有些函数来说 找到原函数往往很困难 有些原函数不能用初等 函数表示出来 例如 a 为有限数 dxe a x 0 2 还有的被积函数尽管能用初等函数的有限形式表示出来 但表达式过于复杂 也不便 使用 特别是在实际问题中 更多的函数是用表格或图形表示的 对这种函数 更无法用 牛顿 莱布尼兹公式求积分 存在很大局限性 因此 有必要研究用数值方法求定积分 的问题 这种数值积分方法也是微分方程 积分方程数值解法的基础 数值积分的基本思想是构造一个简单函数 例如多项式 来近似代替被积分函 n Px 数 然后通过求 求得的近似值 f x b a n dxxP b a dxxf 常用的定积分的近似计算方法有梯形公式法和辛普森公式法 3 2 1 常用求积公式 设 b a dxxfI 3 1 插值型求积公式就是构造插值多项式 使 n Px b a n dxxPII 称 3 1 为插值求积公式 4 浙江科技学院本科毕业设计 论文 14 由拉格朗日插值公式 0 k n k kn xfxlxP 代入 3 1 则有 00 k n k b a kk n k b a k xfdxxldxxfxlI 记 0 k n k k b a kk xfIdxxl 其中为求积系数 为求积节点 k k x 3 2 1 1 梯形公式 构造以 a b 为结点的线性插值多项式 1 bf ab ax af ba bx xP 2 1 2 1 2 1 22 1 bfafabab ab bf ba ba af dxax ab bf dxbx ba af dxbf ab ax af ba bx dxxPT b a b a b a b a 3 2 称 3 2 为梯形公式 4 梯形求积公式的基本思想是 在每个小区间上采用线性函数 近似被积函数 f x 3 2 1 2 辛普森公式 如果我们不用线性函数而改用二次多项式来近似 那么就可以得到另一种在实 f x 际中常用的辛普森求积公式 以 a b 为节点 构造二次插值多项式 2 ab c 2 bf cbab cxax cf bcac bxax af baca bxcx xP 浙江科技学院本科毕业设计 论文 15 210 2 bfcfaf dx cbab cxax bf dx bcac bxax cfdx baca bxcx af dxxPS b a b a b a b a 其中 6 1 2 3 1 1 1 1 22 33 2 0 ab abbccb ab ab baca dxbcxcbx baca dxbxcx baca dx baca bxcx b a b a b a 类似地 可求得 6 4 1 abdx bcac bxax b a 6 1 2 abdx cbab cxax b a 得三点公式 4 6 bfcfaf ab S 3 3 称 3 3 为辛普森生 Simpson 求积公式 4 我们也可以类似地构造五点 七点 n 点求积公式 在实际运算中 以梯形和辛 普森 Simpson 公式为常用 3 2 2 复化求积公式 如果积分区间比较大 直接地使用上述求积公式 精度难以保证 在高等教学的学习 中 曾经介绍过梯形法求积公式 它的几何意义是将积分区间 a b 分成 n 个小区间 对 用分段线性插值 然后积分 类似地 也可以对用分段抛物插值 f x f x 浙江科技学院本科毕业设计 论文 16 通常采取的办法是 1 等分求积区间 比如取步长 分 a b 为 n 等分 分点为 ba h n 0 0 1 2 k xxkh kn 2 在区间 上使用以上求积公式求得 1 kk xx k I 3 取和值 作为整个区间上的积分近似值 1 0 n k k II 这种求积方法称为复化求积方法 2 3 2 2 1 复化梯型公式 2 由 3 2 2 1 kkk xfxf h I 3 4 2 2 2 1 1 1 1 0 1 0 bfxfaf h xfxf h IT n k k kk n k n k kn 3 2 2 2 复化辛普森公式 2 将对分 中点记为 1 kk xx 2 1 2 1 kk k xx x 由 3 3 4 6 1 2 1 k k kk xfxfxf h I 3 5 4 6 1 1 0 1 0 2 1 k k k n k n k kn xfxfxf h IS 为编程序方便 我们将辛普森公式 3 5 写成 n k k k xfxf bfaf n ab S 1 2 2 3 2 1 3 2 3 变步长梯形方法 3 前面介绍的复化求积方法对提高精度是行之有效的 但是使用复化求积公式之前须给 出合适的步长 步长取得太大精度难以保证 步长太小则会导致计算量的增加 而事先给 出一个合适的步长往往又是十分困难的 浙江科技学院本科毕业设计 论文 17 下面介绍一种变步长积分法 其基本思想是将区间逐次对分进行计算 用前后两次计 算的结果进行估计 若合乎精度要求 就停止计算 否则再次对分 重复以上计算过程 直至达到精度要求为止 设将区间 a b n 等分 共有 n 1 个分点 按复化梯形公式计算 需要计算 n 1 个 n T 的值 如果将求积区间再次对分 若仍然直接用复化梯形公式计算二分后的积分值 f x 则需要计算 2n 1 个的值 注意到 的全部分点中有 n 1 个是二分前原有的 2n T f x 2n T 点 重复计算这些 老分点 上的函数值显然是个浪费 因此 我们将二分前后两个积分 值联系起来加以考虑 注意到每个小区间经过二分再增加一个新分点后 用 1 kk xx 2 1 k x 复化梯形公式求得该区间上的积分值为 2 4 1 2 1 k k kk xfxfxf h I 因此有 1 0 1 0 1 0 1 1 1 0 2 22 1 2 4 2 4 2 1 2 1 2 1 n k k n n k k n k kk k k k n k n xf h T xf h xfxf h xfxfxf h T 3 6 3 6 式的前一项是二分前的积分值 在求时可作为已知值使用 而它的后一 n T 2n T 项只涉及二分时新增加的分点 所要计算 f 值的次数为 n 可见递推公式 3 6 由于 2 1 k x 避免了老节点的重复计算 而使计算量减少了一半 3 2 4 求积公式的误差 现在考虑求积公式的误差 因为是的 n 次插值多项式 所以当本身 n Px f x f x 就是次数不超过 n 的多项式时 这时 在不考虑舍入误差时 求积公式 n f xPx 是精确成立的 我们来看一下它的舍入误差 此时 0 k n k k b a n b a xfdxxPdxxfI 浙江科技学院本科毕业设计 论文 18 取 则 1f x ab n k k 0 若的舍入误差小于 则 k f x 0 0 0 abxfxf xfxfII kk n k k k n k kk n k k 可见 舍入误差对数值积分的影响不大 因此 在讨论求积公式的误差时 主要考虑 截断误差 3 2 4 1 梯形公式的截断误差 3 在求积区间 a b 上 设 M2 0 常数 在 a 点附近将展成泰勒 2 Mxf f x 级数 baaxfaxafafxf 2 2 1 32 2 6 1 2 1 2 1 abfafabafab dxaxfaxafafdxxfI b a b a 另一方面 注意到 2 2 1 abfabafafbf 32 2 4 1 2 1 2 1 2 2 2 abfabafafab abfabafaf ab bfaf ab T 3 12 1 abfTI 梯形公式的截断误差为 3 2 12 ab M TI 若取 对每个小区间用梯形公式 n ab h 1 kk xx 浙江科技学院本科毕业设计 论文 19 2 2 1 0 3 1 0 1 0 12 12 1 12 hf ab h n ab f h f IITI n k n k n k kkn 复化梯形公式的截断误差限为 2 2 3 12 hM ab TI n 3 2 4 2 辛普森公式的截断误差 3 假设在区间 a b 上 M4 0 常数 在附近 将展成 4 4 Mxf 2 ba c f x 泰勒级数 bacxfcxcf cxcfcxcfcfxf 4 4 3 2 4 1 3 1 2 1 5 4 3 4 4 3 2 2 60 1 3 1 4 1 3 1 2 1 cb fabcfabcf dxcxfcxcf cxcfcxcfcfdxxfI b a b c 另一方面 4 4 22 2 24 1 2 6 1 22 2 ab f ab cf abcfab cfcfaf 4 4 22 2 4 1 2 3 1 22 2 ab f ab cf abcfab cfcfbf 从而 浙江科技学院本科毕业设计 论文 20 5 4 3 4 4 2 2 36 1 2 3 1 2 12 1 2 4 2 6 4 6 ab f ab cfabcf ab f ab cfcfcf ab bfcfaf ab S 5 4 5 4 5 4 2 90 1 2 36 1 2 60 1 ab f ab f ab fSI 截断误差限为 5 4 290 abM SI 取 得复化辛普森公式的截断误差 n ab h 4 4 5 4 1 1 2 180 2 90 1 h f ab h f IISI n k n k kkn 截断误差限为 4 4 2180 h M ab SI n 3 2 5 应用 2 选取不同的 n 由公式 1 2 0 4 1 dx x 分别利用上述的梯形求积公式和辛普森求积公式计算积分的近似值 从而 1 2 0 4 1 dx x 求得的近似值 并给出其误差 浙江科技学院本科毕业设计 论文 21 解 由编程计算得到如下数据表 表表 3 1 nquanError 43 13117647058823529410 0104161830015579443 83 13898849449108900930 0026041590987042291 123 14043524684685065150 0011574067429425870 163 14094161204138889460 0006510415484043438 203 14117598695412846340 0004166666356647750 243 14130330174832387340 0002893518414693650 283 14138006855989701290 0002125850298962256 323 14142989317497443320 0001627604148188053 nsimpsonError 43 14159250245870691442134859298445070663490136635031 511310863240412947902950521775622680330 10 7 83 14159265122482218974665307283818498578450284518262 364971048715990310441317898412666554193 10 9 123 14159265338214828075850499177714648429147079188682 076449577041383915023563999056986074883 10 10 163 14159265355283627939271421231154251794086237072773 69569590699291709679603662563070286474 10 11 203 14159265358010514912785425032645269144165440617959 6880893347891329530501927555149931956 10 12 243 14159265358654871162257301393977061438013383458463 2445268400703693397322698170355647905 10 12 283 14159265358850655731166078997881603197253338791501 2866811509825933006868522246360114601 10 12 323 14159265358921578178840022303357898154107427521795 774566742431602459239026560951241572 10 13 计算积分 1 2 0 4 1 Idx x 显然有 1 0 4arct3 1415926Igx 由计算机运行结果可知当取 n 8 用复化梯形公式有 8 3 13898849T 误差为 0 00260416 取 n 4 用辛普森公式 4 3 14159250S 误差为 7 1 511311 10 比较与两个结果 它们却需要调用 f 9 次 工作量基本相同 但精度差别却很 8 T 4 S 大 这项计算结果表明 复化辛普森公式是一种精度较高的求积公式 3 对上题中误差的实际观察 选取等 观察值的增1000 10000 100000 n n 加所导致的值的变化情况 直到的增加所导致的的变化小于给定误差界 sns 比较同一个值下梯形公式和辛普森公式计算结果的差别 对两种方法的精度差n 别获得一个感性认识 解 编程由 mathematica 运行得到的结果如下 表表 3 2 浙江科技学院本科毕业设计 论文 22 nquanError 10003 14159248692312657181 666666666667 10 7 20003 14159261192312657184 16666666667 10 8 30003 14159263507127471991 85185185185 10 8 40003 14159264317312657181 04166666667 10 8 50003 14159264692312657186 6666666667 10 9 60003 14159264896016360884 6296296296 10 9 70003 14159265018843269423 4013605442 10 9 80003 14159265098562657182 6041666667 10 9 90003 14159265153218006972 0576131687 10 9 100003 14159265192312657181 6666666667 10 9 110003 14159265221238277011 3774104683 10 9 120003 14159265243238583111 1574074074 10 9 130003 14159265260359994469 861932939 10 10 140003 14159265273945310248 503401361 10 10 150003 14159265284905249777 407407407 10 10 160003 14159265293875157186 510416667 10 10 170003 14159265301309196975 767012687 10 10 180003 14159265307538994635 144032922 10 10 190003 14159265312811272144 616805171 10 10 200003 14159265317312657184 166666667 10 10 nsimpsonError 10003 14159265358979323846202334359696351608101905579446 200396825393681161503435807 10 22 20003 14159265358979323846263369515946320196454109035369 6881200396822326283090215 10 24 30003 14159265358979323846264253274495893402068392538718 505345439501764854739880 10 25 40003 14159265358979323846264323190262726415778671733971 513768756200393826820354 10 25 50003 14159265358979323846264334