微分方程数值解【稻谷书苑】

《微分方程数值解【稻谷书苑】》由会员分享，可在线阅读，更多相关《微分方程数值解【稻谷书苑】（74页珍藏版）》请在金锄头文库上搜索。

1、微分方程数值解微分方程数值解 邓斌邓斌 合肥工业大学数学学院合肥工业大学数学学院 1详细课资 “计算数学”就是研究在计算机上解决数学问题的理论和数值方法。u今天的数值计算方法，无论从形式到内容，还是从工具到效果，已远非半世纪前Von Neumann、Lax等先驱们所处的环境和条件了，计算机技术和应用软件的发展，让计算数学展开了双翼。许多迅速发展的其他学科和社会进步给计算数学的发展开拓出 更为广阔的新天地。2详细课资u随着计算机软件硬件的不断更新和计算方法的迅速发展，科学计算科学计算与实验实验以及理论研究理论研究成为现代科学研究的三大主要手段。u科学计算还能解决实验及理论无法解决的问题，并由此发

2、现一些新的物理现象，加深人们对物理机理的理解和认识，促进科学的发展。 u作为三种科学研究手段之一的科学计算是一门工具性、方法性、边缘性的新学科，发展迅速，它的物质基础是计算机(包括其软硬件系统)，其理论基础主要是计算数学。3详细课资一. 引言24详细课资5详细课资 计算数学发展的历史回顾26详细课资 1.从计算物理谈起u计算数学的发展与科学工程计算是紧密相联的，计算数学的发展历史也就是与其他学科结合，利用计算机不断形成新的理论及数值方法并不断形成新的学科的历史，例如：“计算物理”。7详细课资年月美国总统发布命令，可以揭开曼哈顿计划的内幕，部分内容可以解密。故以“计算物理方法”丛书的名义陆续编辑

3、出版8详细课资 我们常说:u计算物理的物质基础是计算机；u计算物理的关键技术是“计算方法”和“程序设计”;u 计算物理发展的原始动力是美国核武器研制的刺激。9详细课资美国从1942年8月13日开始曼哈顿计划，到1945年制造出三颗原子弹：代号为：“三一”，用于试验（7月16日），“瘦子”投于广岛(8月6日)，“胖子”投于长崎(8月9日)。历时三年，涉及到理论物理、爆轰物理、中子物理、金属物理、弹体弹道等大量的数值计算。10详细课资1949年8月苏联第一次原子弹爆炸后，杜鲁门总统在1950年1月31日下令继续研究各种类型的原子弹武器，成立以氢弹之父特勒(E. Teller)为首的氢弹研制小组。直

4、到1952年10月31日爆炸了代号为“麦克”的核试验。11详细课资在研制原子弹和氢弹过程中，许多物理规律必须通过计算机上的计算摸清楚。计算物理、理论物理与实验物理相辅相成相互促进共同发展，形成现代物理学的三大分支。12详细课资由于核武器研制需要，1950年全球只有15台，到了1962年9月仅美国就有16187台计算机。13详细课资60年代中期开始推出小型计算机，70年代末推出个人计算机，80年代中期又推出高性能的超级微机。而计算物理发展所涉及的大规模科学计算和模拟所需要的大型计算机却未得到发展。14详细课资1981年以哈佛大学普雷斯(W. H. Press)为首的11位著名科学家联名上书，向美

5、国国家科学基金会(NSF)呈送“发展计算物理的建议书”，大声疾呼计算物理发展正处于一个危机阶段，是NSF采取实质性行动的时候了。15详细课资 .科学计算16详细课资1983年一个由美国著名数学家拉克斯(P. Lax)为首的不同学科的专家委员会向美国政府提出的报告之中，强调“科学计算是关系到国家安全、经济发展和科技进步的关键性环节，是事关国家命脉的大事。” 17详细课资1984年美国政府大幅度地增加对科学计算经费的支持, 新建成五个国家级超级计算中心（分别在普林斯顿大学、圣地亚哥、伊里诺大学、康奈尔大学、匹兹堡），配备当时最高性能的计算机，建立NSF-net新网络。18详细课资80年代中期我国将

6、“大规模科学与工程计算”列入国家资助重大项目。19详细课资1987年起美国NSF把“科学与工程计算”、“生物工程”“全局性科学”作为三大优先资助的领域。20详细课资1990年美国国家研究委员会发表振兴美国数学：90年代的计划的报告，建议对由计算引发的数学给予特殊的鼓励和资助。21详细课资报告指出由于大存储的高速计算机的使用已导致了科学和技术方面的两大突出进展：22详细课资一是大量用于设计工作的实验被数学模型的研究逐步取代，如航天飞机设计、反应堆设计、人工心瓣膜设计等；二是能获取和存储大量的数据，并能提取隐秘的信息，如计算机层析X射线摄影，核磁共振等。23详细课资 3战略计算u“战略计算”一词首

7、次出现在1995年美国为了确保核库存的性能、安全性、可靠性和更新需要而实施的“加速战略计算创新（ASCI）计划”。24详细课资这是因为美国克林顿总统在1995年8月11日宣布：“美国决定谋求真正的“零当量”全面禁止试验核武器条约”。25详细课资这并不意味着核竞赛的结束，恰恰相反是核武器计划新时代的开始，要求通过逼真的建模和模拟计算来取代传统的反复试验的工程处理方法。26详细课资这主要依赖于先进的数值计算和模拟能力，为此应用程序必须达到高分辩、三维、全物理和全系统的水平。27详细课资1995年8月22日（即美国总统宣布决定后的11天），能源部（DOE）就采购世界上最快的一台计算机（运算速度超过万

8、亿次）交付圣地亚实验室（96年12月安装）。28详细课资年月日，美国的DOE（Department of Energy ）/FNS共同联合组织召开了关于“先进科学计算”的全国会议。会议强调科学模拟的重要性，希望应用科学模拟来攻克复杂的科学与工程难题。29详细课资年月，美国在全国范围内倡议实施“科学模拟计划”(SSP)，提出要加速“燃烧系统”与全球气候系统“这两大应用领域的科学模拟研究。30详细课资 提问：提问：数值计算方法是做什么用的？数值计算方法是做什么用的？研究对象研究对象：数值问题数值问题有限个输入数据（问题的自有限个输入数据（问题的自变量、原始数据）与有限个输出数据（待求解数据）之变量

9、、原始数据）与有限个输出数据（待求解数据）之间函数关系的一个明确无歧义的描述。间函数关系的一个明确无歧义的描述。如一阶微分方程初值问题如一阶微分方程初值问题求函数解析表达式求函数解析表达式数学问题数学问题求函数求函数 在某些点在某些点的近似函数值的近似函数值数值问题数值问题31详细课资程序程序设计设计上机上机计算计算设计高效、可设计高效、可靠的数值方法靠的数值方法数值数值问题问题求解求解近似结果近似结果输出输出重点讨论重点讨论数值问题的来源：数值问题的来源：实际实际问题问题建立数学模型建立数学模型数值数值问题问题32详细课资数数值值方方法法的的设设计计原原则则收敛性：方法的可行性收敛性：方法的

10、可行性稳定性：初始数据等产生的误差对结果的影响稳定性：初始数据等产生的误差对结果的影响便于编程实现：逻辑复杂度要小便于编程实现：逻辑复杂度要小计算量要小：时间复杂度要小，运行时间要短计算量要小：时间复杂度要小，运行时间要短存贮量要尽量小：空间复杂度要小存贮量要尽量小：空间复杂度要小可可靠靠性性分分析析计计算算复复杂杂性性误差估计：运算结果不能产生太大的偏差且误差估计：运算结果不能产生太大的偏差且能够控制误差能够控制误差33详细课资 误误 差差 /* Error */一、一、 误差的来源与分类误差的来源与分类 /* Source & Classification */ 1、从实际问题中抽象出数学

11、模型从实际问题中抽象出数学模型 模型误差模型误差 /* Modeling Error */ 2、通过观测得到模型中某些参数（或物理量）的值通过观测得到模型中某些参数（或物理量）的值 观测误差观测误差 /* Measurement Error */ 3、数学模型与数值算法之间的误差求近似解数学模型与数值算法之间的误差求近似解 方法误差方法误差 (截断误差截断误差 /* Truncation Error */ ) 4、由于机器字长有限，原始数据和计算过程会产生新的误差由于机器字长有限，原始数据和计算过程会产生新的误差 舍入误差舍入误差 /* Roundoff Error */34详细课资二、二、

12、误差分析的基本概念误差分析的基本概念 /* Basic Concepts */ 设设 为真值（精确值），为真值（精确值）， 为为 的一个近似值的一个近似值 称称 为近似值为近似值 的绝对误差，简称误差。的绝对误差，简称误差。 注：注：误差可正可负，常常是无限位的误差可正可负，常常是无限位的绝对绝对误差限误差限/* accuracy */ 绝对值的上界绝对值的上界如：如：绝对绝对误差还不能完全表示近似值的好坏误差还不能完全表示近似值的好坏（绝对误差（绝对误差/* absolute error */）35详细课资近似值近似值 的误差的误差 与准确值与准确值 的比值：的比值：称为近似值称为近似值 的

13、相对误差，记的相对误差，记作作注：注：实际计算时，相对误差通常取实际计算时，相对误差通常取因为因为（相对误差（相对误差/* relative error */ ）36详细课资 一个算法如果输入数据有扰动（即误差），而计算一个算法如果输入数据有扰动（即误差），而计算过程中舍入误差不增长过程中舍入误差不增长, ,则称此算法是数值稳定的，否则则称此算法是数值稳定的，否则此算法就称为不稳定的。此算法就称为不稳定的。 （数值稳定性（数值稳定性/ /* Numerical Stability */） 对数学问题本身如果输入数据有微小扰动，引起对数学问题本身如果输入数据有微小扰动，引起输出数据（即问题真解）

14、的很大扰动，这就是病态问输出数据（即问题真解）的很大扰动，这就是病态问题。题。（病态问题（病态问题/ /* ill-posed problem */） 它是数学问题本身性质所决定的，与算法无关，它是数学问题本身性质所决定的，与算法无关，也就是说对病态问题，用任何算法（或方法）直接计也就是说对病态问题，用任何算法（或方法）直接计算都将产生不稳定性。算都将产生不稳定性。 37详细课资此公式此公式精确精确成立成立记为记为则初始误差则初始误差? ! !What happened?!例例 计算计算 公式一：公式一：38详细课资考察第考察第n步的误差步的误差我们有责任改变。我们有责任改变。造成这种情况的是

15、造成这种情况的是不稳定的算法不稳定的算法 /* unstable algorithm */迅速积累，误差呈递增趋势。迅速积累，误差呈递增趋势。初始的小扰动初始的小扰动 公式二：公式二：注意此公式与公式一注意此公式与公式一在理论上等价。在理论上等价。方法：先估计一个方法：先估计一个IN , ,再反推要求的再反推要求的In ( n N )。可取可取39详细课资取取40详细课资考察反推一步的误差：考察反推一步的误差：以此类推，对以此类推，对 n N 有：有：误差逐步递减误差逐步递减, 这样的算法称为这样的算法称为稳定的算法稳定的算法 /* stable algorithm */ 在我们今后的讨论中，

16、在我们今后的讨论中，误差误差将不可回避，将不可回避， 算法的算法的稳定性稳定性将会是一个非常重要的话题。将会是一个非常重要的话题。41详细课资例：蝴蝶效应例：蝴蝶效应 纽约的一只蝴蝶翅膀一拍，风和日纽约的一只蝴蝶翅膀一拍，风和日 丽的北京就刮起台风来了？！丽的北京就刮起台风来了？！纽约纽约北京北京这是一个这是一个病态问题病态问题42详细课资 蝴蝶效应蝴蝶效应在社会学界用来说明：一个坏的微小的机制，在社会学界用来说明：一个坏的微小的机制，如果不加以及时地引导、调节，会给社会带来非常大如果不加以及时地引导、调节，会给社会带来非常大的危害，戏称为的危害，戏称为“龙卷风龙卷风”或或“风暴风暴”；一个好

17、的微；一个好的微小的机制，只要正确指引，经过一段时间的努力，将小的机制，只要正确指引，经过一段时间的努力，将会产生轰动效应，称为会产生轰动效应，称为“革命革命”。 蝴蝶效应蝴蝶效应是气象学家洛伦兹是气象学家洛伦兹19631963年提出来的。其大意年提出来的。其大意为：一只南美洲亚马孙河流域热带雨林中的蝴蝶，偶为：一只南美洲亚马孙河流域热带雨林中的蝴蝶，偶尔扇动几下翅膀，可能在两周后引起美国德克萨斯引尔扇动几下翅膀，可能在两周后引起美国德克萨斯引 起一场龙卷风。起一场龙卷风。43详细课资几点注意事项几点注意事项 /* Remarks */1、 避免相近二数相减避免相近二数相减例：例：a1 = 0

18、.12345，a2 = 0.12346，各有，各有5位有效数字。位有效数字。 而而 a2 a1 = 0.00001，只剩下，只剩下1位有效数字位有效数字。 几种经验性避免方法：几种经验性避免方法：当当 | x | 1 时：时：44详细课资2、 避免小分母避免小分母 : 分母小会造成浮点溢出分母小会造成浮点溢出 /* over flow */3、避免大数避免大数吃吃小数小数例：例：用单精度计算用单精度计算 的根。的根。精确解为精确解为 算法算法1 1：利用求根公式利用求根公式在计算机内，在计算机内，109存为存为0.1 1010，1存为存为0.1 101。做加法时，做加法时，两加数的指数先向大指

19、数对齐，再将浮点部分相加。即两加数的指数先向大指数对齐，再将浮点部分相加。即1 的指数部分须变为的指数部分须变为1010，则：，则：1 = 0.0000000001 1010，取取单精度时就成为：单精度时就成为： 109+1=0.10000000 1010+0.00000000 1010=0.10000000 1010大数大数吃吃小数小数45详细课资?算法算法2：先解出先解出 再利用再利用注：注：注：注：求和时求和时从小到大从小到大相加，可使和的误差减小。相加，可使和的误差减小。例：按从小到大、以及从大到小的顺序分别计算例：按从小到大、以及从大到小的顺序分别计算1 + 2 + 3 + + 40

20、 + 1094、 先化简再计算，减少步骤，避免误差积累。先化简再计算，减少步骤，避免误差积累。一般来说，计算机处理下列运算的速度为一般来说，计算机处理下列运算的速度为5、选用稳定的算法。选用稳定的算法。46详细课资二.计算数学的两大分支47详细课资1.微分方程及其数值解 计算机解决实际问题的步骤建立数学模型选择数值方法编写程序上机计算 48详细课资现实世界中绝大多数事物的内外联系是及其现实世界中绝大多数事物的内外联系是及其复杂的，其状态随着时间、地点、条件的不复杂的，其状态随着时间、地点、条件的不同而不同，我们只能通过对问题进行简化和同而不同，我们只能通过对问题进行简化和作某些假定，从中找出其

21、状态和状态的变化作某些假定，从中找出其状态和状态的变化规律之间的关系，也即一个或一些函数与它规律之间的关系，也即一个或一些函数与它们的导数之间的关系，这种关系的数学表达们的导数之间的关系，这种关系的数学表达就是微分方程。就是微分方程。 49详细课资偏微分方程数值解主要是有限差分法和偏微分方程数值解主要是有限差分法和有限元法。有限元法。 50详细课资偏微分方程发展史： (1)十八世纪初, Taylor: (2)十九世纪中期，三类偏微分方程： (3)十九世纪末到二十世纪初，其它方程： 高阶方程： KDV方程： 51详细课资如果能找到一个（或一族）具有所要求阶连续导数的解析函数，如果能找到一个（或一

22、族）具有所要求阶连续导数的解析函数，将它代入微分方程（组）中，恰好使得方程（组）的所有条件将它代入微分方程（组）中，恰好使得方程（组）的所有条件都得到满足，我们就将它称为这个方程（组）的解析解（也称都得到满足，我们就将它称为这个方程（组）的解析解（也称古典解）。古典解）。“微分方程的真解微分方程的真解”或或“微分方程的解微分方程的解”就是指解就是指解析解。寻找解析解的过程称为求解微分方程。析解。寻找解析解的过程称为求解微分方程。 微分方程的解在数学意义上的存在性可以在非常一般的条微分方程的解在数学意义上的存在性可以在非常一般的条件下得到证明，这已有许多重要的结论。但从实际上讲，件下得到证明，这

23、已有许多重要的结论。但从实际上讲，人们需要并不是解在数学中的存在性，而是关心某个定义人们需要并不是解在数学中的存在性，而是关心某个定义范围内，对应某些特定的自变量的解的取值或是近似值范围内，对应某些特定的自变量的解的取值或是近似值-这这样一组数值称为这个微分方程在该范围内的数值解，寻找样一组数值称为这个微分方程在该范围内的数值解，寻找数值解的过程称为数值求解微分方程。数值解的过程称为数值求解微分方程。 52详细课资为什么要研究数值求解方法呢？为什么要研究数值求解方法呢？ 1）在实际问题中我们所能获取的或感兴趣的，往往只是一个）在实际问题中我们所能获取的或感兴趣的，往往只是一个特定点上的数据。如

24、空间的温度分布只能一个点一个点地测定，特定点上的数据。如空间的温度分布只能一个点一个点地测定，火箭升空传回的控制信息只能以某个确定的时间为间隔，一个火箭升空传回的控制信息只能以某个确定的时间为间隔，一个个地发送和接受，如此等等。这些离散点上的函数值对于解决个地发送和接受，如此等等。这些离散点上的函数值对于解决实际问题，已经足够了，寻找解析解的一般形式未必必要。实际问题，已经足够了，寻找解析解的一般形式未必必要。 53详细课资2）在很多情况下，寻找解析解也并无可能。现实问题中归结的）在很多情况下，寻找解析解也并无可能。现实问题中归结的微分方程不满足解析解的存在条件的比比皆是，方程中出现的微分方程

25、不满足解析解的存在条件的比比皆是，方程中出现的有些函数连续性都无法保证，它们并不存在前述意义的解析解。有些函数连续性都无法保证，它们并不存在前述意义的解析解。于是，求数值解便成了在这种情况下解决问题的重要手段了。于是，求数值解便成了在这种情况下解决问题的重要手段了。 54详细课资3）即使微分方程的解析解存在，以并不意味可以将它表示为初）即使微分方程的解析解存在，以并不意味可以将它表示为初等函数，如多项式、对数函数、指数函数三角函数及它们的不等函数，如多项式、对数函数、指数函数三角函数及它们的不定积分的有限组合形式定积分的有限组合形式显式解。显式解。 事实上，有显式解的微分方程只占解析解存在的微

26、分方程中事实上，有显式解的微分方程只占解析解存在的微分方程中的非常小的一部分。的非常小的一部分。 55详细课资56详细课资由于数字电子计算机只能存储有限个数据和作有限次运算，所以任何一种适用于计算机解题的方法，都必须把连续问题离散化，最终化成有限形式的代数方程组。 57详细课资数值求解方法概述数值求解方法概述1 1）区域剖分区域剖分 把整个定义域分成若干个小块，以便对每把整个定义域分成若干个小块，以便对每小块上的点或片求出近似值，这样按一定规律小块上的点或片求出近似值，这样按一定规律对定义域分切的过程称为区域剖分。对定义域分切的过程称为区域剖分。2 2）微分方程的离散微分方程的离散 区域剖分完

27、毕后，依据原来的微分方程去区域剖分完毕后，依据原来的微分方程去形成关于这些离散点或片的函数值的递推公式形成关于这些离散点或片的函数值的递推公式或方程。这是它们的未知量已不是一个连续函或方程。这是它们的未知量已不是一个连续函数，而成了若干个离散的未知值的某种组合了，数，而成了若干个离散的未知值的某种组合了，这个步骤称为微分方程离散。这个步骤称为微分方程离散。58详细课资3 3）初始和边界条件处理初始和边界条件处理 离散后系统是一个递推公式，那它需要若离散后系统是一个递推公式，那它需要若干个初值才能启动。若是一个方程组，那它所干个初值才能启动。若是一个方程组，那它所含的方程个数一般少于未知量的个数

28、，要想求含的方程个数一般少于未知量的个数，要想求解还需要补充若干个方程。这些需要补充的初解还需要补充若干个方程。这些需要补充的初值和方程往往可以通过微分方程的初始条件和值和方程往往可以通过微分方程的初始条件和边界条件来得到，这就是初始和边界条件处理边界条件来得到，这就是初始和边界条件处理过程。过程。 59详细课资4 4）离散系统的性态研究离散系统的性态研究 我们主要研究：这个系统是否可解，即解我们主要研究：这个系统是否可解，即解的存在性、唯一性问题；它与精确解的差距有的存在性、唯一性问题；它与精确解的差距有多大，这个差距当区域剖分的尺寸趋于零时，多大，这个差距当区域剖分的尺寸趋于零时，是否也会

29、趋于零，趋于零的速度多快，即解的是否也会趋于零，趋于零的速度多快，即解的收敛性和收敛速度问题；当外界对数据有所干收敛性和收敛速度问题；当外界对数据有所干扰时，所得的解是否会严重背离离散系统的固扰时，所得的解是否会严重背离离散系统的固有的解，即解的稳定性问题。有的解，即解的稳定性问题。60详细课资u上述问题说道底是一个上述问题说道底是一个误差误差分析分析问题，因为如果从实际问题，因为如果从实际问题到得出数值解的每一步问题到得出数值解的每一步都没有任何误差的话（当然，都没有任何误差的话（当然，这是不可能的），那么数值这是不可能的），那么数值解就应该是离散点上的精确解就应该是离散点上的精确值，也就不

30、用煞费苦心去讨值，也就不用煞费苦心去讨论上面的问题了。论上面的问题了。61详细课资数值求解微分方程过程示意数值求解微分方程过程示意微微分分方方程程区域剖分区域剖分离散系统的离散系统的 性态研究性态研究递推计算或解线递推计算或解线性代数方程组性代数方程组微分方程离散微分方程离散初始和边界条件处理初始和边界条件处理解的存在性、唯一性解的存在性、唯一性解的收敛性和收敛速度解的收敛性和收敛速度解的稳定性解的稳定性得到数值解得到数值解62详细课资求解过程中产生的误差求解过程中产生的误差现实问题现实问题数学模型数学模型离散格式离散格式模型模型误差误差建模建模离散离散舍入舍入误差误差观测观测模型模型截断截断

31、误差误差数值解数值解计算计算63详细课资2.数值代数u数值线性代数又称矩阵计算，它是科学与工程计算的核心。可以毫不夸张地讲，大部分科学与工程问题最终都要归结为一个矩阵计算问题，其中具有挑战性的问题是大规模矩阵计算问题。64详细课资数值线性代数主要包括如下三个大问题： （1）求解线性方程组的问题。（2）矩阵特征值问题。即给定一个矩阵 ，求它的部分或全部特征值以及对应的特征向量（3）线性最小二乘问题。 65详细课资求解线性方程组的问题。若 (1.1) 为相容方程组(即 bR(A) ) , 则一定有解 . 但当其解有无数个时 , 通常我们关心一个范数极小的解 .若 (1.1) 为矛盾方程组 , 通常

32、我们求解使残差向量 b-Ax 范数最小的那个解 , 即最小二乘解 .最常用的范数就是p-范数。若x=x1,x2,.,xn，那么xp=(|x1| p+|x2| p+.+|xn| p)1/p66详细课资u为了统一解决上述问题 , 广义逆理论应运而生 , 并最终解决了这些问题 .67详细课资u根据广义逆的相关理论 , 矛盾线性方程组的最小二乘解与 A 的1,3-逆有关 , 相容线性方程组的极小范数解与 A 的 1,4-逆有关 , 矛盾线性方程组的极小范数最小二乘解则与 A 的 M-P 广义逆有关 .68详细课资则此时的线性方程组会怎样?69详细课资注意对这种无限矩阵注意对这种无限矩阵,结合律不成立结

33、合律不成立!例例:70详细课资所以对这种无限区域的问题,有限差分,有限元方法都不在适用.通常我们会用数学的技巧将其转化,推出边界元或无限元方法.71详细课资数值方法和数值软件过去50年的主要进展附录 72详细课资 1940-1970ufloating point arithmetic; Fortran ; finite differences; finite elements; FFT; simplex algorithm; Monte Carlo; spline function73详细课资 1970-2000uquasi-Newton iterations; adaptivity; stiff ODE solvers; software libraries; Matlab; multigrid; sparse and iterative linear algebra; spectral methods; interior point methods74详细课资