精选优质文档-倾情为你奉上常微分方程的几种数值解法数学与应用数学 肖振华 指导教师 张秀艳【摘要】自然界与工程技术中的很多现象，可以归结为微分方程定解问题。其中，常微分方程求解是微分方程的重要基础内容。但是，对于许多的微分方程，往往很难得到甚至不存在精确的解析表达式，这时候，数值解提供了一个很好的解决思路。，针对于此，本文对常微分方程数值解法进行了简单研究，主要讨论了一些常用的数值解法，如欧拉法、改进的欧拉法、RungeKutta方法、Adams预估校正法以及勒让德谱方法等，通过具体的算例，结合MATLAB求解画图，初步给出了一般常微分方程数值解法的求解过程。同时，通过对各种方法的误差分析，让大家对各种方法的特点和适用范围有一个直观的感受。【关键词】 常微分方程 数值解法 MATLAB 误差分析【Abstract】 Many phenomena in nature and engineering can be attributed to the definite solution of the problem for differential equations. Among them, the ordinary differential equation solving is an important foundation for the content of the differential equations. However, many of the differential equations are often difficult to obtain accurate analytical expression .At this time, the numerical solution provides a good idea. For the Numerical Solution of Ordinary Differential Equations in this article, we focuses on some commonly used numerical solution, such as the Euler method, improved Euler method, Runge-Kutta method, Adams predictor corrector method as well as newer spectral methods. Through specific examples, combined with MATLAB solving and drawing, we initially know the solution process of general numerical solution of ordinary differential equations . At the same time, according to the error analysis of various methods , everyone has an intuitive feel of the characteristics and scope of the various methods.【Keywords】Ordinary Differential Equations Numerical Solution MATLAB error analysis目录1 前言1.1 常微分方程概述方程是一个在数学中非常熟悉的名词，在初等数学里，我们将我们要研究的问题作为一个或几个未知量，通过观察事物的规律，得出这些未知量与已知量之间的等式关系，这样就得到了一个简单的方程或方程组当然，这只是一个很浅显粗略的定义。在数学上，物质的运动和变化规律是通过函数关系来表示的，在一些复杂的现象中，比如火箭的运动轨迹，我们要求的未知量就变成了满足特定条件的一个或一些未知函数（不再是简单的数值）。有时候，我们需要用到导数很微分的知识，即这些未知函数的导数与自变量满足某种关系，这种方程，我们称之为微分方程，即可以这样定义，凡是含有参数，未知函数和未知函数导数（或微分）的方程称之为微分方程。未知函数是一元函数的微分方程称为常微分方程（ODE），未知函数是多元函数的微分方程称为偏微分方程（PDE），我们要研究的，就是常微分方程。1.2 常微分方程解与数值解法微分方程的解，就是找出一个代入方程使之成为恒等式的函数。.若微分方程的解中含有任意的个数与方程的阶数相同，且任意常数之间不能合并，则称此解为该方程的通解(或一般解).当通解中的各任意常数都取特定值时所得到的解，称为方程的特解。在实际问题中，这些函数往往还满足一些特定条件，这个称之为定解问题，如常微分方程的初值问题和边值问题，满足条件的解就是特解。在理想的状况下，我们解这个方程如果能够得出通解的解析表达式 ，那么代入数值很容易得到问题所要求的特解，但事实上，很多常微分方程得不到通解的解析表达式，或者表达式过于复杂。而且，常微分方程的特解是否存在，有几个解，牵涉到解的存在唯一性定理。在实际应用中，我们对解的精度要求往往并没有那么严格。在这种情况下，数值解应运而生，通俗的讲，数值解就是所求函数在特定点的函数值，通过函数在特定点的取值，我们可以反过来观察这个函数的某些性质，所以，求常微分方程的数值解是很有意义的，事实上，解常微分方程很大一部分就是解数值解。 下面我们要讨论几种常见的常微分方程数值解法，目的是展示常微分方程数值解法的一般过程。以下章节，由于本人水平有限，只对数值解法的思路和过程进行讨论，而对于解的存在唯一性，方法精度的具体证明，不做过多讨论。一般而言，构造数值解法都是基于两个想法：一种是基于数值积分的构造方法，一种是基于泰勒展开的构造法。下面我们要讨论的是几种常见的数值解法，如欧拉法、改进的欧拉法、龙格库塔法、阿达姆斯预估校正法以及基于勒让德多项式的谱方法。这些方法既有基于数值积分的构造方法，如欧拉法，也有基于泰勒展开的构造法，如龙格库塔法。下面的章节，我们会一一介绍这些方法的思路，并通过一些简单的算例，让大家了解这些解法的一般过程。为了让大家看到各种数值解法的精确性和有效性，我们可能会尽量选取一些通解存在的算例来求数值解，这样做并不意味着我们实际中求数值解时要先求出精确解，而只是出于对方法对比的需要，望不致引起混淆。2 欧拉法和改进的欧拉法2.1 欧拉法 2.1.1 欧拉法介绍首先，我们考虑如下的一阶常微分方程初值问题 （2-1）事实上，对于更复杂的常微分方程组或者高阶常微分方程，只需要将看做向量，（2-1）就成了一个一阶常微分方程组，而高阶常微分方程也可以通过降阶化成一个一阶常微分方程组。欧拉方法是解常微分方程初值问题最简单最古老的一种数值方法，其基本思路就是把（2-1）中的导数项用差商逼近，从而将一个微分方程转化为一个代数方程，以便求解。设在中取等距节点，因为在节点点上，由（2-1）可得： ， （2-2） 又由差商的定义可得： （2-3） 所以有 （2-4）用的近似值代入（2-4），则有计算的欧拉公式 （2-5）2.1.2欧拉法误差分析从欧拉公式中可以看出，右端的都是近似的，所以用它计算出来的会有累计误差，累计误差比较复杂，为简化分析，我们考虑局部截断误差，即认为是精确的前提下来估计，记为，泰勒展开有 （2-6）联立（2-5），（2-6）即得=+=，根据数值算法精度的定义，如果一个数值方法的局部截断误差=则称这个算法具有P阶精度，所以，欧拉方法具有一阶精度或者称欧拉方法为一阶方法。2.2 改进的欧拉方法2.2.1 改进的欧拉法介绍用数值积分离散化问题（1），两边做积分有： （2-7）对右端积分使用梯形积分公式可得： （2-8）同欧拉方法，用的近似值代入（2-7），联立（2-8）得到改进的欧拉方法计算公式： （2-9）一般来说，如果求解时，算法右端不包含，称为显性计算公式，如果包含，则求解时还需要解方程，这种称为隐式计算公式。显然公式（2-9）是一个隐式计算公式，事实上，改进的欧拉方法是用欧拉方法先求一个预测值，再用这个预测值来计算，即： （2-10）2.2.2改进的欧拉法误差分析同欧拉法误差分析类似，用泰勒展开容易知道改进的欧拉方法具有二阶精度，证明略。2.3 算例2.3.1（一阶常微分方程）求解初值问题 解析：在MATLAB中求解这个方程y=dsolve(Dy=y-2*x/y,y(0)=1,x) 得y =(2*x+1)(1/2) 它的解析解为，下面我们分别用欧拉方法和改进的欧拉方法来求其数值解。欧拉方法：创建M文件euler1.m,内容如下：function x,y=euler1(fun,x0,xfinal,y0,n)if nargin5,n=50;endh=(xfinal-x0)/n;x(1)=x0;y(1)=y0;for i=1:n x(i+1)=x(i)+h; y(i+1)=y(i)+h*feval(fun,x(i),y(i);End再定义函数方程组中的函数f1，创建f1.m文件，内容如下：function f=f1(x,y)f=y-2*x/y在MATLAB中输入：x,y=euler1(f1,0,1,1,20)输出f,x,y的值，将数值解跟精确解画图表示，输入：plot(x,y,r*-,x,sqrt(1+2*x),g+-)； xlabel(x);ylabel(y);title(y=y-2x/y);legend(数值解,精确解)得到图形，保存为euler1.fig，图形如下：改进的欧拉方法：创建M文件eulerprove1.m,内容如下：function x,y=eulerprove1(fun,x0,xfinal,y0,n)if nargin5,n=50;endh=(xfinal-x0)/n;y(1)=y0,x(1)=x0;for i=1:n x(i+1)=x(i)+h; y1=y(i)+h*feval(fun,x(i),y(i)/2; y2=h*feval(fun,x(i+1),y1