数学实验第二章

第二章方程求解 本章的主要内容 1 迭代算法及编程 2 迭代算法的问题及改进 3 迭代算法的MATLAB实现 4 与解方程相关的数学建模 2010 10 21 2 1导言 若干世纪以来 工程师和数学家花了大量的时间用于求解方程 组 研究了各种各样的方程求解方法 这一讨论和求解推动了数学学科和科学技术的发展 本章主要目的是介绍方程求解的实用方法 利用现代计算技术和计算机的综合优势求出符合精度要求的那些解 2 2引例 求方程的非负实根 构造函数 由于连续 并且在和处函数值符号相反 可知f x 在区间 0 1 内必有零点 也就是方程在 0 1 内必然存在根 进一步 以方程的两边绘出函数图像 root 取一个根的估计值 代入方程两边 考察两边的接近程度 验证这个估计值是否为方程的根 猜测x 0 代入方程右边为1 3 左右差为0 3333 再取右边表达式的值1 3作为新的估计值 代入方程右边 左边为1 3 右边为0 4652 左右差为0 1319 再取右边表达式的值0 4652作为新的x估计值 代入方程右边 左边为0 4652 右边为0 5308 左右差为0 0656 如果 不断地用上一次右边表达式ex 3的值作为方程根的近似值 而方程两边不断接近 那末 就可以考虑用迭代的方法来求解 参照 参照 x 在计算机中 从当前的问题出发 我们只需根据前一项算出后一项 然后计算两项之差 并与精度要求比较即可 因此 只设置一个变量 首先 x 前一项 然后 利用x算出后一项 从数学角度来看 由方程可以给出一个递推数列 如果数列收敛 或前后两项无限接近 则有 即极限值就是方程的根 x 如果 x 的绝对值小于10 5 前项 后项 则可把看成是root的近似值 x 如果 x 的绝对值 10 5 前项 后项 则重新赋值 重复上述过程 ex 3 ex 3 ex 3 第二个根 可用x log 3 x 一般地 应考虑 的反函数 由表2 1中结果可知 随着迭代次数的增加 方程左右两边的差距在不断缩小 到第9次 小数点后前二位有效数字相同 继续下去 可以得到方程满足精度要求的近似根 1 例中 如果近似根的精度要求为小数点后10位有效数值 它是多少 2 迭代初始值xo分别取不同值 例如0 1 2 迭代结果有何区别 3 方程的非负实根是否唯一 如果不止一个如何求出其他根 以后 可通过理论分析 或画图来确定其他根的大致位置 然后通过迭代法尝试求解 想 2 3方程和方程组方程和方程组都有线性和非线性之分 其中单个的线性方程太简单了 这里略去 下面我们从线性方程组开始 编程回答 2 3 1线性方程组线性方程组 可以写成矩阵的形式 由线性代数知识可知 关于其解的理论如下 解可能出现三种情形 无解 有唯一解和有无穷多组解 这主要取决于系数矩阵A的秩与增广矩阵 A b 的秩是否相等 当两个秩相等时 与变量个数的关系 具体地 若则无解 为变量个数 则 有唯一一组解 若 若 则 有无穷多组解 求解方法可以划分为两类 直接消去方法 这主要是指高斯消去法 迭代数值解法 关于线性方程组的消去法和迭代算法 简介如下 1 直接消去法理论上 经过有限次算术运算能够求出方程组的精确解 实际上由于计算存在舍入误差 得到的只是近似解 高斯消去法可以用LU分解来表示 如果方程组的系数矩阵A满足一定条件 则A可以分解为一个下三角阵L与一个上三角阵U的乘积 A LU 于是由得到 由于上三角阵与下三角阵的逆矩阵很容易求得 因此 这是一种很好的思想 例如 A 0122103 00210101 01020002 2 迭代法常用的有雅比迭代法和高斯 塞德尔迭代法等 将同解的变形为 可以证明 在B满足谱半径 当时 收敛 其极限就是原方程组的解 x0 bx f x 使得x bx f 谱半径 矩阵A 其各个特征值的最大绝对值称为A的谱半径 取根的估计值 代入上式左右两端 考察左右两端是否相等 接近程度是否符合精度要求 如果不符合要求 则将bx f的值作为新的 根 估计值 又代入上式左右两端进行比较 这一往复的过程可以简化为迭代公式 变量x 例如 雅可比迭代法 由得到 其中 高斯 塞德尔迭代法 将上述矩阵分解的迭代公式写成 其中 设A矩阵可分解为A L D U其中D为对角矩阵 L为下三角矩阵 U为上三角矩阵 用两种不同的迭代算法求解以下线性方程组 以w 11 0251 051 075 作为方程解的初始估计值 5 100 15 10 0 15 1 00 15 x1 x2 x3 x4 4 3 3 8 3 1 4 9 做 2 3 2非线性方程从中学到大学 对数学课程来讲 解方程可以说是经常接触的事情 早在19世纪初 伽罗瓦等人就已证明了高于四次的代数方程不存在类似于二次方程的求根公式 因此 在数学的教科书上 只介绍了一些特殊的代数方程的解法 我们今天针对一般非线性方程 提出一套行之有效的求解方法 考虑一般方程可能难以甚至根本无法找到其精确解 因此我们只好去找近似解 但应该有一定的精度要求 本节的目的就在于探讨能够求解任意一个方程全部根或者指定范围的部分近似根的方法 做 1 对于方程tanx x 用迭代算法求其非负实数根 这个形式x tanx本身就给出了一个迭代公式 还可以考虑不同的迭代公式 因为tan x sin 2 x 1 cos 2 x 所以有迭代公式x sin 2 x 1 cos 2 x 编程x 1 根的初始估计值 弧度值n 0 记录迭代次数whileabs x tan x 10 5 x tan x n n 1 endx n 2 想办法用其他方法求解tanx x 如几何的方法 并对所得到的结果进行分析 3 考虑解方程 显然无根 同解地变形为x x exp tan x 2 以此为迭代公式编程 可求出根 增根 对以上三个小题的讨论可以发现 用迭代方法求解方程 在方程有根的情况下 迭代算法可能找得到根 也可能找不到根 而在方程无解的情况下 可能找出了增根 因此 对于方程求解 常常需要结合具体方程 函数 的性质 从分析的和数值的乃至几何图形的角度上加以考虑 需要多种不同方法 对于方程f x 0的求解 首先可以经过简单变形把它表达为x g x 例如x x f x 记g x x f x 则求解f x 0就等于求解x g x 也就是求函数g x 的不动点 因此 可以尝试用迭代的方法来解决 称此法为方程求解的不动点法 其迭代过程如下 xn 1 g xn n 0 1 其中x0为迭代初值 通常应取为不动点的估计值 2 3 3一般的迭代算法 下例可以说明不同的迭代公式产生不同的迭代结果 例2 1用迭代方法求解方程 解第一步同解地构造迭代公式 对方程 可以给出多个同解的变形 如 它们的迭代效果是不同的 比如x 1作为初始值 即 求一个x使得abs x 3 x 2 x 1 工程设计所要求的很小的正数 第二步选择迭代初值x 1 对以上三个公式分别进行迭代 和的迭代结果如下表 可以看到它们是收敛的 极限近似为1 8393 表2 2迭代产生的数列 初值 序号n1234567891 44221 65371 75321 79951 82091 83081 83541 83751 83853 00001 44442 17161 67251 95541 77301 88221 81361 8554 通过改变初值 比如或 以及其他的诸多初值 和的迭代结果都显示出良好的收敛性 在迭代到第23次以后 结果几乎都为1 8393 而的迭代结果始终令人失望地无限增大 编程x 1 n 0 whileabs x 3 x 2 x 1 10 5 x x 3 x 2 1 n n 1 ifabs x 3 x 2 x 1 9999breakendendx n 结果 仅迭代10次 x inf 迭代失败 由此可知 迭代公式的选取是很重要的 关于第一个迭代公式的编程 第三步迭代失败的改进 加速迭代收敛方法 设采用迭代公式x 失败 考虑与x的加权平均 进行迭代 为参数 x h x 取值为介于0 1之间的数 h x 介于x 0 5 的选择应有利于加速迭代 之间 数学工作者通过对大量成功和失败的迭代案例进行分析 认识到 在满足条件下 取 可以取得一个合适的 这个条件的意义 在a附近 h x 的变化较小 详见 数值分析 教材 考虑到a是未知的 在实际迭代过程中采用xn代替a 因此得迭代公式 对于 进行改进 产生新的迭代公式 加速迭代公式如下 表2 3迭代函数产生的数列 序号n123456789初值0 1 0000 0 50000 6667 1 1481 0 63710 0516 0 9109 0 40512 3240初值1006 81364 71103 34612 49962 03901 86581 83991 83931 8393 改进后的迭代法编程x 1 迭代初值n 0 记录迭代次数whileabs x 3 x 2 x 1 10 5 迭代条件x 2 x 3 x 2 1 3 x 2 2 x 1 n n 1 endx 数值解n 迭代次数 运行的结果 x 1 8393n 13 a b 2 点的迭代法迭代法是构造与方程等价的方程得到迭代函数然后求迭代函数的不动点 注意到形式不唯一 其迭代差异可能很大 如何构造函数 能保证迭代算法的收敛性 对此 有如下结论 定理设在 a b 上连续 且满足 A 对任意的 B 存在常数使得对于任意的 则对任意初值 迭代过程产生的序列收敛在上的不动点 解释一下定理的几何意义由A 在方框内 由B 可以推出即曲线上任意两点的割线斜率 其绝对值不大于1 显然 这是一个很强的充分条件 返回 a b 通过上例可知 求解方程的迭代函数构造形成多样 不同迭代函数的收敛性及对初值的敏感性可能千差万别 因此在选择迭代函数及初值上 值得我们多多注意 在遇到迭代失败的时候 可以利用如上介绍的加速迭代收敛方法加以解决 1 方程有多个实根 能否将其全部找出来 用什么方法 随后介绍的图解法2 以和为例 考察对于迭代收敛的情况 以上加速迭代收敛方法是否能够起到加速收敛的作用 给出充分体现加速迭代的例子 迭代函数 想 以上分析可以引进到求解方程组的情形 2 3 4方程求根系列方法 考虑求方程的解 1 熟知的因式分解法因式分解是我们常用的方法 它的基本原理是 如果则必有或 这个方法的关键在于分解因式 包括对多项式函数 三角函数和指数函数等分解 例2 2求解以下方程 1 2 解1 通分得解 分解因式得 解出 分解因式得 因 所以 这等价于 即 得到其根为 2 图形放大法由于计算机的广泛应用 可以非常方便地作出函数的图形 曲线 找出曲线与x轴的交点的横坐标值 就可求出的近似根 即的零点 这些值虽然粗糙但直观 方程有多少个根 在何范围 一目了然 并且还可以借助于使用图形局部放大功能 将根定位得更加准确 因此如果你拥有一台计算机 那么图形的方法将是分析方程求根的性态最简捷的方法 不过 不要总是想要得到根的精确值 例2 3求方程所有的根及大致分布范围 x linspace 1 54375 1 54365 100 y x 5 2 x 2 4 plot x y 通过观察 不断提高精度 用图形放大法求方程的根 注意观察结果 它有多少个根 将上式变形为 可以画出两条函数曲线与 两条曲线的交点的横坐标x就是原方程的根 做 提示 3 数值迭代逼近法 利用图形的方法或连续函数的零点存在性定理 可以推知在某一区间内有根 下面就基于此来讨论求根的数值方法 关于区间的迭代有对分法 黄金分割法等 关于点迭代有简单迭代法 单点割线 两点割线法以及使用导数值的牛顿法等等 x 0 1 pi pi 100 4 pi y1 sin x y2 1 x plot x y1 r x y2 g 1 区间迭代法对分法是应用零点存在性定理 每次将区间压缩一半且其中一个区间至少包含一个根 逐步缩短区间 直至最终区间长度满足一定的精度要求为止 黄金分割法与对分法本质上一致 只不过每次压缩区间的比例不是一半 而是压缩比例为0 618 黄金分割比例 例如 解方程exp x 8 sqrt x 1 log x 2 其中x属于 1 2 令f x exp x 8 sqrt x 1 log x 2 有f 1 0 f 2 0 a 1 b 2 whileb a 10 5 ifexp a b 2 8 sqrt a b 2 1 log a b 2 2 0a a b 2 elsesprintf f s a b 2 是一个根 endendsprintf f s a b 2 是一个精确到万分位的根 2 点的迭代法简单迭代法是基于构造与方程等价的方程得到迭代函数然后求迭代函数的不动点 注意到形式不唯一 其迭代差异可能很大 如何构造函数 能保证迭代算法的收敛性 对此 有如下结论 定理设在 a b 上连续 且满足 A 对任意的 B 存在常数使得对于任意的 则对任意初值 迭代过程产生的序列收敛在上的不动点 解释一下定理的几何意义由A 在方框内 a b 下面几个迭代方法是根据的几何性质与解析性质 从不同角度构造出的迭代公式 单点割线法是给定初值x0 点p0 及当前值xn 点pn 将曲线y f x 上当前点和初始点构成的割线 代替曲线 与x轴的交点的横坐标 作为本次迭代的输出也为下次迭代的输入 其迭代公式为 由B 可以推出即任意两点的割线斜率 其绝对值不大于1显然 这是一个很强的充分条件 这个公式来源于 两点割线法是给定初值x0及x1 两个点p0 p1 将曲线上最新的相邻 按迭代先后为序 两点的割线 代替曲线 与x轴的交点的横坐标 作为迭代输出 其迭代公式如下 x0 x1 x2 x0 x1 x2 两点割线法的程序示例 求方程x 4 x 3 x 2 2 0的根1 作图观察根的大致范围x 2 0 1 2 y x 4 x 3 x 2 2 plot x y 2 编程计算近似根 x0 1 y0 x0 4 x0 3 x0 2 2 x1 1 5 y1 x1 4 x1 3 x1 2 2 x2 x1 y1 y1 y0 x1 x0 n 0 whileabs x2 4 x2 3 x2 2 2 10 5 x0 x1 y0 x0 4 x0 3 x0 2 2 x1 x2 y1 x1 4 x1 3 x1 2 2 x2 x1 y1 y1 y0 x1 x0 n n 1 endsprintf f为一个近似根迭代 d次数 x2 n f x x0 p0 x1 p1 x2 牛顿法是给定值x0 点p0 将过曲线上当前点的切线与X轴的交点的横坐标作为本次迭代的输出 其迭代公式如下 设有方程f x 0 由f xn xn 1 xn f xn 有迭代公式xn 1 xn f xn f xn 此图的斜率是负数 牛顿迭代法的例题 解方程f x 0这里f x 5 x 5 x 4 x 2 3 x 200 因此f x 25 x 4 4 x 3 2 x 3 迭代公式为编程 n 0 x 1 y 5 x 5 x 4 x 2 3 x 200 y1 25 x 4 4 x 3 2 x 3 whileabs y 10 5 x x y y1 y 5 x 5 x 4 x 2 3 x 200 y1 25 x 4 4 x 3 2 x 3 n n 1 endsprintf d是方程的一个近似根 d为迭代次数 x n 1 关于方程组的数值解 简单迭代法的迭代函数构造与方程情形一致 2 迭代收敛较慢或不收敛时 可以考虑采用前面介绍的加速迭代收敛方法 注意 2 3 5非线性方程组求解非线性方程组的一般形式如下 对于非线性方程组 也适用迭代解法 为表述简单起见 写成向量形式 其中向量为n个构成的n维向量 如果方程组中所有函数都是变量的线性函数 则称之为线性方程组 它是一类特殊的方程组 有其独特的求解方法 n元非线性函数 sin x y 2 ln z 7 3 x 2 y z 3 1 x y z 5 下面介绍牛顿迭代法求解非线性方程组的数值解法 考虑非线性方程组的向量形式F x 0 其中x和F x 都为m维向量或向量函数 若用x n 表示第n次迭代产生的输出 向量 通过将F x 在当前点x n 处展开成一阶泰勒展式 并用x n 1 代替x 就得到求解方程组的牛顿迭代公式 上例 121 取a F a sin 1 2 2 ln 1 7 3 1 2 2 1 3 1 1 2 1 5 F a cos x 2 y1 z 3 2 y ln 2 3 z 2 111 a 0 540341 32 77263 111 记a x 0 a 111 f1 sin a 1 a 2 2 log a 3 7 f2 3 a 1 2 a 2 a 3 3 1 f3 a 1 a 2 a 3 5 n 0 whileabs f1 abs f2 abs f3 10 5 f1 sin a 1 a 2 2 log a 3 7 f2 3 a 1 2 a 2 a 3 3 1 f3 a 1 a 2 a 3 5 f f1f2f3 g cos a 1 2 a 2 1 a 3 3 2 a 2 log 2 3 a 3 2 111 a a inv g f n n 1 ifn 1000 abs f1 abs f2 abs f3 10 5breakendendan 2 牛顿法中 每项次迭代都需要计算雅可比矩阵的逆矩阵 当函数较为复杂或阶数较高时 给计算带来相当不便 因此出现了各种基于牛顿法的修正方法 如拟牛顿法等等 2 4数值求解软件 MATLAB可求方程及方程组的解析解或数值解 2 4 1任意函数方程与线性方程组 MATLAB命令输入格式 或 事实上 这是求解方程的牛顿迭代公式的推广 注意 1 此处要求在处的雅可比矩阵可逆 其中eqni表示第i个方程 vari表示第i个变量 1 一般方程例如 求解方程 输入 输出 又如求解方程 输入 输出 0 注 该方程应该有无究多个解 这里只得到其中的一个 注意 如果不能求得精确的符号解 可以计算可变精度的数值解 例如 求解方程 输入 输出 0 51775736368245829832278747416629 2 多项式方程 除了用上面求解一般方程的方法外 还可以直接用求解多项式方程的MATLAB函数roots p 其中p是多项式的系数按降幂排列所形成的n 1维向量 它能够给出全部根 包含重根 例如 求解多项式方程 输入 roots p 输出 1 2131 0 9017 0 5753i 0 9017 0 5753i 0 2694 0 9406i 0 2694 0 9406i 0 4168 0 8419i 0 4168 0 8419i 0 8608 0 3344i 0 8608 0 3344i 注意 也可以用求解 精度高达小数点后32位 3 线性方程组 除了使用MATLAB函数solve以外 还可以用其他的MATLAB命令 如果将线性方程组写成矩阵形式AX b 就可以考虑用几种形式之一求解 其中inv A 表示A的逆矩阵 因此要求A为方阵且可逆 pinv A 表示A的广义逆矩阵 A可以为任意矩阵 以上MATLAB函数均可以对任意的线性方程组求解 不管有解 无解 有一个还是有无穷多 它们有何区别 想 以A b为例1 当AX b有唯一解时 给出该唯一解 2 当其有无穷多解时 给出其中零元素最多的一组解 3 当其无解时 给出一个最小二乘 广义 解 对于pinv A b当AX b有无穷多解时 给出其中一个最小范数解 长度最小的向量 其他两种情形与A b相同 Linsolve A b 针对线性齐次方程组 等价于A b 算法不相同 相同 考虑以下线性方程组AX b的求解 A 410 1 15 22 3 b 6 14 3 其中hilb n 表示n阶希尔伯特矩阵 或 1 1 3 1 n 1 21 31 4 1 n 1 1 31 41 5 1 n 2 1 n1 n 1 1 n 2 1 2n 1 sum 分析 1 此例中说明方程组有无穷多解 2 此例中说明方程组有唯一解 3 此例是希尔伯特方程 说明方程组有唯一解 2 4 2非线性方程组 一些非线性方程组仍然可以用solve 函数进行求解 一般方程式组有唯一解 解 例如 输入 输出 也可以用以下fsolve 进行求解 输入格式为 其中FUN表示函数M文件 X0表示初始向量 例如求解下列方程组 1 首先建立方程组的函数M文件 2 执行经下程序 functioneq nx x globalnumber number number 1 eq 1 sin x 1 x 2 2 log x 3 7 eq 2 3 x 1 2 x 2 x 3 3 1 eq 3 x 1 x 2 x 3 5 globalnumber number 0 申明number是全局变量 每次迭代number均自加 定义全局变量 赋初值 3 计算结果 0 59912 39592 0050 29 其中迭代步骤为29次 2 5范例 波音公司飞机最佳定价策略 全球最大的飞机制造公司波音公司自1955年推出波音707开始 成功地开发了一系列的喷气式客机 问题 讨论该公司对一种新型客机最优定价策略的数学模型 y fsolve nx 111 调用求解函数 111 为解的估计值 解方程组的迭代步骤体现在fsolve中 方程组写在nx函数中 返回 2 5 1问题分析 定价策略涉及到诸葛亮多因素 这里考虑以下主要因素 价格 竞争对手的行为 出售客机的数量 波音公司的客机制造量 制造成成本 波音公司的市场占有率等等恩素 2 5 2假设及模型 价格记为p 根据实际情况 对于民航飞机制造商 能够与波音公司杭衡的竞争对手只有一个 因此他们可以在价格上达成一致 具体假设事下 1 型号 为了研究方便 假设只有一种型号飞机 2 销售量 其销售量只受飞机价格p的影响 预测以此价格出售 该型号飞机全球销售量为N N应该受到诸多因素的影响 假设其中价格是最主要的因素 根据市场历史的销售规律和需求曲线 假设该公司销售部门预测得到 3 市场占有率 既然在价格上达成一致 即价格的变化是同步的 因此 不同定价不会影响波音公司的市场占有率 因此场占有率是常数 记为h 4 制造成数量 假设制造量等于销售量 记为x既然可以预测该型号飞机全球队销售量 结合波音公司的市场占有率 可以得到 5 制造成本 根据波音产品分析部门的估计 制造成本为 6 利润 假设利润等于销售收入去掉成本 并且公司的最优策略原则为利润R p 最大 利润函数为 由以上简化的分析及假设得到波音公司飞机最佳定价策略的数学模型如下 其中 2 5 3模型求解 我们采用图形的方法求解 具体用MATLAB作出目标函数曲线图 得到一个直观的印象 最优定价策略下价格p大致在6到7之间 再用图形放大方法 进一步估计出 如图2 3 MATLAB程序如下 作函数曲线图的基本程序 注意 1 根据图形的具体情况 不断修改上面程序中的最长一条语句 就可以不断地放大图形 将最优解的范围限制得越来越小 直至找出满意的近似解 2 以上市场占有率h 0 5 对于市场占有率h的其他取值 可以类似地进行 2 5 4进一步思考的几个问题 1 h取其他值时的最优价格 并进行比较 2 这是一个最值问题 由高等数学的知识 可以利用导数求驻点然后求最值 给出用此方法得到最佳价格p 精确到小数点后四位 的求解过程及MAT 利润曲线R P 1780 8336 1780 8336 1780 8336 1780 8335 1780 8334 1780 8334 1780 8334 1780 8333 1780 8332 6 2856 28526 28546 28566 28586 2866 28626 28646 28666 28686 287 价格p 图2 3图形放大法求解波音公司飞机最佳定价策略 LAB程序 3 如果模型假设中 在预测该型号飞机的全球销售量时 使用的不是二次函数 而是其他符合市场规律的曲线 具体考虑几种不同曲线 并进行计算和比较 4 以上问题的6条假设中 哪些较为合理 哪些不太合理 应该如何修改 5 在将此模型推向实际应用时 哪些因素是关键的 哪些因素处理和参数的获取是很困难的 2 6实验 2 6 1实验一 油价与船速的优化问题 油价的上涨 将影响大型海船确定合理的航行速度 以优化航行收入 直观地 油耗的多少直接影响船厂速成的快慢 因而直接影响航行时间的长短 进而影响支会船员人工费用数量 过去有一些经验表明 1 油耗正比于船速的立方 2 在最省油船速成的基础上改变速20 的速度 则引起50 的油耗的变化 作为一个例子 某中型海船厂 每天油耗40吨 减少20 的航速 省油50 即20吨 每吨油价250美元 由此每天减少耗油费用5000美元 而航行时间的增加将增加对船员支付的费用 如何最优化 算例 航程L 1536海哩 标准最省油航速20节 油耗每天50吨 航行时间8天 最低航速10节 本次航行总收入为84600美元 油价250美元 吨 是固定开支1000美元 试确定最佳航速 v 10 0 1 50 船速向量l 24 v 日路程向量t 1536 l 总行程天数zf 1000 t 航程的固定费用ryf 250 160 v 3 t 航程的油料费用plot v zf ryf 对不同航速 描绘航速 费用曲线 10节费用最少plot v 84600 zf ryf 对不同航速 描绘航速 利润曲线 10节费用最多 2 6 2实验二 确定培训项目的总费用 单位时间内接受培训的总人数是一定的 每批培训有一定的固定成本 接受培训的人员在培训结束后 未找到工作 待业 期间 培训单位需要支付培训人员一定的工资 试建立切合实际的数学模型 确定出每批培训的人数及周期 参考文献 1 萧树铁主编 姜启源 何表 高立编著 数学实验 高等教育出版社 1999 2 姜启源 数学模型 高等教育出版社 第二版 1993 3 任善强 雷呜 数学模型 重庆大学出版社 1996 4 DQuinney AnIntroductiontotheNumericalSoLutionofDifferentiaLEquations JohnWiley SonsInc 1985 5 G Fulford Modellingwithdifferentialanddifferenceequations Cambridge CambridgeUniversi tyPress 1997 6 peitgen Heinz Otto Chaosandfractals newfrontiersofscience NewYork Springer Verlag 1992 本章