研究生数值分析-第6章 微分方程数值解法

《微分方程数值解法》【摘要】自然界与工程技术中的很多现象，可以归结为微分方程定解问题。

其中，常微分方程求解是微分方程的重要基础内容。

但是，对于许多的微分方程，往往很难得到甚至不存在精确的解析表达式，这时候，数值解提供了一个很好的解决思路。

，针对于此，本文对常微分方程数值解法进行了简单研究，主要讨论了一些常用的数值解法，如欧拉法、改进的欧拉法、Runge —Kutta 方法、Adams 预估校正法以及勒让德谱方法等，通过具体的算例，结合MA TLAB 求解画图，初步给出了一般常微分方程数值解法的求解过程。

同时，通过对各种方法的误差分析，让大家对各种方法的特点和适用范围有一个直观的感受。

【关键词】 常微分方程 数值解法 MA TLAB 误差分析引言在我国高校，《微分方程数值解法》作为对数学基础知识要求较高且应用非常广泛的一门课程，不仅在数学专业，其他的理工科专业的本科及研究生教育中开设这门课程．近四十年来，《微分方程数值解法》不论在理论上还是在方法上都获得了很大的发展．同时，由于微分方程是描述物理、化学和生物现象的数学模型基础，且它的一些最新应用已经扩展到经济、金融预测、图像处理及其他领域 在实际应用中，通过相应的微分方程模型解决具体问题，采用数值方法求得方程的近似解，使具体问题迎刃而解。

2 欧拉法和改进的欧拉法2.1 欧拉法2.1.1 欧拉法介绍首先，我们考虑如下的一阶常微分方程初值问题 ⎩⎨⎧==00)(),('y x y y x f y（2--1）事实上，对于更复杂的常微分方程组或者高阶常微分方程，只需要将x 看做向量，（2--1）就成了一个一阶常微分方程组，而高阶常微分方程也可以通过降阶化成一个一阶常微分方程组。

欧拉方法是解常微分方程初值问题最简单最古老的一种数值方法，其基本思路就是把（2--1）中的导数项'y 用差商逼近，从而将一个微分方程转化为一个代数方程，以便求解。

设在[]b a ,中取等距节点h ，因为在节点n x 点上，由（2--1）可得：))(()(',n n n x y x f x y =，（2--2）又由差商的定义可得：hx y x y x y n n n )()(('1-≈+）（2--3） 所以有 ))(,()()(1n n n n x y x hf x y x y +≈+ （2--4）用)(k x y 的近似值k y )1,(+=n n k 代入（2--4），则有计算1+n y 的欧拉公式))(,(1n n n n x y x hf y y +=+ （2--5）2.1.2欧拉法误差分析从欧拉公式中可以看出，右端的n y 都是近似的，所以用它计算出来的1+n y 会有累计误差，累计误差比较复杂，为简化分析，我们考虑局部截断误差，即认为n y 是精确的前提下来估计11)(++-n n y x y ，记为1+n ε，泰勒展开有)()()(''2)(')()(1321++<<+++=n n n n n x x h O y h x hy x y x y ξξ（2--6）联立（2--5），（2--6）即得1+n ε=)(''22ξy h +)(3h O =)(2h O ，根据数值算法精度的定义，如果一个数值方法的局部截断误差1+n ε=)(1+p h O 则称这个算法具有P 阶精度，所以，欧拉方法具有一阶精度或者称欧拉方法为一阶方法。

数值分析复习资料一、重点公式第一章 非线性方程和方程组的数值解法 1）二分法的基本原理，误差：~12k b ax α+--<2）迭代法收敛阶：1lim0i pi ic εε+→∞=≠，若1p =则要求01c <<3）单点迭代收敛定理：定理一：若当[],x a b ∈时，[](),x a b ϕ∈且'()1x l ϕ≤<，[],x a b ∀∈，则迭代格式收敛于唯一的根；定理二：设()x ϕ满足：①[],x a b ∈时，[](),x a b ϕ∈， ②[]121212,,, ()(),01x x a b x x l x x l ϕϕ∀∈-≤-<<有 则对任意初值[]0,x a b ∈迭代收敛，且：110111i i iii x x x llx x x lαα+-≤---≤-- 定理三：设()x ϕ在α的邻域内具有连续的一阶导数，且'()1ϕα<，则迭代格式具有局部收敛性；定理四：假设()x ϕ在根α的邻域内充分可导，则迭代格式1()i i x x ϕ+=是P 阶收敛的 ()()()0,1,,1,()0j P j P ϕαϕα==-≠ （Taylor 展开证明）4）Newton 迭代法：1'()()i i i i f x x x f x +=-，平方收敛 5）Newton 迭代法收敛定理：设()f x 在有根区间[],a b 上有二阶导数，且满足： ①：()()0f a f b <； ②：[]'()0,,f x x a b ≠∈；③：[]'',,f x a b ∈不变号④：初值[]0,x a b ∈使得''()()0f x f x <；则Newton 迭代法收敛于根α。

6）多点迭代法：1111111()()()()()()()()()i i i i i i i i i i i i i i i f x f x f x x x x x f x f x f x f x f x f x x x -+-----=-=+----收敛阶：P =7）Newton 迭代法求重根（收敛仍为线性收敛），对Newton 法进行修改 ①：已知根的重数r ，1'()()i i i i f x x x rf x +=-（平方收敛） ②：未知根的重数：1''()(),()()()i i i i u x f x x x u x u x f x +=-=，α为()f x 的重根，则α为()u x 的单根。

许多实际问题的数学模型是微分方程或微分方程的定解问题。

如物体运动、电路振荡、化学反映及生物群体的变化等。

常微分方程可分为线性、非线性、高阶方程与方程组等类；线性方程包含于非线性类中，高阶方程可化为一阶方程组。

若方程组中的所有未知量视作一个向量，则方程组可写成向量形式的单个方程。

因此研究一阶微分方程的初值问题⎪⎩⎪⎨⎧=≤≤=0)(),(y a y bx a y x f dxdy， (9-1) 的数值解法具有典型性。

常微分方程的解能用初等函数、特殊函数或它们的级数与积分表达的很少。

用解析方法只能求出线性常系数等特殊类型的方程的解。

对非线性方程来说，解析方法一般是无能为力的，即使某些解具有解析表达式，这个表达式也可能非常复杂而不便计算。

因此研究微分方程的数值解法是非常必要的。

只有保证问题(9-1)的解存在唯一的前提下，研究其数值解法或者说寻求其数值解才有意义。

由常微分方程的理论知，如果(9-1)中的),(y x f 满足条件（1）),(y x f 在区域} ),({+∞<<∞-≤≤=y b x a y x D ，上连续； （2）),(y x f 在上关于满足Lipschitz 条件，即存在常数，使得y y L y x f y x f -≤-),(),(则初值问题(9-1)在区间],[b a 上存在惟一的连续解)(x y y =。

在下面的讨论中，我们总假定方程满足以上两个条件。

所谓数值解法，就是求问题(9-1)的解)(x y y =在若干点b x x x x a N =<<<<= 210处的近似值),,2,1(N n y n =的方法。

),,2,1(N n y n =称为问题(9-1)的数值解，n n x x h -=+1称为由到1+n x 的步长。

今后如无特别说明，我们总假定步长为常量。

建立数值解法，首先要将微分方程离散化，一般采用以下几种方法： (1) 用差商近似导数在问题(9-1)中，若用向前差商hx y x y n n )()(1-+代替)(n x y '，则得)1,,1,0( ))(,()()(1-=≈-+N n x y x f hx y x y n n n n n)(n x y 用其近似值代替，所得结果作为)(1+n x y 的近似值，记为1+n y ，则有 1(,) (0,1,,1)n n n n y y hf x y n N +=+=-这样，问题(9-1)的近似解可通过求解下述问题100(,) (0,1,,1)()n n n n y y hf x y n N y y x +=+=-⎧⎨=⎩(9-2)得到，按式(9-2)由初值经过步迭代，可逐次算出N y y y ,,21。

微分方程数值解法微分方程数值解法微分方程数值解法【1】摘要：本文结合数例详细阐述了最基本的解决常微分方程初值问题的数值法，即Euler方法、改进Euler法，并进行了对比，总结了它们各自的优点和缺点，为我们深入探究微分方程的其他解法打下了坚实的基础。

关键词:常微分方程数值解法 Euler方法改进Euler法1、Euler方法由微分方程的相关概念可知，初值问题的解就是一条过点的积分曲线，并且在该曲线上任一点处的切线斜率等于函数的值。

根据数值解法的基本思想，我们取等距节点，其中h为步长，在点处，以为斜率作直线交直线于点。

如果步长比较小，那么所作直线与曲线的偏差不会太大，所以可用的近似值，即：，再从点出发，以为斜率作直线，作为的近似值，即：重复上述步骤，就能逐步求出准确解在各节点处的近似值。

一般地，若为的近似值，则过点以为斜率的直线为：从而的近似值为：此公式就是Euler公式。

因为Euler方法的思想是用折线近似代替曲线，所以Euler方法又称Euler折线法。

Euler方法是初值问题数值解中最简单的一种方法，由于它的精度不高，当步数增多时，由于误差的积累，用Euler方法作出的折线可能会越来越偏离曲线。

举例说明：解： ,精确解为：1.2 -0.96 -1 0.041.4 -0.84 -0.933 0.9331.6 -0.64 -0.8 0.161.8 -0.36 -0.6 0.242.0 0 -0.333 0.332.2 0.44 0 0.44通过上表可以比较明显地看出误差随着计算在积累。

2、改进Euler法方法构造在常微分方程初值问题 ,对其从到进行定积分得:用梯形公式将右端的定积分进行近似计算得：用和来分别代替和得计算格式:这就是改进的Euler法。

解：解得：由于 ,是线形函数可以从隐式格式中解出问题的精确解是误差0.2 2.421403 2.422222 0.000813 0.021400.4 2.891825 2.893827 0.00200 0.051830.6 3.422119 3.425789 0.00367 0.094112.0 10.38906 10.43878 0.04872 1.1973通过比较上表的第四列与第五列就能非常明显看出改进Euler方法精度比Euler方法精度高。

微分方程的数值解法微分方程是描述自然界中众多现象和规律的重要数学工具。

然而，许多微分方程是很难或者无法直接求解的，因此需要使用数值解法来近似求解。

本文将介绍几种常见的微分方程数值解法。

1. 欧拉方法欧拉方法是最简单的数值解法之一。

它将微分方程转化为差分方程，通过计算离散点上的导数来逼近原方程的解。

欧拉方法的基本思想是利用当前点的导数值来估计下一个点的函数值。

具体步骤如下：首先，将自变量区间等分为一系列的小区间。

然后，根据微分方程的初始条件，在起始点确定初始函数值。

接下来，根据导数的定义，计算每个小区间上函数值的斜率。

最后，根据初始函数值和斜率，递推计算得到每个小区间上的函数值。

2. 龙格-库塔方法龙格-库塔方法是一种常用的高阶精度数值解法。

它通过进行多次逼近和修正来提高近似解的准确性。

相比于欧拉方法，龙格-库塔方法在同样的步长下可以获得更精确的解。

具体步骤如下：首先，确定在每个小区间上的步长。

然后，根据微分方程的初始条件，在起始点确定初始函数值。

接下来，根据当前点的导数值，使用权重系数计算多个中间点的函数值。

最后，根据所有中间点的函数值，计算出当前点的函数值。

3. 改进欧拉方法（改进的欧拉-克罗默法）改进欧拉方法是一种中阶精度数值解法，介于欧拉方法和龙格-库塔方法之间。

它通过使用两公式递推来提高精度，并减少计算量。

改进欧拉方法相对于欧拉方法而言，增加了一个估计项，从而减小了局部截断误差。

具体步骤如下：首先，确定在每个小区间上的步长。

然后，根据微分方程的初始条件，在起始点确定初始函数值。

接下来，利用欧拉方法计算出中间点的函数值。

最后，利用中间点的函数值和斜率，计算出当前点的函数值。

总结：微分方程的数值解法为我们研究和解决实际问题提供了有力的工具。

本文介绍了欧拉方法、龙格-库塔方法和改进欧拉方法这几种常见的数值解法。

选择合适的数值解法取决于微分方程的性质以及对解的精确性要求。

在实际应用中，我们应该根据具体情况选择最合适的数值解法，并注意控制步长以尽可能减小误差。

第六章 常微分方程数值解法——RK 4法、AB 4法******（学号） *****（姓名）上机题目要求见教材P307,23题。

一、算法原理题目要求采用RK 4法和AB 4法求解最简单的常微分方程初值问题(,),()y f x y a x by a η'=≤≤⎧⎨=⎩ (1)为求解式(1)，采用离散化方法，就是寻求解)(x y 在区间],[b a 上的一系列点<<<<<n x x x x 321上的近似值 ,,,,21n y y y 。

记1(1,2,)i i i h x x i -=-=表示相邻两个节点的间距，称为步长。

求微分方程数值解的主要问题：(1) 如何将微分方程(,)y f x y '=离散化,并建立求其数值解的递推公式; (2) 递推公式的局部截断误差、数值数n y 与精确解)(n x y 的误差估计； (3) 递推公式的稳定性与收敛性. a) Runge-Kutta 方法基本思想：通过在1[,]i i x x +多预报几个点求斜率，并将其加权平均作为k *的近似值，以此构造更高精度的计算公式。

如果每步计算四次函数 的值，完全类似的，可以导出局部截断误差为)(5h O 的四阶Runge-Kutta 公式(RK 4)：1123412132431(22),6(,),(,),221(,),22(,).n n n n n n n n n n y y k k k k k f x y h h k f x y k h k f x h y k k f x h y hk +⎧=++++⎪⎪=⎪⎪⎪=++⎨⎪⎪=++⎪⎪=++⎪⎩ (2)b) Adams 显式公式Runge-Kutta 方法是单步法，计算1+n y 时，只用到n y ， 而已知信息1-n y 、2-n y 等没有被直接利用。

可以设想如果充分利用已知信息1-n y ，2-n y ，…来计算1+n y ，那么不但有可能提高精度，而且大大减少了计算量，这就是构造所谓线性多步法的基本思想。

微分方程是数学中的一种重要的方程类型，广泛应用于物理、工程、经济等领域。

解微分方程有各种方法，其中数值解法是一种重要而实用的方法。

微分方程的数值解法是通过数值计算来求解微分方程的近似解。

它的基本思想是将微分方程转化为差分方程，并用计算机进行迭代计算，从而求得微分方程的数值解。

数值解法的关键在于如何将微分方程转化为差分方程。

常见的方法有欧拉方法、改进欧拉方法、龙格－库塔方法等。

这些方法都是基于泰勒级数展开的原理进行推导的。

以欧拉方法为例，其基本思路是将微分方程中的导数用差商的方式近似表示，然后通过迭代计算，逐步逼近微分方程的解。

欧拉方法的具体步骤如下：首先确定微分方程的初始条件，即给定t0时刻的函数值y0，然后选取一定的步长ℎ，利用微分方程的导数计算差商y′=dy，进而根据差商dt得到下一个时刻的函数值y n+1=y n+ℎy′。

通过不断迭代计算，即可得到微分方程在一定时间区间内的数值解。

数值解法的另一个重要问题是误差控制。

由于数值计算本身的误差以及近似方法的误差，数值解法所得到的结果通常与真实解存在误差。

为了控制误差，常用的方法有缩小步长ℎ、提高近似方法的阶数等。

此外，还可以通过与解析解进行比较，评估数值解的准确性。

微分方程的数值解法具有以下几点优势。

首先，微分方程的解析解通常较难求得，而数值解法可以给出一个近似解，提供了一种有效的解决方案。

其次，数值解法可以利用计算机的高速运算能力，进行大规模复杂微分方程的求解。

此外，数值解法还可以在实际问题中进行仿真和优化，即通过调整参数来求解微分方程，从而得到最优解。

尽管微分方程的数值解法具有广泛的应用前景，但也存在一些问题和挑战。

首先，数值解法的稳定性和收敛性需要深入研究和分析。

其次，数值解法的计算量通常较大，对计算机运算能力和存储空间的要求较高。

此外，数值解法还需要对问题进行适当的离散化处理，从而可能引入一定的误差。

综上所述，“微分方程的数值解法”是一种重要而实用的方法，可以有效地求解微分方程的近似解。

微分方程数值解4.1当常微分方程能解析求解时，可利用Matlab符号工具箱中的功能找到精确解. 见下例求解方程,,,. 键入： yyy，，,20syms x y %定义符号变量diff_equ= ‘D2y+2*Dy-y=0’; %D2y表示,,,Dy= y,yy=dsolve (diff_equ, ‘x’) %定义x为自变量 y=cl*exp ((2^(1/2)-1)*x+c2*exp (-(2^(1/2)+1)*x)%表达式中含c1与c2，表示通解.%初始条件为y (0)=0,,y(0)=1时，按如下方式调用 y=dsolve (diff_equ,‘y (0)=0’, ‘Dy (0)=1’, ‘x’) y=1/4*2^(1/2)*exp ((2^(1/2)-1)*x)—1/4*2^(1/2)*exp (-(2^(1/2)+1)*x)%画出函数y=y (x)的图形ezplot (y,[-2,2])图形具体形式请上机试之.在方程无法获得解析解的情况下，可方便地获得数值解. 下面的例子说明用Matlab求数值解的方法及应注意的问题.[例1] 求解范德堡（vander pol）方程2dxdx2,,，,,(1)0xx 2dtdt求解高阶方程，必须等价地变换为一阶微分方程组，对本例，通过定义两个新的变量，实现这一变换yxydxdt1,2/,, 则令 dydty1/2,2dydtyyy2/(11)*21,,,,编写求解程序分为两部分，第一部分为待求解的方程，存盘的文件名为,待求解方程的函数名.m,，第二部分为求解主程序，本例中取名为main1.m.首先编写待求解方程的文件. 文件存盘名为“vdpol.m”. M,function yprime=vdpol(,)tyyprime (1)=y (2);; (1(1)^2)*(2)(1),,yyymu=2 yprime=[yprime (1);yprime (2)]; yprime (2)=mu*说明函数yprime=vdpol中. 定义为自变量，的形式取决于求解方程的阶数，本(,)tyyt 例中，,为解向量，为导数向量. yprime， y(2)yyyy,[(1),(2)],(1)(1)(1),y yprime,,，函数返回vander pol方程的导数列向量. 因为所求结果为方程数值解，(2)(1),y所以各向量维数只有在主程序求解时定下精度后才能确定.主程序定名为main1.m，你可用你所喜欢的其它名子，但vdpol.m除外. clear functions%调试程序时，放置这一语句是必要的. 它清除前边已编译的存在于内存中的废弃程序[]=ode23 (‘vdpol’,[0,30],[1,0]); ty,y1=y (:,1); %解曲线.y2=y (:,2); %解曲线的导数.polt ( ‘_ _’) tyty,1,,2,说明龙格_库塔的2阶与4阶改进型求解公式的实现，其指令分别为：[]=ode23 (‘f’,tsx,0,options) tx,[]=ode45 (‘f’,tsx,0,options) tx,其中可由系统依据精度要求自动设定，亦可由使用者依据实际需要自己确定，分别说明之. ts（1）若令tstttf,[0,1,,]，则输出在指定时刻tttf0,1,,给出，当tstktf,0::时，输出在区间[0,]ttf的等分点上给出，为步长. k（2）若tsttft,[0,],0为自变量初值，tf为终值，此时，options决定自变量的维数，t中的时间点不是等间隔的，这是为了保证所需的相对精度，积分算法改变了步长. 用于t,3,6设定误差限的参数options可缺省，此时系统设定相对误差为，绝对误差为，若1010自行设定误差限，可用如下语句：options=odeset (‘reltol’,, ‘abstol’,) rtat这里的与分别为设定的相对与绝对误差. rtat须注意的是无论用哪种方法确定ttf0,的取值方式，必须由使用者确定且应与相匹配. x0t,y0[1,0],ts,[0,30]y(0)1,y(0)0,为初始条件，本例中，因为，这意味着解曲线，，x0一般说，当解nnn个未知函数的方程组时，为维向量，共含有个初始条件. x0两个输出参数是列向量xx与矩阵，它们具有相同的行数，而矩阵的列数等于方程t组的个数，本例中y(:,1)y(:,2)的列数为2，其中，为自变量上各点函数值，为上各ytt点导数值.最后，提请读者注意的是：ode45也不总是比ode23好，在很多时候，低阶算法更有效，有关微分方程数值解法的更进一步信息，请参考数值分析方面的书籍. 有些参考书提供了一些关于算法选择和如何处理那些时间常数变化范围大的病态方程的非常实用的算法. 4.2 -设有一阶方程与初始条件,yfxy,(,), （4.1） ,yxy(),00,其中适当光滑，关于满足Lipschitz条件，即存在使 fxy(,)LyfxyfxyLyy(,)(,),,,1212则（4.1）式的解存在且惟一.关于yyx,()的解析解一般难以求到或根本无解析解，因而，实际问题中，通常，采用差分的方法. 在一系列离散点xxx,,,,yyy,,,,上寻求其数值近似解. 12n12n相邻两个节点间的间距xxnhn,，,,1,2,hxx,,称为步长，一般地取等步长，则hn0nnn，11、欧拉方法在区间[,]xx上用差商 nn，1yxyx()(),nn，1 h代替（4.1）式中,[,]xxxxxy，对fxy(,)中在上取值还是，而形成向前欧拉公式nn，1nn，1与向后欧拉公式.（1）向前欧拉公式xfxy(,)取左端点，得如下公式 nyxyxhfxyx()()(,()),，（4.2） nnnn，1从yxy(),x点出发，由初值代入（4.2）求得 000yyhfxy,，(,) （4.3） 1000反复利用（4.2），有yyhfxyn,，,(,) 0,1,2, （4.4） nnnn，1其几何意义如图4.1所示.y 图中yyx,()为方程（4.1）的精确 P P43P 2解曲线，其上任意点(,)xy处切线斜率为误差 P 1yyx,() 32Pxy(,)fxy(,). 从初值点出发，用该 P000 0y 0点斜率fxy(,)xx,作一直线段，在 001yyx,() yx() 3处得到Pxy(,)y，由（4.2）式确定， 1111y 3再从Pxy(,)fxy(,)出发，以为斜率 11111作直线段，在xx,Pxy(,)处得到， 2222xxxxx x O03124PPP， 012作为积分曲线yx()的近似，用图4.1 yyx,()n，1这一过程继续下去，形成折线表示在xy处的精确值，为解的近似值，不难得到 n，1n，12h32,,()()()()yxyyxOhOh,,，, nnn，，112P，1这一误差称为局部截断误差. 若一种算法局部截断误差为Oh()，则称该算法具有阶P精度，所以向前欧拉公式具有1阶精度.（2）向后欧拉公式若[,]xxxx中取中的，则有如下公式： fxy(,)nn，1n，1yyhfxyn,，,(,) 0,1,2, （4.5） nnnn，，，111称式（4.5）为向后欧拉公式，因为此式中y未知，故称其为隐式公式，无法用其直n，1接计算y，一般用向前欧拉公式产生初值. n，1(0)yyhfxyn,，,(,) 0,1,2, 11nnnn，，再按下式迭代(1)()kk，yyhfxykn,，,,(,),0,1,,0,1, nnnn，，，111其误差估计如下2h32,,()()()()yxyyxohoh,,，, nnn，，112精度亦为1阶，将向前欧拉公式（4.4）与向后欧拉公式（4.5）及它们的误差的几何说明作一对比，是十分有益的，见图4.2.y 为讨论局部截断误差，在图4.2中设点APxy(,)落在积分曲线yyx,()上，按式 nnnyyx,() （4.4）及式（4.5）分别得 ,P点为与， ABn，1 B且P AB,yyx,()点一定在积分曲线上相应 n点的上、下两边,所以将式（4.4）与（4.5） , 平均之，一定能得到更好的结果.xxx （3）梯形公式 nn，1 将向前与向后欧拉公式加以平均得到所图4.2 谓梯形公式hyyfxyfxyn,，，,[(,)(,)] 0,1,2, （4.6） nnnnnn，，，11123其局部截断误差为Oh()，具有2阶精度.（4）改进的欧拉公式为使计算简单，又免去迭代的繁复，将公式（4.6）简化为两步yyhfxy,，(,)nnnn，1h （4.7） yyfxyfxyn,，，,[(,)(,)], 0,1,2,nnnnn，，11n，12或写为h,yykk,，，()nn，112,2,1nn （4.8） n,0,1,2,kfxy,(,),,211nn，kfxyhk,，(,),,最后指出，上述欧拉方法可推广至微分方程组，如,yfxyz,(,,),,,zgxyz,(,,), ,yxy(),00,,zxz(),,00向前欧拉公式为yyhfxyz,，(,,),nnnnn，1 n,0,1,2, ,zzhgxyz,，(,,),nnnnn，12、龙格_库塔方法由微分中值定理,[()()]/(),01yxyxhyxh,,，,,,, nnn，1又因为,,yxhfxhyxh()(,())，,，，,,,yfxy,(,)，所以 nnn从而有yxyxhfxhyxh()()(),(),，，，,, （4.9） nnnn，1令[,]xx，称其为区间上的平均斜率，由（4.9）可知，给kfxhyxh,，，(,()),,nn，1nn出一种平均斜率，可相应导出一种算法. 向前欧拉公式中，精度低. 改进欧kfxy,(,)nn1拉公式中取[,]xxkfxyfxy,，((,)(,))，精度提高，下面，我们在区间内nn，1nnn，1n，12多取几个点，将其斜率加权平均，就能构造出精度更高的计算公式，公式的推导不再具体给出，只开列具体结果.（1）2阶龙格_库塔公式yyhkk,，，(),,,nn，11122,kfxy,(,) （4.10） 1nn,,21nnkfxahyhka,，，,,(,),0,1,,,1,其中,,,，,,,1,,1a，由于4个未知数只有3个方程，所以解不惟一，若令1222a1，即得改进的欧拉公式，具有2阶精度. ,,,，,,,,1a122（3）4阶龙格_库塔公式只给出精典格式中最常用的一种.h,yykkkk,，，，，(22)nn，11234,6,kfxy,(,),1nn,hh, （4.11） kfxyk,，，(,),21nn22,hh,kfxyk,，，(,)32nn,22,kfxhyhk,，，(,),43nn,其计算精度为4阶4.31、模型与问题[例2] 单摆运动图4.3中一根长的细线，一端固定，另一端悬挂质量为 lm的小球，在重力作用下，小球处于竖直的平衡位置. 现使小球偏离平衡位置一个小的角度，然后使其自由运动，在不 ,考虑空气阻力情形下，小球将沿弧线作周期一定的简谐运动.为平衡位置，在小球摆动过程中，当与平衡位置夹 ,,0角为,mgsin,时，小球所受重力在其动运轨迹的分量为 , , l（负号表示力的方向使减少），由牛顿第二定律可得微分方 ,程,,mltmg,,()sin,, （4.12）设小球初始偏离角度为,，且初速为0，式（4.12）的初 0,始条件为,,,,(0),(0)0,, （4.13） 0 mg 当,不大时，，式（4.12）化为线性常系数微sin,,,0分方程图4.3g,,,,，,0 （4.14） l解得g,,()costt, （4.15） 0ll简谐运动的周期为. T,2,g现在的问题是：当,较大时，仍用近似，误差太大，式（4.12）又无解析解，,sin,0试用数值方法在,,30,10,,两种情况下求解，画出的图形，与近似解（4.15）,()t00比较，这里设. l,25cm[例3] 捕食与被捕食当鲨鱼捕食小鱼，简记为乙捕食甲，在时刻，小鱼的数量为，鲨鱼的数量为，xt()yt()t当甲独立生存时它的（相对）增长率与种群数量成正比，即有,r，为增长率，xtrxt()(),而乙的存在使甲的增长率r减少，设减少率与乙的数量成正比，而得微分方程,xtxtraytrxaxy()()(()),,,, （4.16）比例系数a反映捕食者掠取食饵的能力.乙离开甲无法生存，设乙独自存在时死亡率为,，ytdyt()(),,，甲为乙提供食物，d使乙的死亡率降低，而促其数量增长，这一作用与甲的数量成正比，于是yt()满足 d,ytydbxtdybxy()(()),,，,,，（4.17）比例系数反映甲对乙的供养能力，设若甲，乙的初始数量分别为 bxxyy(0);(0),, （4.18） 00则微分方程（4.16），（4.17）及初始条件（4.18）确定了甲，乙数量xt()、yt()随时间变化而演变的过程，但该方程无解析解，试用数值解讨论以下问题：（1）设rdabxy,,,,,,1,0.5,0.1,0.02,25,2，求方程（4.16），（4.17）在00条件（4.18）下的数值解，画出xtyt(),()的图形及相图(,)xy，观察解xtyt(),()的周期变化，近似确定解的周期和的最大、小值，近似计算在一个周期内的平均值. xy,xy,（2）从式（4.16）和（4.17）消去得到 dtdxxray(),, （4.19） dyydbx(),，解方程（4.19），得到的解即为相轨线，说明这是封闭曲线，即解确为周期函数.（3）将方程（4.17）改写为,1yxtd()[],，（4.20） by在一个周期内积分，得到xt()yt()一周期内的平均值，类似可得一周期内的平均值，将近似计算的结果与理论值比较.2、实验（1）方程（单摆问题）,,mlmg,,,,sin, ,,,,,(0),(0)0,,0,无解析解，为求其数值解，先将其化为等价的一阶方程组. 令,xx,,,,,，方程化为 12,,xx,12,g,,21,,xxsin ,l,102,,xx(0),(0)0,,,其中，gl,,9.8,25,,为（弧度）与（弧度）两种情况，具100.1745,300.5236,0体编程如下：先以danbai.m为文件名存放待解方程. 键入： function xdot =danbai (t,x) %x=[x (1),x (2)],g=9.8;1=25; %x (1)为解向量, x (2)是导数 xdot (1)=x (2); %xdot(1)=,,(1)= x,xdot (2)=-g/1*sin (x (1)); %xdot (2)=,, ,xdot=[xdot (1);xdot (2)]; %必须将导数向量以列向量形式给出.再以主程序(求数值解)调用待求方程，主程序用main2.m为文件名存盘，其代码如下：clear functions[t,x]=ode23 (‘danbai’,[0,10],[0.1745,0])%只计算,,100.1745,（弧度）的情形.01g,,()cos,twtw,,%对近似解，周期 T2,,01gw=sqrt (9.8/25);y=0.1745*cos (w*t);[t,x (:,1),y] %显示数据，无分号. hold on %欲在同一幅图中画近似解. plot (t,x (:,1), ‘b’) %画数值解，绿色plot (t,y, ‘g*’) %用*号，红色画近似解. Hold off（2）食饵_捕食者方程,xtrxaxy(),,,ytdybxy(),,，可化为如下形式，，xrayx,0，，，，,,, ,,,,,,0,，dbxy，，，，y，，初始条件xxyy(0),(0),,表示为 00xx(0)，，，，0 ,,,,,yy(0)，，，，0以shier.m存盘如下代码function xdot=shier (t,x)r=1;d=0.5;a=0.1;b=0.02;xdot=diag ([r-a*x (2),-d+b*x (1)])*x;，，xxx(1)，，，，%x=,,,xdot= ,,,,,,,xy(2)，，，，y，，以main3.m存盘如下代码.clear functionsts=0:0.1:15;x0=[25,2];[t,x]=ode45 (‘shier’,ts,x0);[t,x] %显示数据t,x,y plot (t,x)gird %加网格线gtext (‘x(t)’),gtext (‘y(t)’), %用点鼠标方式pause,figure (2) %将x(t),x(t)放至指定点 12plot (x (:,1),x (:,2)) %以x (1)与x (2)为坐标点画相图xlabel (‘x’),ylabel (‘y’)可以猜测，xt()xt()(,)xx与是周期函数，相图是封闭曲线，从数值解可近似定出2121周期为10.7,x2.0,x的最大和最小值分别为与的最大和最小值分别为和，99.328.42.012为求xx与在一个周期的平均值，只需键入： 12y1=x (1:108,1); %x周期为10.7. 1x1p=trapz (y1)*0.1/10.7, %trapz (y1)返回按 y2=x (1:108,2); %梯形法对y1的积107yiyi1()1(1)，，x2p=trapz (y2)*0.1/10.7, %分值,x ,i12i,可得xx,,25,10 12对方程dxxray(),,，,dbxray,dxdy,化为 dyydbx(),，xy积分得解dbxray,,()()xeyec, （4.21）即为原方程组的相轨线，其中c由初始条件确定. 为说明上述相轨线是封闭的，令dbxray,,fxxegyye();(),,设其最大值分别为fg,xy,，若满足 mm00fxffyg(),(),, 00mmdr则有,,xy,fgc,,（令，可解出），又当时，相轨线（4.21）fg,,0,0xy,,,00mm00ba有意义. 当fgc,0,,cfgPxy(,)时，相轨线退化为一个点，对时，相轨线如mmmm00图4.4（c），而图（a），（b）分别为fx()与gy()的图像，下面讨论如何由（a），（b）画（c）.fxf(),nm y gyg(),gv()nm fx() y2 gm fm Pxy(,)qq 00 1 yP q02 q3 y1x y x x xxxy xyxy012001221 （a）（b）（c）图4.4设cpgpf,,,(0)yy,，xx,fxp(),，若令，则有，由（a）知，，使mm012qxy(,)xxx,,qxy(,)fxfxp()(),,，且于是相轨线应过与，对11010222012 fxgxpggyg()()(),,,xxx,(,)fxp(),gyq(),，，由，令，又由（b）知，mm12 存在yyy,,qxy(,)yy,gygyq()(),,qxy(,)使，且，于是相轨线又过与10231121242yyyyy,(),,xqq,两点，所以对间每个，轨线总要通过纵坐标为的两点，于是1210212相轨线是一条封闭曲线.相轨线封闭等价于xtyt(),()是周期函数，记周期为，为求其在一个周期内的平均T值(,)xy，由1yxtd()(),， by两边在一个周期内积分有((0)())yyT,：T11ln()ln(0)yTydTd,，，xxtdt,,，,() ,,,0TTbbb，，同理ry, a从而xxyy,,,00即xy的平均值恰为相轨线中心点的坐标，提请读者注意的是：这里的与与xtyt(),()p00初始条件中的xy,不是一件事. 00注意到在生态学上的意义，上述结果表明，捕食者的数量与食饵的增长率成rdab,,,正比，与它掠取食饵的能力成反比，食饵的数量与捕食者死亡率成正比，与它供养捕食者的能力成反比.3、练习内容（1）编写改进欧拉公式求微分方程数值解的程序，并用其与ode23求下列微分方程数值解，对二者作出比较.a）22,yxyy,,,,(0)0或y(0)1,.222 b),,,xyxyxny，，,,()02,,,21 yxsin,n,,（Bessel方程，这里令，其精确解为）. yy()2,(),,,x222,c),,,yyxyy，,,,cos0,(0)1,(0)0.（2）倒圆锥形容器，上底面直径为1.2m，容器的高亦为1.2m，在锥尖的地方开有一直径为3cm的小孔，容器装满水后，下方小孔开启，由水利学知识可知当水面高度为时，h水从小孔中流出的速度为为重力加速度，若孔口收缩系数为0.6（即若一个面vghg,2,积单位的小孔向外出水时，水柱截面积为0.6），问水从小孔中流完需多少时间？2分钟时，水面高度是多少？（3）一只小船渡过宽为的河流，目标是起点正对着的另一岸上点，已知河水ABd流速vv与船在静水中的速度之比为. k12（a）建立小船航线的方程，求其解析解.（b）设dvv,,,100m,1m/s,2m/s，用数值解法求渡河所需时间，任意时刻小12 船的位置及航行曲线，作图并与解析解比较.（c）若流速v为0,0.5,2 (m/s)，结果将如何？ 1（4）研究种群竞争模型. 当甲、乙两个种群各自生存时，数量演变服从下面规律xyxtrxytry()(1),()(1),,,, 12nn12其中，rr,nn,xtyt(),()分别为时刻甲，乙两个种群的数量，为其固有增长率，为它t1212们的最大容量，而当这两个种群在同一环境中生存时，由于乙消耗有限资源而对甲的增长产生影响，将甲的方程修改为xyxrxs,,,(1) （4.22） 11nn12这里s的含意是：对于供养甲的资源而言，单位数量乙（相对n）的消耗率为单位数12量甲（相对n）消耗的s倍，类似地，甲的存在亦影响乙的增长，乙的方程应改为 11xyytrys()(1),,, （4.23） 22nn12给定种群的初始值为xxyy(0),(0),, （4.24） 00及参数rrssnn,,,,,后，方程（4.22）与（4.23）确定了两种群的变化规律，因其解析121212解不存在，试用数值解法研究以下问题：（a）设rrnnssxy,,,,,,,,1,100,0.5,2,10，计算，画出xtyt(),()12121200 它们的图形及相图(,)xy，说明时间充分大以后xtyt(),()的变化趋势（人们今天看到的t已经是自然界长期演变的结果）.（b）改变rrxy,,,，但s与s不变（保持ss,,1,1），分析所得结果，若12001212ss,,,,1.5(1),0.7(1)，再分析结果，由此你得到什么结论，请用各参数生态学上的12含义作出解释.（c）试验当ss,,,,0.8(1),0.7(1)ss,,,,1.5(1),1.7(1)时会有什么结果；当1212时又会出现什么结果，能解释这些结果吗？。

微分方程数值解法微分方程是数学中的重要概念，它描述了物理系统中变量之间的关系。

解微分方程是许多科学领域中常见的问题，其中又可以分为解析解和数值解两种方法。

本文将重点介绍微分方程的数值解法，并详细讨论其中的常用方法和应用。

一、微分方程的数值解法概述微分方程的解析解往往较为复杂，难以直接求解。

在实际问题中，我们通常利用计算机进行数值计算，以获得方程的数值解。

数值解法的基本思想是将微分方程转化为一组离散的数值问题，通过逼近连续函数来获得数值解。

二、常见的数值解法1. 欧拉法欧拉法是最基础的数值解法之一，其核心思想是将微分方程转化为差分方程，通过逼近连续函数来获得数值解。

欧拉法的基本形式为：yn+1 = yn + h·f(xn, yn)其中，yn表示第n个时间步的数值解，h为时间步长，f为微分方程右端的函数。

欧拉法的精度较低，但计算简单，适用于初步估计或简单系统的求解。

2. 改进的欧拉法（Heun法）改进的欧拉法（Heun法）是对欧拉法的改进，其关键在于求解下一个时间步的近似值时，利用了两个斜率的平均值。

Heun法的基本形式为：yn+1 = yn + (h/2)·(k1 + k2)k1 = f(xn, yn)，k2 = f(xn+h, yn+h·k1)Heun法较欧拉法的精度更高，但计算量较大。

3. 龙格-库塔法（RK方法）龙格-库塔法是一类常用的数值解法，包含了多个不同阶数的方法。

其中，最常用的是经典四阶龙格-库塔法（RK4法），其基本形式为：k1 = f(xn, yn)k2 = f(xn + h/2, yn + (h/2)·k1)k3 = f(xn + h/2, yn + (h/2)·k2)k4 = f(xn + h, yn + h·k3)yn+1 = yn + (h/6)·(k1 + 2k2 + 2k3 + k4)RK4法实现较为复杂，但精度较高，适用于解决大多数常微分方程问题。

第一章 非线性方程和方程组的数值解法 1）二分法的基本原理，误差：~12k b ax α+--<2）迭代法收敛阶：1lim0i pi ic εε+→∞=≠，若1p =则要求01c <<3）单点迭代收敛定理：定理一：若当[],x a b ∈时，[](),x a b ϕ∈且'()1x l ϕ≤<，[],x a b ∀∈，则迭代格式收敛于唯一的根；定理二：设()x ϕ满足：①[],x a b ∈时，[](),x a b ϕ∈， ②[]121212,,, ()(),01x x a b x x l x x l ϕϕ∀∈-≤-<<有 则对任意初值[]0,x a b ∈迭代收敛，且：110111i i iii x x x ll x x x lαα+-≤---≤--定理三：设()x ϕ在α的邻域内具有连续的一阶导数，且'()1ϕα<，则迭代格式具有局部收敛性；定理四：假设()x ϕ在根α的邻域内充分可导，则迭代格式1()i i x x ϕ+=是P 阶收敛的 ()()()0,1,,1,()0j P j P ϕαϕα==-≠L （Taylor 展开证明）4）Newton 迭代法：1'()()i i i i f x x x f x +=-，平方收敛 5）Newton 迭代法收敛定理：设()f x 在有根区间[],a b 上有二阶导数，且满足： ①：()()0f a f b <； ②：[]'()0,,f x x a b ≠∈；③：[]'',,f x a b ∈不变号④：初值[]0,x a b ∈使得''()()0f x f x <；则Newton 迭代法收敛于根α。

6）多点迭代法：1111111()()()()()()()()()i i i i i i i i i i i i i i i f x f x f x x x x x f x f x f x f x f x f x x x -+-----=-=+----收敛阶：P =7）Newton 迭代法求重根（收敛仍为线性收敛），对Newton 法进行修改 ①：已知根的重数r ，1'()()i i i i f x x x rf x +=-（平方收敛） ②：未知根的重数：1''()(),()()()i i i i u x f x x x u x u x f x +=-=，α为()f x 的重根，则α为()u x 的单根。

第六章常微分方程的数值解法第六章常微分方程的数值解法在自然科学研究和工程技术领域中，常常会遇到常微分方程的求解问题。

传统的数学分析方法仅能给出一些简单的、常系数的、经典的线性方程的解析表达式，不能处理复杂的、变系数的、非线性方程，对于这些方面的问题，只能求诸于近似解法和数值解法。

而且在许多实际问题中，确确实实并不总是需要精确的解析解，往往只需获得近似的解或者解在若干个点上的数值即可。

在高等数学课程中介绍过的级数解法和逐步逼近法，能够给出解的近似表达式，这一类方法称为近似解法。

还有一类方法是通过计算机来求解微分方程的数值解，给出解在一些离散点上的近似值，这一类方法称作为数值方法。

本章主要介绍常微分方程初值问题的数值解法，包括Euler 方法、Runge-Kutta 方法、线性多步法以及微分方程组与高阶微分方程的数值解法。

同时，对于求解常微分方程的边值问题中比较常用的打靶法与有限差分法作了一个简单的介绍。

§1 基本概念1.1 常微分方程初值问题的一般提法常微分方程初值问题的一般提法是求解满足如下条件的函数，，b x a x y ≤≤)(=<<=α)(),(a y bx a y x f dxdy， (1.1) 其中),(y x f 是已知函数，α是给定的数值。

通常假定上面所给出的函数),(y x f 在给定的区域},),{(+∞<≤≤=yb x a y x D 上面满足如下条件：(1) 函数),(y x f 在区域D 上面连续；(2) 函数),(y x f 在区域D 上关于变量y 满足Lipschitz(李普希茨)条件：212121,),(),(y y b x a y y L y x f y x f ?≤≤?≤?，， (1.2)其中常数L 称为Lipschitz(李普希茨)常数。

由常微分方程的基本理论可以知道，假如(1.1)中的),(y x f 满足上面两个条件，则常微分方程初值问题(1.1)对于任意给定的初始值α都存在着唯一的解，，b x a x y ≤≤)(并且该唯一解在区间[a,b]上是连续可微的。