您现在的位置：数值计算原理
数值计算原理
在线支付满98免快递费
作　者：，，编著
ISBN：6出版时间：重印）页数：461
包装：简裝本开本：20cm字数：
《数值计算原理》简介：　　本书的内容是现代科学计算中常用的数值计算方法及其原理，包括数值逼近，插值与拟合，数值积分，线性与非线性方程组数值解法，矩阵特征值与特征向量计算，常微分方程初值问题、刚性问题与边值问题数值方法，以及并行算法概述等。本书是为学过少量《计算方法》的理工科研究生学习《数值分析》而编写的教材。内容较新，起点较高，叙述严谨，系统性强，偏重数值计算一般原理。每章附有习题及数值试验题，附录介绍了Matlab软件以便于读者使用。本书可作为理工科研究生《数值分析》课程的教材或参考书，也可供从事科学与工程计算的科技人员学习参考。
暂缺《数值计算原理》作者简介
《数值计算原理》目录：第1章数值计算原理与计算精确度1数值计算的一般原理1-1数学问题与数值计算1-2数值问题与算法1-3数值计算的共同思想和方法2数值计算中的精确度分析2-1误差来源与误差估计问题2-2算法的数值稳定性2-3病态问题与条件数3并行算法及其基本概念3-1并行算法及其分类3-2并行算法基本概念3-3并行算法设计与二分技术评注习题数值实验题第2章数值逼近与数值积分1函数逼近的基本概念1-1数值逼近与函数空间1-2范数与赋范空间1-3函数逼近与插值1-4内积与正交多项式2多项式逼近2-1最佳平方逼近与勒让德展开2-2曲线拟合的最小二乘法2-3最佳一致逼近与切比雪夫展开3多项式插值与样条插值,3-1多项式插值及其病态性质3-2三次样条插值3-3B-样条函数4有理逼近4-1有理逼近与连分式4-2有理插值4-3帕德逼近5高斯型求积公式5-1代数精确度与高斯型求积公式5-2高斯-勒让德求积公式5-3高斯-切比雪夫求积公式5-4固定部分节点的高斯型求积公式6积分方程数值解7奇异积分与振荡函数积分计算7-1反常积分的计算7-2无穷区间积分7-3振荡函数积分8计算多重积分的蒙特卡罗方法8-1蒙特卡罗方法及其收敛性8-2误差估计8-3方差缩减法8-4分层抽样法8-5等分布序列评注习题数值实验题第3章线性代数方程组的数值解法1引言.线性代数的一些基础知识1-1引言1-2向量空间和内积1-3矩阵空间和矩阵的一些性质1-4向量的范数1-5矩阵的范数1-6初等矩阵2Gauss消去法和矩阵的三角分解2-1Gauss顺序消去法2-2矩阵的三角分解.直接三角分解解法2-3选主元的消去法和三角分解2-4对称正定方程组3矩阵的条件数与病态方程组3-1矩阵的条件数与扰动方程组的误差界3-2病态方程组的解法4大型稀疏方程组的直接方法4-1稀疏矩阵及其存储4-2稀疏方程组的直接方法介绍4-3带状方程组的三角分解方法4-4三对角和块三对角方程组的追赶法和循环约化方法5迭代法的一般概念5-1向量序列和矩阵序列的极限5-2迭代法的构造5-3迭代法的收敛性和收敛速度5-4J法和GS法的收敛性6超松弛迭代法6-1超松弛迭代法和对称超松弛迭代祛6-2超松弛迭代法的收敛性6-3块迭代方法6-4模型问题的红黑排序7极小化方法7-1与方程组等价的变分问题?-2最速下降法7-3共轭梯度法7-4预处理共轭梯度方法7-5多项式预处理评注习题数值实验题第4章非线性方程组数值解法1引言1-1非线性方程组求解问题1-2几类典型非线性问题2向量值函数的导数及其性质2-1连续与可导2-2导数性质与中值定理3压缩映射与不动点迭代法3-1压缩映射与不动点定理3-2不动点迭代法及其收敛性4牛顿法与牛顿型迭代法4-1牛顿法及其收敛性4-2牛顿法的变形与离散牛顿法4-3牛顿松弛型迭代法5拟牛顿法与Broyden方法5-1拟牛顿法基本思想5-2秩1拟牛顿法与Broyden方法6延拓法6-1延拓法基本思想6-2数值延拓法6-3参数微分法7并行多分裂方法7-1线性多分裂方法7-2非线性多分裂方法8非线性最小二乘问题数值方法评注习题数值实验题第5章矩阵特征值问题的计算方法1特征值问题的性质和估计1-1特征值问题的性质1-2特征值的估计1-3特征值的扰动2正交变换和矩阵分解2-1Householder变换2-2Givens变换2-3矩阵的QR分解2-4矩阵的Schur分解2-5正交相似变换化矩阵为Hessenberg形3幂迭代法和逆幂迭代法3-1幂迭代法3-2加速技术(Aitken方法)3-3收缩方法3-4逆幂迭代法4QR算法4-1QR迭代的基本算法和性质4-2Hessenberg矩阵的QR方法4-3带有原点位移的QR方法4-4双重步QR方法5对称矩阵特征值问题的计算5-1对称QR方法5-2Rayleigh商加速和Rayleigh商迭代5-3Lanczos方法评注习题数值实验题第6章常微分方程数值方法1初值问题数值方法1-1数值方法概述1-2局部截断误差与相容性1-3收敛性与稳定性1-4绝对稳定性与绝对稳定域2刚性微分方程及其数值方法的稳定性概念2-1刚性方程组2-2稳定性概念的扩充3解刚性方程的线性多步法3-1吉尔方法及其改进3-2含二阶导数的线性多步法3-3隐性问题与迭代法4隐式龙格-库塔法4-1龙格-库塔法的一般结构4-2基于数值求积公式的隐式RK方法4-3稳定性函数与隐式RK方法的A-稳定性4-4对角隐式RK方法5非线性方法6边值问题数值方法6-1打靶法6-2差分法评注习题数值实验题附录A数学软件Matlab入门附录BMatlab的工具箱附录C其他数学软件工具概览参考文献索引
商品问答（0条）
暂时没有问答
同类图书热卖榜
&12.2042折
&26.1045折
&12.6045折
&11.2043折
&18.2065折
&20.3026折
&16.0050折
&20.2023折
淘书推荐特价好书
&32.2046折
&16.8043折
&36.0045折
&10.3015折
合作单位：您所在位置： &
&nbsp&&nbsp&nbsp&&nbsp
解刚性问题一类指数拟合方法.pdf 53页
本文档一共被下载：
次 ,您可全文免费在线阅读后下载本文档。
下载提示
1.本站不保证该用户上传的文档完整性，不预览、不比对内容而直接下载产生的反悔问题本站不予受理。
2.该文档所得收入(下载+内容+预览三)归上传者、原创者。
3.登录后可充值，立即自动返金币，充值渠道很便利
解刚性问题一类指数拟合方法
你可能关注的文档：
··········
··········
刚性常常是实际科学研究中严重干扰数值解稳定和精度的一个重要因素，而刚
性微分方程数值积分方法的研究也已经成为了数值积分方法中一个最为活跃的研究
根据 Lambert
提出的在积分的局部区间上用一个有理函数来近似地表示刚性问
题微分方程的解的基本思想，我们在解的每个积分局部区间上构造了一个指数拟合
的函数,使其近似逼近微分方程的解曲线。
围绕指数拟合这一主题，本文构造了一类通过变量代换改进的指数拟合的方法。
这类变量代换的实质相当于在刚性微分方程解的快变区间作用一个衰减函数，使之
用常规的低阶算法就能很方便准确的求解。通过与普通积分算法的比较分析，作用
变量代换之后的方法无论在收敛速度还是在保持解的稳定性方面均高于未作用变量
代换之前的方法。这类变量代换构造的改进方法实现简单，计算方便，理论分析和
数值试验证明,在不降低原有积分方法相容性和收敛性的前提下,将该代换用于普通
的数值积分算法之后，能取得比普通的数值算法更好的收敛性,稳定性及计算精度。
此外，我们对这类指数拟合方法进行了误差分析和相应的改进，并在此基础上
构造了矩阵指数拟合的显式欧拉方法。实验结果表明，矩阵指数的拟合方法在收敛
速度和计算精度上均要远远高于直接用变量代换改进的方法。基于Lawson 将刚性方
程换成非刚性方程的思想,我们也类似地得到了广义Runge-Kutta 方法，并以此为基
础，构建了基于指数拟合的 Runge-Kutta 方法，理论分析证明，该方法具有良好的
关键词：刚性问题 变量代换 指数拟合 矩阵指数
factor,which
algebraic stability and computational accuracy in the science study 。The algebraic integral
methods of the stiff differential equations have been one of the most active study aspect in
the algebraic integral methods 。
On the basis of the theory of Lambert,approximately to approach the real value of the
stiff differentical equation by using an rational function，we construct an exponentially
fitted method to approach the value of the differential equation in the every local intergral
exponentially
exponentially
differential
substitution。The essential of this variable substitution is equivalent to effect an attenuation
function on the fast change area of the value of the differential equation,then
we can gain the result more conveniently
and more exactly by using a routine low-rank
integral method 。compared to the normal m
正在加载中，请稍后...常微分方程初值问题数值解法朱欲辉（浙江海洋学院 数理信息学院, 浙江 舟山 316004） [摘要]：在常微分方程的课程中讨论的都是对一些典型方程求解析解的方法. 然而在生产实 际和科学研究中所遇到的问题往往很复杂, 在很多情况下都不可能给出解的解析表达式. 本篇文章详细介绍了常微分方程初值问题的一些数值方法, 导出了若干种数值方法, 如 Euler法、 改进的Euler法、 Runge－Kutta法以及线性多步法中的Adams显隐式公式和预测校正 公式, 并且对其稳定性及收敛性作了理论分析. 最后给出了数值例子, 分别用不同的方法 计算出近似解, 从得出的结果对比各种方法的优缺点. [关键词]：常微分方程; 初值问题; 数值方法; 收敛性; 稳定性; 误差估计Numerical Method for Initial-Value ProblemsZhu Yuhui(School of Mathematics, Physics, and Information Science, Zhejiang Ocean University, Zhoushan, Zhejiang 316004)[Abstract]: In the course about ordinary differential equations, the methods for analyticsolutions of some typical equations are often discussed. However, in scientific research, the problems are very complex and the analytic solutions about these problems can’t be expressed explicitly. In this paper, some numerical methods for the initial-value problems are introduced. these methods include Euler method, improved Euler method, Runge-Kutta method and some linear multistep method (e.g. Adams formula and predicted-corrected formula). The stability and convergence about the methods are presented. Some numerical examples are give to demonstrate the effectiveness and accuracy of theoretical analysis.[Keywords]: Ordinary d Initial- NC S Error estimate1 前言自然界和工程技术中的很多现象, 例如自动控制系统的运行、 电力系统的运行、 飞行器 的运动、 化学反应的过程、 生态平衡的某些问题等, 都可以抽象成为一个常微分方程初值问 题. 其真解通常难以通过解析的方法来获得, 至今有许多类型的微分方程还不能给出解的 解析表达式, 一般只能用数值的方法进行计算. 有关这一问题的研究早在十八世纪就已经 开始了, 特别是计算机的普遍应用, 许多微分方程问题都获得了数值解, 从而能使人们认 识解的种种性质及其数值特征, 为工程技术等实际问题提供了定量的依据. 关于常微分方程初值问题的数值计算方法, 许多学者己经做了大量的工作. 1768年, Euler提出了关于常微分方程初值问题的方法, 1840年, Cauchy第一次对初值问题进行了仔 细的分析, 早期的常微分方程数值解的问题来源于天体力学. 在1846年, 当Adams还是一 个学生的时候, 和Le Verrier一起根据天王星轨道中出现的己知位置, 预测了它下一次出现 的位置. 1883年, Adams提出了Adams一Bashforth和Adams一Moulton方法. Rull(1895年)、 Heun(1900年)和Kutta(1901年)提出Runge.Kutta方法. 二十世纪五十年代, Dahlquist建立了常微分方程数值解法的稳定性理论, 线性多步法 是常微分方程初值问题的一种数值方法. 由于通常的数值方法, 其绝对稳定区域是有限的, 不适用于求解刚性常微分的初值问题. 刚性微分方程常常出现于航空、航天、热核反应、自 动控制、 电子网络及化学动力学等一系列与国防和现代化建设密切相关的高科技领域, 具有 无容置疑的重要性. 因此, 刚性微分方程的研究工作早在二十世纪五十年代就开始了, 1965年, 在爱丁堡举行的IFIP会议后, 更进一步地认识刚性方程的普遍性和重要性. 自从 六十年代初, 许多数值分析家致力于探讨刚性问题的数值方法及其理论, 注意到刚性问题 对传统数值积分方法所带来的挑战. 这一时期, 人们的研究主要集中在算法的线性稳定性 上, 就是基于试验方程 y，? ? y，（? ? C） 数值解的稳定性研究. 在此领域发表了大量的 论文, 取得了许多重要的理论成果. 例如, 1963年, Dahlquist给出A稳定性理论, 1967年, Widlund给出 A(? ) ―稳定性理论, 1969年, Gear将 A ―稳定性减弱, 给出刚性(Stiff)稳定性 理论, 并找到了当 k ? 6 的k步k阶的刚性稳定方法, 1969年Dill找到刚性稳定的7阶和8阶以 及1970年Jain找到刚性稳定的9阶到11阶, 但可用性没有检验. 这些稳定性理论和概念都是 在线性试验方程的框架下推导出的, 从严格的数学意义上来说, 这些理论只适用于常系数 线性自治系统. 但从实用的观点来说, 这些理论无疑是合理和必要的, 对刚性问题的算法设计具有重要的指导意义. 在八十至九十年代, 国内也有一些学者研究线性理论, 主要有 匡蛟勋、陈果良、项家祥、李寿佛、黄乘明、李庆扬和费景高等. 线性理论虽然对一般问题具有指导作用, 但其不能作为非线性刚性问题算法的稳定性 理论研究基础. 为了将线性理论推广到非线性问题中, 人们开始对非线性模型问题进行研 究.但是, 早期文献主要致力于数值方法基于经典Lipschitz条件下的经典收敛理论, 即认为 良好的稳定性加上经典相容性和经典相容阶就足以描述方法的整体误差性态. 直到1974 年, Prothero和Robinson首先注意到算法的经典误差估计由于受刚性问题巨大参数的影响而严重 失真, 产生阶降低现象, 这时人们认识到经典收敛理论对于非线性刚性问题以及线性模型 的不足. 于是, 1975年, Dahlquist和Butcher分别提出了单支方法和线性多步法的G一稳定概 念和B稳定概念. 这两个概念填补了非线性稳定性分析理论, 引起了计算数学家们的极大关 注, 在上述理论的基础上, 1975年至1979年, Burrage和Butcher提出了AN一稳定性与BN一 稳定性概念, 并相应地建立了基本的B一稳定及代数稳定理论. 年, Frank, schneid和ueberhuber建立Runge一kutta方法的B一收敛理论. B稳定与B一收敛理论统称B一理 论, 它是常微分方程数值解法研究领域的巨大成就之一, 是刚性问题算法理论的突破性进 展, 标志着刚性问题研究从线性向非线性情形深入发展. 国内也有众多学者致力于B一理 论的研究, 如李寿佛、曹学年等. 1989年, 李寿佛将Dahlquist的G一稳定概念推广到更一般 的(C,P,Q)代数稳定, 克服了G稳定的线性多步法不能超过二阶的限制. 对于一般线性方法, 李 寿 佛 建 立 了 一 般 线 性 方 法 的 (K,P,Q) 稳 定 性 理 论 及 (K,P,Q) 弱 代 数 稳 定 准 则 和 多 步 Runge-Kutta法的一系列代数准则. 此外, Dahquist, Butcher和Hairer分别深刻地揭示了单支 方法、一般线性方法和Runge.Kutta方法线性与非线性稳定性之间的内在关系. 为了求 解刚 性微分方程, 不少文献中构造含有稳定参数的线性多步方法, 利用适当选择稳定参数来扩 大方法的稳定区域. 所有改进的思想, 都是通过构造一些特殊的显式或隐式线性多步法,使其具有增大的稳定域, 或使 A(? ) 一稳定的 ? 角增大. 八十年代, 就成为国内外学者所 研究的一个课题, 学者主要有Rodabaugh和Thompson、Feinberg、李旺尧、李寿佛、包雪松、 徐洪义、刘发旺、匡蛟勋、项家祥、蒋立新、李庆扬、谢敬东和李林忠等.当前国内外研究 刚性问题的一个主要趋势就是在B一理论指导下寻找更为有效的新算法. 另一个发展趋势就 是力图突破单边lipschitz常数和内积范数的局限, 建立比B一理论更为普遍的定量分析收敛 理论. 近年来, 刚性延迟系统的算法研究成为刚性问题的另一个热点研究领域, 张诚坚将 Burrage等人创立的针对刚性常微分系统的B一理论拓展到非线性刚性延迟系统.常微分方程的数值算法发展到今天己有了线性多步法、 龙格一库塔法和在此基础上发展 起来的单支方法、分块方法、循环方法、外推法、混合方法、二阶导数法以及各种常用的估 校正算法. 其中经常用到的线性多步法公式有Euler公式、Heun公式、中点公式、Milne公 式、Adams公式、simpson公式、Hamming公式, Gear方法、Adams预估一校正法和Mile预估 一Hamming校正法公式等, 此外还包含许多迄今尚末探明的新公式. Burage曾将线性多步 法和Runge―Kutta法比作大海中的两座小岛, 在浩瀚的汪洋之中, 还有许多到现在没有发 现的新方法. 本篇文章详细介绍了常微分方程初值问题的一些数值方法, 导出了若干种数值方法, 如Euler法、改进的欧拉法、显式龙格-库塔法、隐式龙格-库塔法以及线性多步法中的Adams 显隐式公式和预测校正公式, 并且对其稳定性及收敛性作了理论分析. 最后给出了数值例 子, 进行了计算机程序算法的分析与实现, 以计算机的速度优势来弥补计算量大的不足. 从得出的结果对比各种方法的优缺点.2 常微分方程初值问题的数值解法2.1 数值方法的基本思想考虑一阶常微分方程初值问题dx ? f ( x, y ), x ? [a, b], dyy( x0 ) ? y0的数值解法, 其中 f 是 x 和 y 的己知函数, y0 是给定的初始值.(2.1)对于常微分方程初值问题(2.1)数值解法, 就是要算出精确解 y ( x) 在区间 [ a, b] 上的一 系列离散节点 a ? x0 ? x1 ? ? ? xn? 1 ? xn ? b 的函数值 y( x0 ) , y ( x1 ) , ?, y( xn ) 的近似 值 y0 , y1 , ? , yn .相邻两个节点的间距 h ? xi ?1 ? xi 称为步长, 这时节点也可以表示为xi ? x0 ? ih (i ? 1, 2,?, n) . 数值解法需要把连续性的问题加以离散化, 从而求出离散节点的数值解. 通常微分方程初值问题(2.1)的数值方法可以分为两类:(1) 单步法-----计算 y ( x) 在 x ? xn?1 处的值仅取决于 xn 处的应变量及其导数值. (2) 多步法-----计算 y ( x) 在 x ? xn?1 处的值需要应变量及其导数在 xn ?1 之前的多个网个节 点出的值.2.2 Euler 方法 2.2.1 Euler 公式若将函数 y ( x) 在点 xn 处的导数 y' ( xn ) 用两点式代替, 即y ' ( xn ) ?y ( xn ?1 ) ? y ( xn ) , h再用 yn 近似地代替 y( xn ) , 则初值问题(2.1)变为? yn ?1 ? yn ? hf ( xn , yn ) ? ? y0 ? y ( x0 ), n ? 0,1, 2,?(2.2)式就是著名的欧拉(Euler)公式.(2.2)2.2.2 梯形公式欧 拉 法 形 式 简 单 精 度 低 , 为 了 提 高 精 度 , 对 方程 y ? f ( x, y) 的 两端 在 区 间'[ xn , xn?1 ] 上积分得y( xn?1 ) ? y( xn ) ? ?改用梯形方法计算其积分项, 即xn?1 xnf [ x, y( x)]dx(2.3)?xn?1xnf [ x, y ( x)]dx ?xn ?1 ? xn [ f ( xn , y ( xn )) ? f ( xn ?1 , y ( xn ?1 ))] 2代入式(2.3), 并用 yn 近似代替式中 y( xn ) 即可得到梯形公式h yn ?1 ? yn ? [ f ( xn , yn ) ? f ( xn ?1 , yn ?1 )] 2度高一个数值方法.(2.4)由于数值积分的梯形公式比矩形公式精度高, 所以梯形公式(2.4)比欧拉公式(2.2)的精式(2.4)的右端含有未知的 yn ?1 , 它是关于 yn ?1 的函数方程, 这类方法称为隐式方法.2.2.3 改进的欧拉公式梯形公式实际计算时要进行多次迭代, 因而计算量较大. 在实用上, 对于梯形公式 (2.4)只迭代一次, 即先用欧拉公式算出 yn ?1 的预估值 yn?1 , 再用梯形公式(2.4)进行一次 迭代得到校正值 yn ?1 , 即? y n?1 ? yn ? hf ( xn , yn ) ? ? h ? yn?1 ? yn ? [ f ( xn , yn ) ? f ( xn ?1 , y n?1 )] ? 22.2.4 欧拉法的局部截断误差(2.5)衡量求解公式好坏的一个主要标准是求解公式的精度, 因此引入局部截断误差和阶数 概念. 定 义 2.1 在 yn 准 确 的 前 提 下 , 即 yn ? y( xn ) 时 , 用 数 值 方 法 计 算 yn ?1 的 误 差Rn?1 ? y( xn?1 ) ? yn?1 , 称为该数值方法计算 yn?1 时的局部截断误差.对于欧拉公式, 假定 yn ? y( xn ) , 则有yn?1 ? y( xn ) ? h[ f ( xn , y( xn ))] ? y( xn ) ? hy' ( xn )而将真解 y ( x) 在 xn 处按二阶泰勒展开式有h2 ''' yn?1 ? y ( xn ) ? hy ( xn ) ? y (? ), ? ? ( xn , xn ?1 ) 2!'因此有h 2 '' y ( xn ?1 ) ? yn ?1 ? y (? ) 2!定义 2.2 若数值方法的局部截断误差为 O(hp ?1) , 则称这种数值方法的阶数是 P. 步长(h&1)越小, P 越高, 则局部截断误差越小, 计算精度越高.2.3Runge－Kutta 方法2.3.1 Runge－Kutta 方法的基本思想欧拉公式可改写成? yn ?1 ? yn ? hK1 ? ? K1 ? f ( xn , yn )则 yn ?1 的表达式与 y ( xn ?1 ) 的泰勒展开式的前两项完全相同, 即局部截断误差为 O(h2 ) . 改进的欧拉公式又可改写成h ? ? yn ?1 ? yn ? 2 ( K1 ? K 2 ) ? . K1 ? f ( xn , yn ) ? ? K ? f ( x , y ? hk ) n ?1 n 1 ? 2 ?上述两组公式在形式上有一个共同点: 都是用 f ( x, y) 在某些点上值的线形组合得出y( xn?1 ) 的近似值 yn?1 , 而且增加计算 f ( x, y) 的次数, 可提高截断误差的阶. 如欧拉公式,每步计算一次 f ( x, y) 的值, 为一阶方法. 改进的欧拉公式需计算两次 f ( x, y) 的值, 它是 二阶方法. 它的局部截断误差为 O(h3 ) . 于是可考虑用函数 f ( x, y) 在若干点上的函数值的线形组合来构造近似公式, 构造时 要求近似公式在 ( xn , yn ) 处的泰勒展开式与解 y ( x) 在 xn 处的泰勒展开式的前面几项重合, 从而使近似公式达到所需要的阶数. 既避免求偏导, 又提高了计算方法精度的阶数. 或者 说, 在 [ xn , xn?1 ] 这一步内多预报几个点的斜率值, 然后将其加权平均作为平均斜率. 则可 构造出更高精度的计算格式, 这就是 Runge－Kutta 方法的基本思想.2.3.2 Runge－Kutta 方法的构造一般地, Runge－Kutta 方法设近似公式为 （下面的公式修改了）p ? yn ?1 ? yn ? h? ci K i ? i ?1 ? ? K1 ? f ( xn , yn ) ? ? i ?1 ? Ki ? f ( xn ? ai h, yn ? h? bij K j ) ? j ?1 ??i ? 2,3,?, p?其中 ai , bij , ci 都是参数, 确定它们的原则是使近似公式在 ym 处的泰勒展开式与 y ( x) 在xn 处的泰勒展开式的前面的项尽可能多的重合, 这样就使近似公式有尽可能高的精度. 以此我们可以通过一个复杂的计算过程得出常用的的三阶和四阶 Runge－Kutta 公式h ? ? yn ?1 ? yn ? 6 ( K1 ? 4 K 2 ? K 3 ) ? K1 ? f ( xn , yn ) ? ? ? K ? f (x ? h , y ? h K ) 2 n n 1 ? 2 2 ? K ? f ( x ? h, y ? hK ? 2hK ) n n 1 2 ? 3和h ? ? yn ?1 ? yn ? 6 ( K1 ? 2 K 2 ? 2 K 3 ? K 4 ) ? K1 ? f ( xn , yn ) ? ? h h ? K 2 ? f ( xn ? , yn ? K1 ) ? 2 2 ? ? h h K 3 ? f ( xn ? , yn ? K 2 ) ? 2 2 ? ? ? K 4 ? f ( xn ? h, yn ? hK 3 )式(2.6)称为经典 Runge－Kutta 方法.(2.6)2.4 线性多步法在逐步推进的求解过程中, 计算 yn ?1 之前事实上已经求出了一系列的近似值 y0 , y1 ,? , yn . 如果充分利用前面多步的信息来预测 yn?1 , 则可期望获得较高的精度, 这就是构造多步法的基本思想.线性 k 步方法的一般公式为k ?1 j ?0 k ?1yn?1 ? ? a j yn? j ? h ? b j f ( xn? j , yn? j )j ??1(2.7)其中 a j , bj 均为与 n 无关的常数, ak ?1 ? bk ?1 ? 0 . 当 b?1 ? 0 时为显格式; 当 b?1 ? 0 时 为隐格式. 特别当 k ? 1, a0 ? b0 ? 1, b?1 ? 0 时为 Euler 公式； k ? 1, a0 ? 1, b0 ? b?1 ? 当 为梯形公式.1 时 2定义 2.3 称Rn?1 ? y( xn?1 ) ? [ yn?1 ? ? a j yn? j ? h ? b j f ( xn? j , yn? j )]j ?0 j ??1k ?1k ?1为 k 步公式(2.7)在 xn ?1 处的局部截断误差. 当Rn?1 ? O(h p?1 )时称式(2.7)是 p 阶的.应用方程 y' ( x) ? f ( x, y( x)) 可知局部截断误差也可写成为Rn?1 ? y( xn?1 ) ? [ yn?1 ? ? a j yn? j ? h ? b j y ' ( xn? j )]j ?0 j ??1k ?1k ?1定义 2.4如果线性 k 步方法(2.7)至少是 1 阶的, 则称是相容的; 如果线性 k步法(2.7)是 p 阶的, 则称是 p 阶相容的. 2.4.1 Adams 外插法将微分方程dx ? f ( x, y ), x ? [a, b], dy的两端从 xn 到 xn ?1 进行积分, 得到y( xn?1 ) ? y( xi ) ? ?我们用插值多项式代替右端的被积函数.xn?1xnf ( x, y( x))dx(2.8)Adams 外插法选取 k 个点 xn , xn ?1,? ,?, xn?k ?1 作为插值基点构造 f ( x, y) 的 k-1 阶多项式 Adams 外插法的计算公式为yn?1 ? yn ? h? a j ? j f mj ?0k ?1其中 a j 满足如下代数递推式:1 1 1 a j ? a j ?1 ? a j ?2 ? ? ? a0 ? 1 , 2 3 j ?1j ? 0,1,?根据此递推公式, 可逐个的计算 a j? j ? 0,1,?? ,表 2.1表 2.1 给出了 a j 的部分数值:j0 11 1/22 5/123 3/84 251/7205 95/288?aj?2.4.2 Adams 内插法根据插值理论知道, 插值节点的选择直接影响着插值公式的精度, 同样次数的内插公 式的精度要比外插公式的高.仍假定已按某种方法求得问题(2.1)的解 y ( x) 在 xi (i ? 0,1,?n) 处的数值 yi , 并选取插 值 节 点 xn? k ? ,? , x n ? ,x n? p p 是 正 整 数 , 用 Lagrange 型 插 值 多 项 式 Ln,k ?1 ( x) 构 造 , , 1( p)y ' ? f ( x, y)可以导出解初值问题(2.1)的Adams内插公式:y ( xn ?1 ) ? y ( xn ) ?当 p ? 0 时上式就退化为内插公式. 内插公式(2.9)除了包含 f 在 xn?k ?1 ,?, xn 处的已知值外, 还包含了在点 xn ,?, xn? p , 处的未知值.因此内插公式(2.9)只给出了未知量 yn ,?, yn? p 的关系式, 实际计算时仍需要xn?1?L( p) n , k ?1( x)dx(2.9)xn解联立方程组. p ? 1 的内插公式是最适用的, Ln,k ?1 ( x) 采用Newton向后插值公式得到(1)Adams内插公式yn?1 ? yn ? h? a* j ? j f mj ?0k其中系数 a j 定义为0 ? ?? ? a* j ? (?1) j ? ? ?d? ?1 ? j ?*j ? 0,1,?其中 a j 满足如下代数递推式:*1 1 1 * ?1, j ? 0 a* j ? a* j ?1 ? a* j ?2 ? ? ? a 0 ?? 2 3 j ?1 ?0, j ? 0根据此递推公式, 可逐个的计算 a* j? j ? 0,1,?? ,表 2.2表 2.2 给出了 a j 的部分数值:*j0 11 -1/22 -1/123 -1/244 -19/7205 -3/160? ?a* j2.5 算法的收敛性和稳定性 2.5.1 相容性初值问题(2.1):dx ? f ( x, y ), x ? [a, b], dyy( x0 ) ? y0的显式单步法的一般形式为yn?1 ? yn ? h? ( xn , yn , h) ,其中 h ?n ? 0,1,?, N ? 1(2.10)b?a , xn ? a ? nh . 这样我们用差分方程初值问题(2.10)的解 yn 作为问题(2.1)的 N解 y( xn ) 在 x ? xn 处的近似值, 即 y( xn ) ? yn .因此, 只有在问题(2.1)的解使得y ( x ? h) ? y ( x ) ? ? ( x, y ( x), h) h逼近y '? f ( x, y( x)) ? 0时, 才有可能使(2.10)的解逼近于问题（2.1）的解. 从而, 我们期望对任一固定的 x ?[a, b], 都有lim[h ?0y ( x ? h) ? y ( x ) ? ? ( x, y ( x), h)] ? 0 h假设增量函数 ? ( x, y, h) 关于 h 连续, 则有? ( x, y, h) ? f ( x, y)定义 2.5 若关系式? ( x, y, h) ? f ( x, y)成立, 则称单步法(2.10)与微分方程初值问题(2.1)相容. 容易验证, Euler 法与改进 Euler 法均满足相容性条件. 事实上, 对 Euler 法, 增量函数 为? ( x, y, h) ? f ( x, y)自然满足相容性条件, 改进 Euler 法的增量函数为? ( x, y, h) ? [ f ( x, y ) ? f ( x ? h, y ? hf ( x, y))]因为 f ( x, y) , 从而有1 2? ( x, y, h) ? [ f ( x, y ) ? f ( x, y )] ? f ( x, y )所以 Euler 法与改进 Euler 法均与初值问题（2.1）相容. 一般的, 如果显示单步法有 p 阶精 度 ( p ? 0) , 则其局部误差为1 2y( x ? h) ? [ y( x) ? h? ( x, y( x), h)] ? O(h p?1 )上式两端同除以 h, 得 ? n?1 ? 0 令 h ? 0 , 如果 ? ( x, y, h) , 则有y' ( x) ? ? ( x, y, h) ? 0即? ( x, y, h) ? f ( x, y)所以 p ? 0 的显示单步法均与初值问题(2.1)相容.由此各阶 RK 方法与初值问题(2.1)相容.2.5.2 收敛性定义 2.6 一种数值方法称为是收敛的, 如果对于任意初值 y0 及任意 x ? (a, b], 都有lim yn ? y ( x)h ?0( x ? a ? nh)其中 y ( x) 为初值问题(2.1)的准确解. 按定义, 数值方法的收敛性需根据该方法的整体截断误差 ? n ?1 来判定. 已知 Euler 方法 的整体截断误差有估计式? n ?1 ?当h ? 0,hM L (b ?a ) e ? 1 ? O ( h) 2L? n?1 ? 0 , 故 Euler 方法收敛.定理 2.1 设显式单步法具有 p 阶精度, 其增量函数 ? ( x, y, h) 关于 y 满足 Lipschitz 条件, 问题(2.1)是精确的, 既 y( x0 ) ? y0 , 则显式单步法的整体截断误差为? n?1 ? y( xn?1 ) ? yn?1 ? O(h p ) .定理 2.2 设增量函数 ? ( x, y, h) 在区域 S 上连续, 且关于 y 满足 Lipschitz 条件时, 则显式单步法收敛的充分必要条件是相容性条件成立, 即? ( x, y, h) ? f ( x, y) .2.5.3 稳定性? 定义 2.7 如果存在正常数 h0 及 C, 使得对任意的初始出发值 y0 , y0 , 单步法(2.10)的 ? 相应精确解 yn , yn , 对所有的 0 ? h ? h0 , 恒有yn ? ? n ? C y0 ? ? 0 , y y则称单步法是稳定的.nh ? b ? a定理 2.3 若 ? ( x, y, h) 对于 x ? [a, b] , h ? (0, h0 ] 以及一切实数 y, 关于 y 满足 Lipschitz 条件, 则单步法(2.7)是稳定的. 定义 2.8 对给定的微分方程和给定的步长 h, 如果由单步法计算 yn 时有大小为 ? 的误 差, 即计算得出 ? n ? yn ? ? , 而引起其后之值 ym (m ? n) 的变化小于 ? , 则说该单步法 y是绝对稳定的.一般只限于典型微分方程y' ? ? y考虑数值方法的绝对稳定性, 其中 ? 为复常数(我们仅限于 ? 为实数的情形). 若对于所有 的 ? h ? (? , ? ) , 单步法都绝对稳定, 则称 (? , ? ) 为绝对稳定区间.根据以上定义我们可以 得出 Euler 方法的绝对稳定区间为(-2,0), 梯形公式的绝对稳定区间为 (??, 0) , Runge－ Kutta 的绝对稳定区间为(-2.78,0).3 实例分析例 3.1 分别用 Euler 法、改进 Euler 法、Runge－Kutta 解初值问题? y ' ? ?2 xy 2 (0 ? x ? 1.2) ? y (0) ? 1 ?的数值解, 取 h=0.1, N=10, 并与真实解作出比较.由于该简单方程可以用数学方法求得其精确描述式 y ?1 , 所以可以据此检验近 1 ? x2似数值解同真实解的误差情况. 对于其他一些结构复杂的常微分方程的数值解实现方法也 是一样的. （1）使用欧拉算法进行一般求解 经过欧拉算法程序计算得出 yi 在各点 xi 处的近似数值解及各自同精确解的误差, 如下 表3.1 表3.1x[i ] 各分点0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8y[i ] (数值解)1........622018y ( x[i ]) (精确解)0........609756y( x[i]) ? y[i] (误差)0........0122620.9 1.0 1.1 1.20....4077830....4098360....002053从实验结果看误差已经不算太小, 况且这还仅仅是一个用于实验的比较简单的常微分 方程, 对于实际工程中应用的结构复杂的方程其求解结果的误差要远比此大得多, 由于还 存在着局部截断误差和整体截断误差, 因此可采取一定措施来抑制减少误差, 尽量使结果 精确. （2）使用改进欧拉算法进行一般求解 经过改进欧拉算法程序计算得出 yi 在各点 xi 处的近似数值解及各自同精确解的误差, 如下表3.2 表3.2x[i ] 各分点0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2y[i ] (数值解)0............410859y ( x[i ]) (精确解)0............409836y( x[i]) ? y[i] (误差)0............001023可以看出, 这种经过改进的欧拉算法所存在的误差已大为减少, 误差的减少主要是 由于先利用了欧拉公式对 yi ?1 的值进行了预估, 然后又利用梯形公式对预估值作了校正, 从而在预估―校正的过程中减少了误差. 此算法有一定的优点, 在一些实际而且简单的工 程计算中可以直接单独应用.（3）使用经典Runge－Kutta法进行一般求解 经过经典Runge－Kutta法程序计算得出 yi 在各点 xi 处的近似数值解及各自同精确解的 误差, 如下表3.3 表3.3x[i ] 各分点0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2y[i ] (数值解)0............409837y ( x[i ]) (精确解)0............409836y( x[i]) ? y[i] (误差)0............000001从表3.3记录的程序运行结果来看, 所计算出来的常微分方程的数值解是非常精确的, 用浮点数进行运算时由近似所引入的误差几乎不会对计算结果产生影响, 这样高的精度是 很令人满意的. 无论是从计算速度上还是从计算精度要求上, 都能很好地满足工程计算的 需要.例 3.2 分别用 Euler 法、改进 Euler 法、Runge－Kutta 解初值问题? y ' ? ?2 xy 2 (0 ? x ? 1.2) ? y (0) ? 1 ?的数值解, 取等分数分别为 N=12, 24, 36, … ,120, 分别计算出与真实解的最大误差.（1） 使用欧拉算法进行一般求解经过计算我们可以得出初值问题在取不通步长时, 数值解与真实解的最大误差如表 3.4 表 3.4 N(等分数) 12 24 36 48 60 最大误差 0.....004878 N(等分数) 72 84 96 108 120 最大误差 0.....002416从表 3.4 记录的程序运行结果来看当等分数 N 变大时, 它的误差正在减小, 根据定义 2.6 我 们可以证得该方法是收敛的. （2）使用改进欧拉算法进行一般求解 经过计算我们可以得出初值问题在取不通步长时, 数值解与真实解的最大误差如表 3.5 表 3.5 N(等分数) 12 24 36 48 60 最大误差 0.....000041 N(等分数) 72 84 96 108 120 最大误差 0.....000010从表 3.5 记录的程序运行结果来看当等分数 N 变大时, 它的误差正在减小, 根据定义 2.6 我 们可以证得该方法是收敛的. （3）使用经典Runge－Kutta法进行一般求解 经过计算我们可以得出初值问题在取不通步长时, 数值解与真实解的最大误差如表 3.6 表 3.6 N(等分数) 最大误差 N(等分数) 最大误差12 24 36 48 600.....00000072 84 96 108 1200.....000000从表 3.6 记录的程序运行结果来看当等分数 N 变大时, 它的误差正在减小, 根据定义 2.6 我 们可以证得该方法是收敛的. （4）收敛速度对比 为对比以上三种方法的收敛速度, 我们计算出了它们的误差精度达到 10?5 的最小等分 数N如下表3.7 表3.7 方法 Euler法 改进Euler法 Runge－Kutta法 N(最小等分数)
误差精度10?510?5 10?5从表3.7记录的程序运行结果来看Runge－Kutta法的收敛速度最快, 改进Euler法其次, 而 Euler法最差. 由此看来Runge－Kutta法是他们当中最理想的数值解法.小结针对工程计算中遇到的常微分方程初值问题的求解, 根据实际情况引如了计算机作为 辅助计算工具, 并根据高等数学的有关知识将欧拉公式、 改进的欧拉公式、 经典龙格－库塔 方法等融入到计算机程序算法之中, 充分利用了计算机的速度优势, 大大减轻了工程技术 人员的劳动强度, 同时也使计算结果更加可靠、 准确。 但在实际应用时, 针对不同的工程函 数选择合适的求解方法需要有较高的要求, 既要考虑到方法的简易, 又要减少计算量, 同 时又不能让截断误差超出指定范围. 一般来说经典龙格－库塔算法精确度高又利于计算机 编程实现, 收敛速度快, 稳定性也很好, 可以考虑作为首选实现算法.参考文献[1]李寿佛. 多步 Runge-Kuta 方法的代数稳定性. 系统仿真学报, ): 51 一 56. [2]李寿佛, 曹学年. 一般多值方法的BH一代数稳定性. 高等学校计算数学学报, ): 67一77. [3]包雪松, 徐洪义, 王长富,翟灿芳. 关于具有增大绝对稳定域的线性多步法. 高等学校 计算数学学报, ): 311一318. [4]匡蛟勋, 项家样. 一类修正的BDF方法. 计算数学, ): 41一418. [5]蒋立新. Gear方法的改进. 湘潭大学自然科学学报,): l3一18. [6]李庆扬, 谢敬东. 刚性常微分方程的一种线性多步法. 清华大学学报(自然科学 版),): 1一11. [7]李林忠, 何轶良. 三阶隐式线性三步方法的A一稳定和A(a)一稳定性. 高等学校计算数 学报,): 158一168. [8]李林忠, 何轶良, 陈小飞. 显式线性三步二阶方法的稳定性讨论. 高等学校计算数学学 报, ): 42一53. [9]李大侃. 常微分方程数值解法. 杭州：浙江大学出版社，1994. [10] AtkinsonKE著. 数值分析引论. 匡蛟勋译. 上海: 上海科学技术出版社, 1986. [11] 李庆扬. 常微分方程数值解法(刚性问题与边值问题. 北京:高等教育出版社, 1992. [12] 关治, 陆金甫. 数值分析基础. 北京: 高等教育出版社, 1998. [13] 林成森. 数值计算方法. 北京: 科学出版社, 2005. [14] 刘师少. 计算方法. 北京: 科学出版社, 2005.致谢本人在撰写论文的过程中, 得到了许多老师和同学的热心帮助. 这次论文的成功撰写, 凝结了我四年大学生涯的心血. 感谢抚养我长大的父母, 感谢所有传授我知识的老师, 感 谢教会我独立自主生活的班主任, 感谢所有关心过我、帮助过我的恩师们. 这里, 我要特 别感谢黄红英老师对我的悉心指导和严格要求, 使我完成了论文的研究工作, 提高了论文 的写作水平和研究问题的能力. 此外, 黄老师详尽的审阅了论文初稿, 给我提出了宝贵的 修改意见, 对英文翻译也进行了逐字逐句的修改与校正. 黄老师严谨的治学态度和诲人不 倦的师者作风使我受益终身, 再次向黄老师表示衷心的感谢.