证明隐式Euler方法稳定性

第六章数值积分

6.1数值积分基本概念

6.1.1引言

在区间上求定积分

(6.1.1)

是一个具有广泛应用的古典问题，从理论上讲，计算定积分可用Newton-Leibniz公式

(6.1.2)

其中F(x)是被积函数f(x)的原函数.但实际上有很多被积函数找不到用解析式子表达的原函数，例如等等，表面看它们并不复杂，

但却无法求得F(x).此外，有的积分即使能找到F(x)表达式，但式子非常复杂，计算也很困难.还有的被积函数是列表函数，也无法用(6.1.2)的公式计算.而数值积分则只需计算f(x) 在节点xi(i=0，1,…,n)上的值，计算方便且适合于在计算机上机械地实现.

本章将介绍常用的数值积分公式及其误差估计、求积公式的代数精确度、收敛性和稳定性以及Romberg求积法与外推原理等.

6.1.2插值求积公式

根据定积分定义，对及都有

(极限存在)若不取极限，则积分I(f)可近似表示为

(6.1.3)

这里称为求积节点，与f无关，称为求积系数，(6.1.3)称为机械求积公式.

为了得到形如(6.1.3)的求积公式，可在上用Lagrange插值多项式

，则得

其中

(6.1.4)

这里求积系数由插值基函数积分得到，它与f(x)无关.如果求积公式(6.1.3)中的系数由(6.1.4)给出，则称(6.1.3)为插值求积公式.此时可由插值余项得到

(6.1.5)

这里ξ∈，(6.1.5)称为插值求积公式余项.

当n=1时，，此时

由(6.1.4)可得

于是

(6.1.6)

称为梯形公式.从几何上看它是梯形A bB(见图6-1)的面积近似曲线y=f(x)下的曲边梯形面积，公式(6.1.6)的余项为

(6.1.7)

6.1.3求积公式的代数精确度

当被积函数即f为次数不超过n的代数多项式时，，故

由(6.1.5)有，它表明插值求积公式(6.1.3)精确成立.对一般机械求积公式(6.1.3)，同样可以根据公式是否对m次多项式精确成立作为确定公式(6.1.3)中系数及节点的一种方法.在此先给出定义.

定义1.1一个求积公式(6.1.3)若对精确成立，而对不精确成立，则称求积公式(6.1.3)具有m次代数精确度.

根据定义，当时公式(6.1.3)精确成立，故有等式

(6.1.8)

而

(6.1.8)是关于系数及节点的方程组，当节点给定时，(6.1.8)取m=n就是关于系数的线性方程组，求此方程组就可求得求积系数.

例如n=1，取，求积公式为

在(6.1.8)中令m=1，可得

解得

它就是梯形公式(6.1.6)的系数，它与用公式(6.1.4)算出的结果完全一样.对梯

形公式(6.1.6)，当时

故求积公式(6.1.6)的代数精确度为一次.

对于具有(n+1)个节点的插值求积公式(6.1.3)，当时，故公式精确成立，它至少有n次代数精确度.反之，若求积公式(6.1.3)至少有n次代数精确度，则它是插值求积公式，即(6.1.3)的求积系数一定可用(6.1.4)求出.实际上，此时对求积公式(6.1.3)精确成立，若取f(x)为插值基函数，即

由(6.1.3)精确成立，可得

这就是(6.1.4)得到的插值求积公式系数.

定理1.1求积公式(6.1.3)是插值求积公式的充分必要条件是(6.1.3)至少具

有n次代数精确度.

定理表明直接利用代数精确度概念，由(6.1.8)可求得插值求积公式.更一般地，含有被积函数的导数的求积公式也同样可用代数精确度定义建立.如下例所示.

例6.1求积公式，已知其余项表达式为

.试确定系数及，使该求积公式具有尽可能高的代数精确度，并给出代数精确度的次数及求积公式余项.

解本题虽用到的值，但仍可用代数精确度定义确定参数及.令，分别代入求积公式.令公式两端相等，则

当

当

当

于是有

解得

，

此时，而上式右端为，两端不等，则求积公式对再令

，

不精确成立，故它的代数精确度为二次.

为求余项可将代入求积公式

当，代入上式得

，即

所以余项

.

6.1.4 求积公式的收敛性与稳定性

则称求积公式(6.1.3)是收敛的.

定义1.2若[

，

定义中n→∞包含了通常都要求用于计算积分(6.1.1)的求积公式(6.1.3)是收敛的.本章后面给出的求积公式都必须先证明其收敛性.

稳定性是研究计算和式

当有误差时，的误差是否增长.现设，误差为

.

定义1.3对任给，只要，就有

则称求积公式(6.1.3)是(数值)稳定的.

定义表明只要被积函数f(x)的误差充分小，积分和式的误差限就可任意小，则(6.1.3)就是数值稳定的.

定理1.2若求积公式(6.1.3)的系数则求积公式(6.1.3)是稳定的.

证明由于，故有

于是对，只要，就有

故求积公式(6.1.3)是稳定的.

讲解：

数值积分就是将求积分转化为求和

（即6.1.3式）这样不管被积函数多么复杂，它都能在计算机上机械实现。把（6.1.3）式称为机械求积公式，为求积节点，为

求积系数，建立求积公式有两种途径，一是利用的插值多项式积分得到，二是根据代数精确度概念，通过解方程得到及。特别当节点给定时，方程（6.1.8）是关于的线性方程组，它是容易求解的。

定理1.1给出了插值求积公式与代数精确度之间的关系。

求积公式收敛性简单说就是当和式收敛于积分值。而稳定

性是研究计算和式的误差积累，即当有误差时，只要误差充分小则和式误差也任意小，这就是稳定的。定理1.2表明只要求积公式（6.1.3）的系数，则求积公式就是稳定的。

6.2梯形公式与Simpson求积公式

6.2.1Newton-Cotes公式与Simpson公式

在插值求积公式中，若求积节点，此时可将求积公式写成

(6.2.1)

称为Newton-Cotes求积公式，其中系数称作Cotes系数.按(6.1.4)式引入变换，则有

(6.2.2)

由于是多项式积分，计算不会有困难.例如n=1时，

，

.这时求积公式就是我们熟悉的梯形公式(6.1.6).n=1～8时的系数见表6-1.从表中看到n=8时出现负数，稳定性没保证,所以一般只用n≤4的公式.

表6-1

当n=2时，由(6.2.2)可得，求积公式为

(6.2.3)

称为Simpson求积公式.

对梯形公式(6.1.6)，已知它的代数精确度为一次，且它的余项已由(6.1.7)给出，若记，则

由于，在上不变号，由积分中值定理得知，使

(6.2.4)

这就是梯形公式(6.1.6)的截断误差.

对Simpson公式(6.2.3)可以证明它的代数精确度是三次，根据定理1.1，

显然(6.2.3)对精确成立，再对，左端为

，右端为

故(6.2.3)对也精确成立.而对，公式(6.2.3)不精确成立.故求积公式(6.2.3)的代数精确度是三次，即(6.2.3)对任何精确成立.若令P(x)满足条件

这里,于是由(6.2.3)有

根据上章例5.6中(5.5.12)式有

上式两边积分，并记，则得

由于在区间上(不变号，故由积分中值定理得

于是有

(6.2.5)

这就是Simpson求积公式(6.2.3)的余项，即截断误差.

讲解：

当求积节点(k=0,1,…n)为等距节点，直接由插值求积得到的求积公式就转

为Newton-Cotes求积公式。但使用时通常只用n=1，2，4三种情况，它们分别称为梯形公式，Simpson公式和Cotes公式，梯形公式

它的局部截断误差直接由插值余项积分得到，由（6.2.4）给出，即

对于n=2的情形给出的Simpson求积公式（6.2.3），它具有3次代数精确，也就是它对任何次数不超3次的多项式精确成立，为使求积公式余项能反映这一性

质，我们不用2次的Lagrange插值近似被积函数f(x)，而用带导数的3点3次

Hermite插值，即造，满足条件：

于是由Hermite插值余项表达式，可得

对上式两端从到b积分，并利用积分中值定理就可得到Smpson求积公式的余项（6.25）对于n=4的Newton-Cotes公式根据（6.2.2）可算出，

，于是有

这里，余项

6.2.2复合梯形公式与复合Simpson公式

直接用梯形公式(6.1.6)及Simpson公式(6.2.3)计算积分I(f)误差较大，为达到精度要求，通常可将分为n个小区间，在小区间上应用梯形公式及Simpson

在每个小区间公式即可达到要求.为此取分点

,

上用梯形公式(6.1.6)，则得

或

(6.2.6)

称为复合梯形公式.根据定积分定义可知

故复合梯形公式(6.2.6)是收敛的，且(6.2.6)的求积系数，

故它也是稳定的.

(6.2.6)的截断误差可由(6.2.4)得到

若，根据连续函数性质，,使

于是得

它表明复合梯形公式(6.2.6)的截断误差阶为

如果在每个小区间上使用Simpson公式(6.2.3)，则得

(6.2.8) 称为复合Simpson公式，它的余项由(6.2.5)可得

即(6.2.9) 它表明.此外，还可证明

故复合Simpson公式(6.2.8)是收敛的，并且，故公式也是稳定的.

例6.2用n=8的复合梯形公式及n=4的复合Simpson公式，计算积分

，并估计误差.

解只要将区间［0，1］分为8等分，用公式(6.2.6)时取n=8，h=0.125，对

复合Simpson公式取n=4,h=0.25.计算各分点的函数值

.

由公式(6.2.6)及(6.2.8)得

为了估计误差，要求的高阶导数，由于

所以

故

由(6.2.7)得

对Simpson公式，由(6.2.9)得

例6.3 计算积分，若用复合梯形公式，问区间［0,1］应分多少等分才能使截断误差不超过?若改用复合Simpson公式，要达到同样精确度，区间［0，1］应分多少等分?

解本题只要根据及的余项表达式(6.2.7)及(6.2.9)即可求出其截断误差

应满足的精度.由于，对复合梯形公式

即.取n=213,即将区间［0,1］分为213等分时，用复

而用复合Simpson公式，则要求

合梯形公式计算误差不超过

.

取n=4，即将区间分为8等分，用n=4的复合即

.

Simpson公式即可达到精确度

.

讲解：

由于插值求积公式只能用n=1,2,4,用它算积分往往达不到精度要求，为提高求积精度我们可以采用将求积区间分为n个等距的小区间

在每个小区间上应用梯形公式或Simpson公式，这就得到

复合梯形公式及复合Simpson公式。它们的余项由每个小区间积分余项相加即可得到，结果分别由（6.2.7）及（6.2.9）给出。这两个求积公式是收敛和稳定的。如何利用公式计算积分近似值并估计误差可见例6.2及例6.3。

6.3外推原理与Romberg求积

6.3.1复合梯形公式递推化与节点加密

在计算机上用等距节点求积公式时，若精度不够可以逐步加密节点.设将区间分为n等分，节点，在区间上梯形公式为

若节点加密一倍，区间长为，记中点为在同一区间

上的复合梯形公式为

于是

(6.3.1)

它表明是在的基础上再加新节点的函数值之和乘新区间长，而

不必用(6.2.6)重新计算，这时有误差估计式

若，则得

(6.3.2)

这也是在计算机上估计梯形公式它表明用，其误差近似

.

误差的近似表达式.若(给定精度)，则.

若在区间［a,b］中做2n等分时，在上用Simpson公式计算，则由(6.2.8)

可知

它恰好是(6.3.2)中I(f)的近似值，即

它表明用(6.3.2)计算I(f)，其精度已由提高到

如果再将区间分半，使分为4个小区间，长度为，则

可由(6.3.1)计算出及，利用复合公式余项(6.2.9)得

如果，则有

(6.3.3)

从而有复合Simpson公式的误差估计

如果用(6.3.3)近似，即

(6.3.4)

则精度可达到.类似做法还可继续下去.这样对区间逐次分半，利

用公式(6.3.1)逐次递推.再由(6.3.2)，(6.3.3)逐次构造出精度愈来愈高的计算积分I(f)的公式，这就是Romberg求积的基本思想.

6.3.2外推法与Romberg求积公式

仍从梯形公式出发，区间中的节点如前所述.此时复合梯形公式可表示

为

当分为2n等分，区间长变为时，记，由于

此处将看作h的函数，将按的幂展开，可得到以下结果.

定理3.1设f在上的各阶导数存在，则复合梯形公式可展成

(6.3.5)

其中为不依赖h的常数.

定理证明见［2］.由(6.3.5)，显然有.在(6.3.5)中若用代替则得

(6.3.6)

用4乘以(6.3.6)减去(6.3.5)除以3，则得

若记

(6.3.7)

显然

(6.3.8)

具有精度.实际上就是复合Simpson公式中的为提高精度，可由及中消去，得

用精度为，如此逐次做下去，可得到

(6.3.9)

用时，精度为，这种将步长h逐次减半，使

逼近，以便精度逐次提高的方法称为外推法，它对于可展成h的幂级数的计算公式的加速收敛是很有效的，这里只将外推法用于计算积分

.

下面若用表示将区间二分k次得到的复合梯形公式，此时分

当k=0,1,…逐次得到即为等分的为等分，步长

，

复合梯形公式，加速一次得序列即为Simpson公式序列.加速m次则得，由(6.3.9)可将它表示为

(6.3.10)

称为Romberg求积公式.计算从k=0，即h=b-出发记，逐次二分得到T表(见表6-2).

当k增加时，先由(6.3.1)根据算出，再由(6.3.10)对m=1,2,…,k计算

.当f充分光滑时可证明

(T表任一列)

(T表对角线)

计算到(精度要求)为止.

例6.4用Romberg求积公式求的近似值，使其具有6位有效数字.

解本题直接用梯形递推公式(6.3.1)及Romberg求积公式(6.3.10)，按T表依次计算

其余计算结果见T表.

故计算停止，I≈0.946 083 1即为所求.

由于

，

例6.5证明等式

试依据的值，用外推算法求π的近似值.

解本题可利用Taylor展开式用外推原理求π的近似值.可令

.由Taylor公式展开得

若记，则其误差为

.

由外推法

其误差为

，

其误差为

，

根据以上公式计算结果如下表所示.

π=3.141 58即为所求.

讲解：

在等距节点的情况下，通过对求积区间的逐次分半，由梯形公式出可逐次提高求积公式精度，这就是Romberg求积的基本思路，由于梯形公式余项只有精度，即，但当节点加密时可组合成

四阶龙格库塔和向前欧拉方法和隐式欧拉法以及改进欧拉法教材

实验九欧拉方法 信息与计算科学金融崔振威201002034031 一、实验目的： 1、掌握欧拉算法设计及程序实现 二、实验内容： 1、p364-9.2.4、p386-9.5.6 三、实验要求： 主程序： 欧拉方法（前项）： function [x,y]=DEEuler(f,a,b,y0,n); %f:一阶常微分方程的一般表达式的右端函数 %a:自变量的取值下限 %b:自变量的取值上限 %y0：函数的初值 %n:积分的步数 if nargin<5,n=50; end h=(b-a)/n; x(1)=a;y(1)=y0; for i=1:n x(i+1)=x(i)+h; y(i+1)=y(i)+h*feval(f,x(i),y(i)); end format short 欧拉方法（后项）： function [x,y]=BAEuler(f,a,b,y0,n); %f:一阶常微分方程的一般表达式的右端函数 %a:自变量的取值下限 %b:自变量的取值上限 %y0：函数的初值 %n:积分的步数 if nargin<5,n=50; end h=(b-a)/n; x(1)=a;y(1)=y0; for i=1:n x(i+1)=x(i)+h; y1=y(i)+h*feval(f,x(i),y(i));

y(i+1)=y(i)+h*feval(f,x(i+1),y1); end format short 梯形算法： function [I,step,h2] = CombineTraprl(f,a,b,eps) %f 被积函数 %a,b 积分上下限 %eps 精度 %I 积分结果 %step 积分的子区间数 if(nargin ==3) eps=1.0e-4; end n=1; h=(b-a)/2; I1=0; I2=(subs(sym(f),findsym(sym(f)),a)+subs(sym(f),findsym(sym(f)),b))/h; while abs(I2-I1)>eps n=n+1; h=(b-a)/n; I1=I2; I2=0; for i=0:n-1 x=a+h*i; x1=x+h; I2=I2+(h/2)*(subs(sym(f),findsym(sym(f)),x)+subs(sym(f),findsym(sym(f)),x1)); end end I=I2; step=n; h2=(b-a)/n; 改进欧拉方法： function [x,y]=MoEuler(f,a,b,y0,n); %f:一阶常微分方程的一般表达式的右端函数 %a:自变量的取值下限 %b:自变量的取值上限 %y0：函数的初值 %n:积分的步数 if nargin<5,n=50; end

差分法欧拉格式浅谈

差分法欧拉格式浅谈 科学计算中常常求解常微分方程的定解这类问题的最简形式是一阶方程的初值问题 ???=='00 )(),(y x y y x f y （1） 这里假定右函数),(y x f 适当光滑，譬如关于y 满足Lipschitz 条件，以保证上述初值问题的解y(x)存在且唯一。 虽然求解常微分方程有各种各样的解析方法，但解析方法只能用来求解一些特殊类型的方程。求解从实际问题当中归结出来的微分方程主要靠数值解法。 差分方法是一类重要的数值解法。这类方法回避解y(x)的函数表达式，而是寻求它在一系列离散节点 <<<<

证明隐式Euler方法稳定性

第六章数值积分 6.1数值积分基本概念 6.1.1引言 在区间上求定积分 (6.1.1) 是一个具有广泛应用的古典问题，从理论上讲，计算定积分可用Newton-Leibniz公式 (6.1.2) 其中F(x)是被积函数f(x)的原函数.但实际上有很多被积函数找不到用解析式子表达的原函数，例如等等，表面看它们并不复杂， 但却无法求得F(x).此外，有的积分即使能找到F(x)表达式，但式子非常复杂，计算也很困难.还有的被积函数是列表函数，也无法用(6.1.2)的公式计算.而数值积分则只需计算f(x) 在节点xi(i=0，1,…,n)上的值，计算方便且适合于在计算机上机械地实现. 本章将介绍常用的数值积分公式及其误差估计、求积公式的代数精确度、收敛性和稳定性以及Romberg求积法与外推原理等. 6.1.2插值求积公式 根据定积分定义，对及都有 (极限存在)若不取极限，则积分I(f)可近似表示为 (6.1.3) 这里称为求积节点，与f无关，称为求积系数，(6.1.3)称为机械求积公式. 为了得到形如(6.1.3)的求积公式，可在上用Lagrange插值多项式 ，则得

其中 (6.1.4) 这里求积系数由插值基函数积分得到，它与f(x)无关.如果求积公式(6.1.3)中的系数由(6.1.4)给出，则称(6.1.3)为插值求积公式.此时可由插值余项得到 (6.1.5) 这里ξ∈，(6.1.5)称为插值求积公式余项. 当n=1时，，此时 由(6.1.4)可得 于是 (6.1.6) 称为梯形公式.从几何上看它是梯形A bB(见图6-1)的面积近似曲线y=f(x)下的曲边梯形面积，公式(6.1.6)的余项为 (6.1.7)

欧拉法仿真

仿真作业 1.普通欧拉公式法 clear; x1=0;x2=0; M=1.0;K=1.0;B=0.5; U=1.0; ST=20; T=0.1 NP=ST/T; for i=1:1:NP F1=x2; F2=(-K/M)*x1-(B/M)*x2+(1/M)*U; x1=x1+T*F1; x2=x2+T*F2; y(i)=x1; end plot((1:NP)*T,y,'m-') hold on 2.梯形公式法 clear; x10=0;x20=0; M=1.0;K=1.0;B=0.5; U=1.0; ST=50; T=0.1

NP=ST/T; for i=1:NP F10=x20; F20=(-K/M)*x10-(B/M)*x20+(1/M)*U; x11=x10+T*F10; x21=x20+T*F20; F11=x21; F21=-K*x11/M-B*x21/M+U/M; x12=x10+T*(F10+F11)/2; x22=x20+T*(F20+F21)/2; x10=x12; x20=x22; y(i)=x12; end plot((1:NP)*T,y,'r-') hold on 3.四阶龙格－库塔法（R-K法） clear; X10=0;X20=0; M=1;K=1;B=0.5;U=1; ST=30;T=0.1; NP=ST/T; for i=1:NP F10=X20;

F20=-K/M*X10-B/M*X20+U/M; X11=X10+T/2*F10; X21=X20+T/2*F20; F11=X21; F21=-K/M*X11-B/M*X21+U/M; X12=X10+T/2*F11; X22=X20+T/2*F21; F12=X22 F22=-K/M*X12-B/M*X22+U/M; X13=X10+T*F12; X23=X20+T*F22; F13=X22; F23=-K/M*X13-B/M*X23+U/M; F1=(F10+2*F11+2*F12+F13)/6; F2=(F20+2*F21+2*F22+F23)/6; X14=X10+T*F1; X24=X20+T*F2; X10=X14; X20=X24; Y(i)=X14; end plot((1:NP)*T,Y,'r-') hold on

常微分方程数值解法-欧拉方法(1)

课题报告 题目：常微分方程的数值解法-欧拉方法 院 (系)：理学院 专业：数学与信息专业 指导教师：岳宗敏 组员：艾佳欢(组长) 邓云娜 柏茜钟岩刘磊 2015 年 5 月 11 日

常微分方程数值解法-欧拉方法 摘要：从常微分方程数值解的基本概念入手，了解最基本的数值解法--欧拉方法。并利用 欧拉方法显式隐式的特点探究如何求解微分方程，以及欧拉方法的误差分析及校正。 关键词：数值解，欧拉方法，误差，校正 ABSTRACT: From the basic concept of numerical solution of ordinary differential equations, and understand the most basic numerical solution of euler method. And by using euler explicitly implicit characteristics and explore how to solve differential equations and the error analysis and correction of euler method. KEYWORDS ：arithmetic solution,Euler's method,error,revise 1.初值问题数值解基本概念 初值问题的数值解法，是通过微分方程离散化而给出解在某些节点上的近似值。 在[]b a ,上引入节点{}),,1(,:1100n k x x h b x x x a x k k k n n k k =-==<<<=-=称 为步长。在多数情况下，采用等步长，即),1,0(,n k kh a x n a b h k =+=-= 。记准确解为)(x y ,记)(k x y 的近似值为k y ,记),(k k y x f 为k f .一阶常微分方程的初值问题 ?? ?=∈=') 2.1()()() 1.1(),())(,()(0 x y a y b a x x y x f x y ， 若f 在{} ?+∞≤≤=y b x a D ,内连续，且满足Lip 条件：0≥?L 使 2121),(),(y y L y x f y x f -≤-，则初值问题的连续可微解)(x y 在],[b a 上唯一存在， 称解)(x y 在节点i x 处的近似值)(i i x y y =为其数值解，该方法称为数值方法。 2. Euler 方法 在工程计算中许多实际问题的数学模型都可以用常微分方程来描述。除少常系数线性微分方程和少数特殊的微分方程可用解析方法求解外，大多数常微分方程难以求得其精确解。因此研究常微分方程的数值解法具有重要的意义。本课题研究的是关于常微分方程初值问题的最简单的数值解法，单步法中的一种---Euler 方法。