精华内容
下载资源
问答
  • 几种常用数值积分方法的比较.doc
  • 它包括以下程序:梯形规则(用于数据和函数),辛普森规则(函数),梯形与辛普森规则... *图片由史蒂文·查普拉(Steven Chapra)和雷蒙德·卡纳尔(Raymond Canale)在其著作《工程师的数值方法》(第六版)中提供*
  • 数值积分方法

    万次阅读 2015-06-19 10:38:47
    数值积分是工程师和科学家经常使用的基本工具,用来计算无法解析求解的定积分的近似解。...那么我们就要通过数值积分方法来计算,数值积分的目的是,通过在有限个采样点上计算f(x) f (x)的值来逼近 f(x)f (x)在区间[a

    数值积分是工程师和科学家经常使用的基本工具,用来计算无法解析求解的定积分的近似解。
    如: Φ(x)=x0t3et1dt 不存在 Φ(x) 的解析解,要求 Φ(5)
    那么我们就要通过数值积分的方法来计算,数值积分的目的是,通过在有限个采样点上计算 f(x) 的值来逼近 f(x) 在区间 [a,b] 上的定积分.
    a=x0<x1<<xM=b . 称形如

    这里写图片描述

    且具有性质 baf(x)dx=Q[f]+E[f] 的公式为数值积分或 面积公式。项 E[f] 称为积分的截断误差,值 {xk}Mk=0 称为 面积节点 {wk}Mk=0 称为

    下面介绍几种常用的数值积分方法。

    基于多项式差值的面积公式

    通过M+1个等距点 {(xk,f(xk))}Mk=0 存在唯一的次数小于等于M的多项式 PM(x) 。当用该多项式来近似 [a,b] 上的 f(x) 时, PM(x) 的积分就近似等于 f(x) 的积分,这类公式称为牛顿-科特斯公式。当使用采样点 x0=axM=b 时,称为闭型牛顿-科特斯公式。

    这里写图片描述

    左/中/右矩形公式、梯形公式

    左/中/右矩形公式

    这里写图片描述

    梯形公式
    这里写图片描述

    图形如下
    这里写图片描述


    辛普森公式

    推导过程

    f(x) [a,b] 上有定义,将区间 [a,b] 分割成 n 等分(取n为偶数),则有 a=x0<x1<<xM=b ,其中

    xi=a+iΔxi=0,1,,n,Δx=ban

    这里我们想通过 (x0,f(x0)),(x1,f(x1)),(x2,f(x2)) 三点抛物线 g(x)=αx2+βx+γ 来取代 f(x) [x0,x2] 的定义,今儿求出它的近似积分值(如图),最后用累加的方式求得 f(x) [a,b] 上的近似积分。

    这里写图片描述

    由假设我们有
    这里写图片描述
    这里写图片描述

    所以可得
    baf(x)dx=x2x0+x4x2++xnxn2f(x)dx

    这里写图片描述

    误差估计

    若令 Sn=Δx3[f(x0)+4f(x1)+2f(x2)+4f(x3)+2f(xn2)+4f(xn1)+f(xn)] f(4)(x)[a,b] 上连续,则我们可以估计出辛普森公式的误差值为

    这里写图片描述

    例题1

    试用辛普森公式估计 10ex2dx ,取 n=6
    解:
    f(x)=ex2,Δx=16

    这里写图片描述


    拉格朗日插值

    分段线性插值

    例如:函数 f(x)=11+x2 如果在区间 [55] 上取 11 个等距节点: xk=5+k(k=0,1,2,...,10) ,由Lagrange插值公式可得到 f(x) 10 L10(x) 。如图所示:

    这里写图片描述

    L10(x) 仅在区间的中部能较好的逼近函数 f(x) ,在其它部位差异较大,而且越接近端点,逼近效果越差。
    可以证明,当节点无限加密时, Ln(x) 也只能在很小的范围内收敛,这一现象称为 Runge现象。它表明通过增加节点来提高逼近程度是不适宜的,因而不采用高次多项式插值。

    推导过程

    已知函数 f(x) 在区间 [xk,xk+1] 的端点上的函数值 yk=f(xk),yk+1=f(xk+1) ,求一个一次函数 y=P1(x) 使得 yk=f(xk),yk+1=f(xk+1) , 其几何意义是已知平面上两点 (xk,yk),(xk+1,yk+1) ,求一条直线过该已知两点。
    由直线的点斜式公式可知:

    P1(x)=yk+yk+1ykxk+1xk(xxk)

    把此式按照 ykyk+1 写成两项:
    P1(x)=xxk+1xkxk+1yk+xxkxk+1xkyk+1


    lk=xxk+1xkxk+1     lk+1=xxkxk+1xk

    并称它们为一次插值基函数。该基函数的特点如下:
    li(xk)={10k=ik!=i

    对于 i=0
    l0(x)=xx1x0x1,x[x0,x1]

    其它点上, l0(x)=0 ;
    对于 i=1,2,,n1
    li(x)={xxi1xixi1xxi+1xixi+1x[xi1,xi]x[xi,xi+1]

    其它点上, li(x)=0 ; 对于 i=n
    ln(x)=xxn1xnxn1,x[xn1,xn]

    其它点上, ln(x)=0 . 于是,
    P1(x)=k=0nyklk(x)

    此表达式与前面的表达式是相同的,这是因为在区间 [xk,xk+1] 上, 只有 lk(x),lk+1(x) 是非零的,其它基函数均为零。
    从而
    P1(x)=yklk(x)+yk+1lk+1(x)

    此形式称之为拉格朗日型插值多项式。其中, 插值基函数与 ykyk+1 无关,而由插值结点 xkxk+1 所决定。一次插值多项式是插值基函数的线性组合, 相应的组合系数是该点的函数值 ykyk+1 .

    这里写图片描述

    误差估计

    根据拉格朗日一次插值函数的余项,可以得到分段线性插值函数的插值误差估计: 对 x[a,b] ,当 x[xk,xk+1] 时,

    R(x)=12f(ξ)(xxk)(xxk+1)

    R(x)h28M ,其中
    h=max0kn1|xk+1xk|,M=maxx(a,b)f(x)

    于是有下面的定理:
    如果 f(x) [a,b] 上二阶连续可微,则分段连续函数 φ(x) 的余项有以下误差估计:
    |R(x)|=|f(x)ϕ(x)|h28M

    其中
    h=max0kn1|xk+1xk|,M=maxx(a,b)f(x)

    于是可以加密插值结点, 缩小插值区间, 使 h 减小, 从而减小插值误差。

    例题1

    已知函数y=f(x)=11+x2在区间 [0,5] 上取等距插值节点(如下表),
    求区间上分段线性插值函数,并利用它求出 f(4.5) 近似值。

    这里写图片描述

    解:
    在每个分段区间 [k,k+1] 上,
    P1(x)=x(k+1)k(k+1)yk+x(k)k+1(k)yk+1=yk(xk1)+yk+1(xk)

    于是
    P1(4.5)=0.05882×(4.55)+0.03846×(4.54)=0.04864

    拉格朗日型二次插值多项式

    推导过程

    已知函数 y=f(x) 在点 xk1,xk,xk+1 上的函数值 yk1=f(xk1),yk=f(xk),yk+1=f(xk+1) , 求一个次数不超过二次的多项式 P2(x) , 使其满足,

    P2(xk1)=yk1,P2(xk)=yk,P2(xk+1)=yk+1

    其几何意义为:
    已知平面上的三个点
    (xk1,yk1),(xk,yk),(xk+1,yk+1) ,
    求一个二次抛物线, 使得该抛物线经过这三点。

    插值基本多项式

    有三个插值结点 xk1,xk,xk+1 构造三个插值基本多项式,要求满足:
    (1) 基本多项式为二次多项式;
    (2) 它们的函数值满足下表:

    这里写图片描述

    因为 lk1(xk)=0,lk1(xk+1)=0 , 故有因子 (xxk)(xxk+1) , 而其已经是一个二次多项式, 仅相差一个常数倍, 可设
    lk1(x)=a(xxk)(xxk+1)

    又因为
    lk1(xk1)=1a(xk1xk)(xk1xk+1)=1

    得 
    a=1(xk1xk)(xk1xk+1)

    从而
    lk1(x)=(xxk)(xxk+1)(xk1xk)(xk1xk+1)

    同理得
    lk(x)=(xxk1)(xxk+1)(xkxk1)(xkxk+1)

    lk+1(x)=(xxk1)(xxk)(xk+1xk1)(xk+1xk)

    拉格朗日型二次插值多项式

    由前述, 拉格朗日型二次插值多项式:

    P2(x)=yk1lk1(x)+yklk(x)+yk+1lk+1(x)

    P2(x) 是三个二次插值多项式的线性组合,因而其是次数不超过二次的多项式,且满足:
    P2(xi)=yi,(i=k1,k,k+1)

    例题2

    已知:

    xi 101520
    yi=lgxi 11.17611.3010

    利用此三值的二次插值多项式求 lg(12) 的近似值。
    解:设 x0=10,x1=15,x2=20 ,则:


    这里写图片描述

    故:

    这里写图片描述
    这里写图片描述

    所以

    这里写图片描述

    误差估计

    我们在 [a,b] 上用多项式 Pn(x) 来近似代替函数 f(x) , 其截断误差记作 Rn(x)=f(x)Pn(x)
    x 在插值结点xi上时 Rn(xi)=f(xi)Pn(xi)=0 ,下面来估计截断误差:

    定理1:
    设函数 y=f(x) n 阶导数y(n)=f(n)(x) [a,b] 上连续,
    y(n+1)=f(n+1)(x) (a,b) 上存在;
    插值结点为: ax0<x1<<xnb ,
    Pn(x) n 次拉格朗日插值多项式;
    则对任意x[a,b]有:

    Rn(x)=1(n+1)!f(n+1)(ξ)ωn+1(x)

    其中 ξ(a,b) , ξ 依赖于x:ωn+1(x)=(xx0)(xx1)(xxn)
    maxaxb|fn+1(x)|Mn+1 ,则
    |Rn(x)|1(n+1)!Mn+1ωn+1(x)

    易知,线性插值的截断误差为:
    R1(x)=12f(ξ)(xx0)(xx1)

    二次插值的截断误差为:
    R2(x)=16f(3)(ξ)(xx0)(xx1)(xx2)

    例题3

    在例2中,用 lg10,lg15lg20 计算 lg12 .
    P2(12)=1.0766 ,
    e=|1.07921.0766|=0.0026
    估计误差:


    这里写图片描述


    三次样条插值法

    三次样条曲线原理

    x:a=x0<x1<<xM=by:y0   y1      yn

    定义

    样条曲线 S(x) 是一个分段定义的公式。给定 n+1 个数据点,共有 n 个区间,三次样条方程满足以下条件:

    a. 在每个分段区间[xi,xi+1]i=0,1,,n1x递增 S(x)=Si(x) 都是一个三次多项式。

    b. 满足 S(xi)=yii=0,1,,n

    c. S(x) ,导数 S(x) ,二阶导数 S(x) [a,b] 区间都是连续的,即 S(x) 曲线是光滑的

    所以n个三次多项式分段可以写作:

    Si(x)=ai+bi(xxi)+ci(xxi)2+di(xxi)3,i=0,1,,n1

    其中 ai,bi,ci,di 代表 4n 个未知系数。

    求解

    已知:
    a. n+1 个数据点 [xi,yi],i=0,1,,n
    b. 每一分段都是三次多项式函数曲线进行逼近
    c. 节点达到二阶连续
    d. 左右两端点处特性(自然边界,固定边界,非节点边界)

    根据定点,求出每段样条曲线方程中的系数,即可得到每段曲线的具体表达式。

    插值和连续性:
    这里写图片描述, 其中 i=0,1,,n1
    微分连续性:
    这里写图片描述, 其中 i=0,1,,n2
    样条曲线的微分式:
    image
    将步长 hi=xi+1xi 带入样条曲线的条件:
    a. 由 Si(xi)=yi(i=0,1,,n1) 推出

    ai=yi

    b. 由 Si(xi+1)=yi+1(i=0,1,,n1) 推出
    ai+hibi+h2ici+h3idi=yi+1

    c. 由 Si(xi+1)=Si+1(xi+1)(i=0,1,,n2) 推出
    Si(xi+1)=bi+2ci(xi+1xi)+3di(xi+1xi)2=bi+2cih+3dih2Si+1(xi+1)=bi+1+2ci(xi+1xi+1)+3di(xi+1xi+1)2=bi+1

    由此可得:
    bi+2hici+3h2idibi+1=0

    d. 由 Si(xi+1)=Si+1(xi+1)(i=0,1,,n2) 推出
    2ci+6hidi2ci+1=0

    mi=Si(xi)=2ci ,则
    a. 2ci+6hidici+1=0 可写为:
    mi+6hidimi+1=0 ,推出

    di=mi+1mi6hi

    b. 将 ci,di 带入 yi+hibi+h2ici+h3idi=yi+1 可得:
    bi=yi+1yihihi2mihi6(mi+1mi)

    c. 将 bi,ci,di 带入 bi+2hici+3hidi=bi+1(i=0,1,,n2) 可得:
    himi+2(hi+hi+1)mi+1+hi+1mi+2=6[yi+2yi+1hi+1yi+1yihi]

    端点条件

    由i的取值范围可知,共有n-1个公式, 但却有n+1个未知量m 。要想求解该方程组,还需另外两个式子。所以需要对两端点 x0 xn 的微分加些限制。 选择不是唯一的,3种比较常用的限制如下。

    a. 自由边界(Natural)

    首尾两端没有受到任何让它们弯曲的力,即 S=0 。具体表示为 m0=0 mn=0
    则要求解的方程组可写为:

    这里写图片描述
    这里写图片描述

    b. 固定边界(Clamped)

    首尾两端点的微分值是被指定的,这里分别定为A和B。则可以推出


    这里写图片描述
    这里写图片描述

    将上述两个公式带入方程组,新的方程组左侧为
    这里写图片描述

    c. 非节点边界(Not-A-Knot)

    指定样条曲线的三次微分匹配,即

    So(x1)=S1(x1)Sn2(xn1)=Sn1(xn1)

    根据 Si(x)=6di di=mi+1mi6hi ,则上述条件变为
    h1(m1m0)=h0(m2m1)hn1(mn1mn2)=hn2(mnmn1)

    新的方程组系数矩阵可写为:

    这里写图片描述

    右下图可以看出不同的端点边界对样条曲线的影响:


    这里写图片描述

    算法总结

    假定有n+1个数据节点

    (x0,y0),(x1,y1),(x2,y2),,(xn,yn),

    a. 计算步长 hi=xi+1xi(i=0,1,,n1)
    b. 将数据节点和指定的首位端点条件带入矩阵方程
    c. 解矩阵方程,求得二次微分值 mi
    d. 计算样条曲线的系数:

    这里写图片描述

    其中 i=0,1,,n1
    e. 在每个子区间 xixxi+1 中,创建方程
    gi(x)=ai+bi(xxi)+ci(xxi)2+di(xxi)3

    展开全文
  • 数值积分的基本思想由积分中值定理可知,在积分区间 内存在一点 ,成立 式的几何意义即为:底为 而高为 的矩形的面积恰等于所求曲边梯形的面积 。因此,要想求出 式左端积分,我们只需要知道三个值: 即可,这里 是...

    2f0ca1960975d55e47b79965e0312674.png

    数值积分的基本思想

    由积分中值定理可知,在积分区间

    内存在一点
    ,成立

    式的几何意义即为:底为
    而高为
    的矩形的面积恰等于所求曲边梯形的面积
    。因此,要想求出
    式左端积分,我们只需要知道三个值:
    即可,这里
    是显然的,问题在于
    的具体位置一般是不清楚的,从而
    未知。我们暂且将
    称为区间
    上的平均高度,我们的目标就是寻求一种求出平均高度
    的算法,这样我们的问题也就解决了。

    首先我们想到的是用区间两端点的“高度”

    的算术平均值作为平均高度
    的近似值,从而导出下述求积公式

    式即为我们我们熟知的
    梯形公式。写出MATLAB函数为
    function

    而若用区间中点

    的“高度”
    近似取代平均高度
    ,则可导出
    中矩形公式(简称 矩形公式

    其MATLAB函数为

    function

    总结:更一般地,我们可以在区间

    上适当选取某些节点
    ,然后用
    的加权平均得到平均高度
    的近似值,这样构造出来的求积公式具有下列通式

    其中

    称为
    求积节点
    称为
    求积系数,也称作伴随节点

    这种数值积分方法通常称为机械求积,其优势在于将积分求值问题归结为被积函数值的计算,很适合在计算机上使用。

    代数精度

    由于数值积分方法是近似方法,为保证精度,自然希望我们的求积公式能够对“尽可能多”的函数准确成立,这便提出了代数精度的概念。

    如果某个求积公式对于次数不超过
    的多项式均能够准确成立,但对于
    次多项式就不准确成立,则称该求积公式具有
    次代数精度
    (或 代数精确度)。

    前面提到的梯形公式

    和矩形公式
    均具有一次代数精度。

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

    在求积公式
    中,若
    其中
    ,则称求积公式
    收敛的。
    对任给
    ,若
    ,只要
    就有
    成立,其中
    为计算
    时产生误差
    后实际得到的值,即
    ,则称求积公式
    稳定的。

    定理

    若求积公式
    中系数
    ,则此求积公式是稳定的。

    牛顿-柯特斯(Newton-Cotes)公式

    将积分区间

    划分为
    等份,此时步长为
    ,选取等距节点
    构造出的插值型求积公式

    称为牛顿-柯特斯(Newton-Cotes)公式,其中

    称为
    柯特斯系数。可以通过下式确定

    • 时,
      这时的求积公式便是之前的梯形公式
    • 时,
      此时对应的求积公式便是
      辛普森(Simpson)公式:
      写成MATLAB函数为
    function
    • 时的牛顿-柯特斯公式则特别称为
      柯特斯公式,其形式为
      这里

    其MATLAB代码为

    function

    复合求积公式

    由于牛顿-柯特斯公式在

    时不具有稳定性,故不可能再通过提高阶的方法来提高求积精度。
    复合求积法便是通过把积分区间分成若干个子区间(通常是等分),再在每个子区间上使用低阶求积公式的方法来提高精度的。

    其实细心的同学可以发现,前面我所展示的MATLAB代码使用的便是对每个求积公式使用复合法完成的。常用的有复合梯形公式复合辛普森公式

    龙贝格(Romberg)求积公式

    前面介绍的复合求积方法可提高求积精度,如若精度仍不够,则可通过将步长逐次减半的方式来提高精度。如对复合梯形公式可导出其递推公式

    其中

    表示在
    基础上步长
    减半后的积分值。

    定理

    ,则有

    其中
    系数
    无关。

    式可以看出,
    阶,若用
    代替
    ,有

    再用

    式再减去
    式后再除以

    这里

    是与
    无关的系数。

    式可看出,用近似积分值
    ,其误差阶为
    ,显然误差阶是提高了。类似这种将计算
    的近似值的误差阶由
    提高到
    的方法称为
    外推算法,也称为 理查德森(Richardson)外推算法
    这是“数值分析”中一个重要的技巧,只要真值与近似值的误差能够表示成
    的幂级数,如
    式所示,都可以使用外推算法,提高精度。

    龙贝格(Romberg)算法

    次外推加速为

    余项为

    此方法常称为理查德森外推加速方法

    设用

    表示二分
    次后求得的梯形值,且以
    表示序列
    次加速值,则由递推公式
    可得

    上次则称为龙贝格求积算法,计算过程如下:

    1. ,求
      。令
      (
      记区间
      的二分次数)。
    2. 求梯形值
      ,即按递推公式
      计算
    3. 求加速值。
    4. (预先给定的精度),则终止计算,并且
      ;否则令
      继续计算。

    MATLAB代码为

    function

    自适应积分方法

    复合求积方法通常适用于被积函数变化不太大的积分,如果在求积区间中被积函数变化很大,有的部分函数值变化剧烈,另一部分却变化平缓。这时如果统一将区间等分再用复合求积公式计算积分将会导致计算量很大,我们想实现的是在满足误差要求的前提下,对变化剧烈部分将区间细分,而平缓部分则可使用大步长,也即针对被积函数在区间上不同情形采用不同的步长,使得再满足精度前提下积分计算工作量尽可能小。其算法技巧是在不同区间上预测被积函数变化的剧烈程度确定相应步长。这就是自适应积分方法的基本思想。

    下面是一个基于复合辛普森法的自适应积分算法的MATLAB代码:

    function

    下面我们计算积分

    ,可输入下列语句调用上述MATLAB函数进行计算
    clear
    

    结果为

    I 

    高斯求积公式

    下面研究带权积分

    为权函数)的求积公式

    下面看定义

    如果求积公式
    具有
    次代数精度,则称其节点
    高斯点,相应公式
    称为
    高斯型求积公式。

    这里主要包括四种高斯型求积公式,即

    • 高斯-勒让德求积公式
    • 高斯-切比雪夫求积公式
    • 高斯-拉盖尔求积公式
    • 高斯-埃尔米特求积公式

    多重积分

    这里主要说一下多重积分计算的思路,归结为一句话,就是:多重积分化累次积分,再由里到外使用数值积分公式进行求积计算。

    展开全文
  • VC++常用数值计算方法

    2008-11-21 10:03:46
    本收共不数值计算中常用的Visual C++子过程近200个,内容包括:解线性代数议程组、插值、数值积分、特殊函数、函数逼近、随机数、排序、特征值问题、数据拟合、方程求根和非线性方程组求解、函数的极值和最优化、...
  • 各种数值积分方法总结(龙贝格积分、高斯积分等)(一重积分的)常用数值积分方法牛顿-科茨(Newton-Cotes)积分公式梯形法则(2点积分)辛普森法则(3点积分和4点积分)辛普森1/3法则(3点积分)辛普森3/8法则(4...

    本文整理了各种类型的数值积分,从简单的梯形积分、辛普森积分到高精度的自适应积分、龙贝格积分和高斯积分;对于多重积分(高维积分),本文先简单介绍一下蒙特卡洛方法(适用于积分维数大于4的积分),后续会发帖专门针对多重积分方法做介绍。
    (如有疏漏,欢迎指正,谢谢~)

    1 (一重积分)常用的数值积分方法

    常用积分方法包括梯形积分、辛普森积分、自适应积分、龙贝格积分和高斯积分等;其中自适应积分、龙贝格积分和高斯积分属于高精度积分。

    1.1 牛顿-科茨(Newton-Cotes)积分公式

    牛顿-科茨(Newton-Cotes)积分公式:
    在这里插入图片描述
    包含了梯形法则、辛普森法则、布尔法则等。

    1.1.1 梯形法则(2点积分)

    梯形法则积分计算公式:

    其中,在这里插入图片描述
    在这里插入图片描述
    于是得:
    在这里插入图片描述
    梯形公式要计算两个点的函数值f(a)和f(b)。所以说是两点积分
    在这里插入图片描述
    梯形法则的误差为:
    在这里插入图片描述

    1.1.2-3 辛普森法则(3点积分和4点积分)

    辛普森法则分为辛普森1/3法则和辛普森3/8法则,就是用二次/三次插值曲线来代替直线进行积分。
    在这里插入图片描述

    1.1.2 辛普森1/3法则(3点积分)

    辛普森1/3法则法则积分计算公式:
    在这里插入图片描述

    1.1.3 辛普森3/8法则(4点积分)

    辛普森3/8法则法则积分计算公式:
    在这里插入图片描述

    1.1.4 布尔法则(5点积分)

    布尔法则与辛普森法则类似,只是更高阶的多项式积分,其公式为:
    在这里插入图片描述

    1.1.5牛顿-科茨积分公式总结

    在这里插入图片描述
    注意:1.对于3点积分与4点积分具有相同的误差阶;5点积分与6点积分具有相同的误差阶。对于点数更多的积分也一样。因此优先使用奇数个点的积分。
    2.对于点数大于等于9的牛顿-科茨积分,求积公式的稳定性得不到保证(具体原因可以找数值积分的数看看,本文源自《现代数值积分》第5章 数值积分与数值微分),积分不能保证收敛,因此实际计算中一般不采用高阶牛顿-科茨积分公式。(处于效率考虑,一般也很少用超过5个点的牛顿-科茨积分公式)。

    (穿插)积分改进思路1

    牛顿-科茨积分公式思想:根据积分点构造近似的插值多项式,进而利用多项式积分表示原积分。
    对于高次的多项式插值,会出现龙格现象(对于高次的多项式插值,插值多项式会出现不收敛的现象,成为龙格现象,想深入了解可以去查询相关资料。并且次数越高,越容易出现不收敛现象)。

    针对龙格现象,更好的选择切比雪夫节点来进行插值(想深入了解可以去查询相关资料,本文部分源自《现代数值积分》第3章 多项式插值与样条插值)。

    在这里插入图片描述
    由于切比雪夫节点的特殊性,对于高阶插值的可以收到较好的效果,因此,对于“牛顿-科茨积分公式总结”中提到的点数大于等于9的牛顿-科茨积分,可以采用切比雪夫计算节点,然后构造插值多项式,再计算积分(此方法实现和推广比较麻烦,此处不做扩展,感兴趣的可以自己下去推导实现)。

    1.2 复合积分(复合梯形积分、复合辛普森积分等)

    上述主要介绍的是单个区间的积分。为了计算更准确,我们可以将区间划分为若干个小的区间,然后对每个小区间进行积分,然后对各个小区间的积分求和。由于对于单个小区间的积分可分为梯形法则,辛普森法则,和布尔法则,因此复合积分的每个小区间的积分可以是这些方法中的任意一种。
    下面主要讲一下 复合梯形积分和复合辛普森积分。

    1.2.1 复合梯形积分

    小区间的划分又分为等长区间和不等长区间两种划分方法。
    不等长区间,复合梯形积分计算公式:
    在这里插入图片描述
    区间不等长情况应用场景很有限。

    等长区间,复合梯形积分计算公式:
    在这里插入图片描述
    等长区间,复合梯形积分计算误差:
    在这里插入图片描述

    1.2.2 复合辛普森积分

    同样的,复合辛普森积分区间的划分也可以分为等长区间和不等长区间两种划分方法。
    由于不等长区间应用很少,因此直接介绍等长区间,复合辛普森1/3积分计算公式:

    在这里插入图片描述
    等长区间,复合辛普森1/3积分计算误差:
    在这里插入图片描述

    由于3点积分(辛普森1/3积分)与4点积分(辛普森3/8积分)具有相同的误差阶,因此一般优先使用复合辛普森1/3积分,复合辛普森3/8积分应用极少,因此本文不做介绍。感兴趣的同学可以下去自己推一下复合辛普森3/8积分计算公式,应该也很简单。

    (穿插)积分改进思路2

    可以看出,将区间划分为更小的区间可以有效的提高积分计算的精度。那么我们是不是可以将积分区间不断地细分,从而来得到更精确的积分结果?理论上是可以的!!
    方法1:(区间折半法)
    a.计算2个区间的积分,用复合辛普森1/3法则(当然,也可以用梯形法则或者布尔法则等都可以,但是推荐复合辛普森1/3法则);
    b.将两个区间等分为4个区间,用复合辛普森1/3法则计算四个子区间的积分和;
    c.比较a和b两步计算的结果相对误差是否接近一个很小的数:(Ia - Ib)/Ib<e。如果相对误差满足要求,则结束,否则继续细分区间,用复合辛普森1/3法则计算积分,直到满足要求为止。

    方法2:改进方法1
    	a.对每个子区间一分为二后,比较两个区间的积分值与一个区间的积分值是否相近;
    	b.如果相近则返回该值;否则对每个子区间执行a。(递归实现)
    

    上面说,理论上是可以的!!理论上是可以的!!但是实际不行!!!
    因为计算机浮点数运算是存在误差的。当区间数量很大的时候,浮点数计算的舍入误差(不理解的可以自行去学习了解一下)变得很大,会限制积分的提高。下面就是一个很好的例子。
    在这里插入图片描述

    1.3常用的高精度积分

    上述区间二分的思想是很有用的。在龙贝格积分和自适应积分中都会找到它的影子。

    1.3.1 龙贝格积分

    龙贝格积分的思想,及推导过程:
    在这里插入图片描述
    在这里插入图片描述
    在这里插入图片描述
    在这里插入图片描述
    在这里插入图片描述
    类推:联合两个误差为O(h^4)的积分,改进后可以得到一个误差为O(h ^ 6)的积分:
    在这里插入图片描述
    由此产生龙贝格积分:
    在这里插入图片描述
    在这里插入图片描述

    (穿插)积分改进思路3

    龙贝格积分是用梯形公式推导得到的,那么我们是不是可以用辛普森公式做类似的推导,也能到达一套新的公式呢?答案是:理论上是可以的!!

    由辛普森1/3公式推导出来的加入误差修正的积分计算式:
    在这里插入图片描述
    那么问题就来了。
    我们的计算机表达的double值是有范围的,而对于辛普森外,采用递归,递归8次就是16^8 = 2 ^ 32。计算机递归计算的结果很容易就产生溢出了。除非采用多精度浮点数计算才行,而采用多精度浮点数可能存在一些计算效率的问题。。。更深入的问题我就不继续深入探讨了(有兴趣的可以试试)。。。好了,总之来说,用辛普森推导出来的方法不是很实用。还是梯形龙贝格积分经典好使。

    1.3.2 自适应积分

    龙贝格积分中有两点:
    a。每次都会将所有的步长折半;(那么我们能不能只对计算误差较大的那个区间进行折半成两个子区间再积分呢?)
    b。辛普森公式推导的方法,由于计算机double表示范围有限没用上。

    下面都会在自适应积分中派上用场。自适应积分思想:(递归方法)
    step1.辛普森1/3法则计算区间[a,b]的积分;
    step2.区间步长折半,复合辛普森1/3法则计算区间[a,b]的积分;
    如果两次误差小于给定的误差限,停止执行,返回计算结果
    在这里插入图片描述
    如果两次误差大于给定的误差限,对两个子区间分别执行step1和step2。

    伪代码:
    在这里插入图片描述

    1.3.3 高斯积分

    高斯积分应用比较广泛的是高斯-勒让德公式,此外还有高斯-切比雪夫公式;以及区间【0,正无穷】的广义积分 高斯-拉盖尔公式;以及区间【负无穷,正无穷】的广义积分 高斯-埃尔米特公式。

    高斯-勒让德公式

    介绍高斯积分的文档就很多了,本文不细讲推导的过程,给出计算方法:
    高斯-勒让德
    多点高斯-勒让德公式:
    在这里插入图片描述
    在这里插入图片描述
    高斯勒让德公式的区间为【-1,1】.对于任意区间为【a,b】;需要做一下简单的变量代换,
    在这里插入图片描述

    高斯积分误差:

    高斯勒让德公式的区间为【-1,1】,误差:
    在这里插入图片描述
    对于任意区间为【a,b】,误差:
    在这里插入图片描述
    在这里插入图片描述
    由上式继续往下推导,可以得出,当区间长度【b-a】> 2(n+1)时,误差会发散。因此对于高斯点为n的高斯积分,积分区间长度不要大于2(n+1),否则积分结果很可能会不准确。
    此外,高斯积分的误差还与被积函数的光顺性有关,被积函数光顺性越差,误差越大。详细可参考《现代数值积分》第5章 数值积分与数值微分,P139页)。

    高斯积分的优点:高斯积分时给定节点数下代数精度最高的求积公式。

    高斯积分的不足:高斯积分每次改变积分点的个数,所以的积分点的函数值都需要重新计算。

    高斯-切比雪夫公式

    高斯点为切比雪夫节点。

    (未完待续。。。)

    广义积分:高斯-拉盖尔公式

    广义积分:高斯-埃尔米特公式

    2 多重积分方法

    二重积分的辛普森公式

    二重积分的蒙特卡洛方法

    多重(高维)积分的蒙特卡洛方法

    展开全文
  • 几种数值积分方法

    千次阅读 2016-10-08 11:44:21
    梯形多步法和辛普森积分

    梯形多步法和辛普森积分
    具体算法以及相关定理参考以下链接:

    http://wenku.baidu.com/link?url=ZAScYogajXHrTTRa5xjpUPtS7OQQXZ_LfXaWNkczTtWf2MJgx0RZFUuca4iRGUcPHAVXOyGOIydwz3j-rriBiyLxRMnio7m4fNFcJOKJqui

    #include <stdio.h>

    #include <string.h>
    #include <stdlib.h>
    #include <math.h>

    #define EPS 1.0E-14        //计算精度
    #define DIV 1000        //分割区间,值越大,精确值越高
    #define ITERATE 20        //二分区间迭代次数,区间分割为2^n,初始化应该小一点,否则会溢出

    #define RECTANGLE 0        //矩形近似
    #define TRAPEZOID 1        //梯形近似
    #define TRAPEZOID_FORMULA 2    //递推梯形公式
    #define SIMPSON_FORMULA 3     //辛普森公式
    #define BOOL_FORMULA 4        //布尔公式

    double function1(double);       
    double function2(double);
    double function3(double);
    void Integral(int, double f(double), double, double);      //矩形, 梯形逼近求定积分公式
    void Trapezoid_Formula(double f(double), double a, double b);    //递推梯形公式
    void Simpson_Formula(double f(double), double a, double b);        //递推辛普森公式
    void Bool_Formula(double f(double), double a, double b);        //递推布尔公式
    void Formula(int, double f(double), double, double);

    double function1(double x)
    {
        double y;
        y = cos(x);
        return y;
    }

    double function2(double x)
    {
        double y;
        y=1/x;
       // y = 2+sin(2 * sqrt(x));
        return y;
    }

    double function3(double x)
    {
        double y;
        y = exp(x);
        return y;
    }

    int main()
    {    
        printf("y=cos(x) , x=[0, 1]\n");
        printf("Rectangle : "); Integral(RECTANGLE, function1, 0, 1);
        printf("Trapezoid : "); Integral(TRAPEZOID, function1, 0, 1);
        Formula(TRAPEZOID_FORMULA, function1, 0, 1);
        Formula(SIMPSON_FORMULA, function1, 0, 1);
        Formula(BOOL_FORMULA, function1, 0, 1);
        printf("\n");
        
        printf("y=1/x , x=[1, 5]\n");
        printf("Rectangle : "); Integral(RECTANGLE, function2, 1, 5);
        printf("Trapezoid : "); Integral(TRAPEZOID, function2, 1, 5);
        Formula(TRAPEZOID_FORMULA, function2, 1, 5);
        Formula(SIMPSON_FORMULA, function2, 1, 5);
        Formula(BOOL_FORMULA, function2, 1, 5);
        printf("\n");
        
        printf("y=exp(x) , x=[0, 1]\n");
        printf("Rectangle : "); Integral(RECTANGLE, function3, 0, 1);
        printf("Trapezoid : "); Integral(TRAPEZOID, function3, 0, 1);
        Formula(TRAPEZOID_FORMULA, function3, 0, 1);
        Formula(SIMPSON_FORMULA, function3, 0, 1);
        Formula(BOOL_FORMULA, function3, 0, 1);
        
        return 0;
    }
    //积分:分割-近似-求和-取极限
    //利用黎曼和求解定积分
    void Integral(int method, double f(double),double a,double b)   //f(double),f(x), double a,float b,区间两点
    {
        int i;
        double x;            
        double sum = 0;      //矩形总面积
        double h;            //矩形宽度
        double t = 0;        //单个矩形面积
        
        h = (b-a) / DIV;
         
        for(i=1; i <= DIV; i++)       
        {
            x = a + h * i;     
            switch(method)
            {
                case RECTANGLE :
                        t = f(x) * h;  
                    break;
                case TRAPEZOID :
                        t = (f(x) + f(x - h)) * h / 2;
                    break;
                default:
                    break;
            }
            sum = sum + t;   //各个矩形面积相加
        }
        
        printf("the value of function is %.10f\n", sum);
    }

    //递推梯形公式
    void Trapezoid_Formula(double f(double), double a, double b)//递推梯形公式
    {
        unsigned int i, j, M, J=1;
        double h;
        double F_2k_1 = 0;
        double T[32];
        double res = 0;
        
        T[0] = (f(a) + f(b)) * (b-a)/2;
        for(i=0; i<ITERATE; i++)
        {
            F_2k_1 = 0;
            J = 1;
            
            for(j=0; j<=i; j++)    //区间划分
                J <<=  1;    //2^J
            h = (b - a) /J;    //步长
            //M = J/2;            //2M表示将积分区域划分的块数
            for(j=1; j <= J; j += 2)    //
            {
                F_2k_1 += f(a + j*h);    //f(2k-1)累加和
            }     
            T[i+1] = T[i]/2 + h * F_2k_1;        //递推公式
            res =  T[i+1];
            
            if(fabs(T[i+1] - T[i]) < EPS)
                if(i < 3) continue;
                else break;
        }
        
        printf("Trapezoid_Formula : the value of function is %.10f\n", res);
        //return res;
    }
    //辛普森公式
    void Simpson_Formula(double f(double), double a, double b)        //辛普森公式
    {
        unsigned int i, j, M, J=1;
        double h;
        double F_2k_1 = 0;
        double T[32], S[32];
        double res_T = 0, res_S = 0;
        
        T[0] = (f(a) + f(b)) * (b-a)/2;
        for(i=0; i<ITERATE; i++)
        {
            F_2k_1 = 0;
            J = 1;
            
            for(j=0; j<=i; j++)    //区间划分
                J <<=  1;    //2^J
            h = (b - a) /J;    //步长
            //M = J/2;            //2M表示将积分区域划分的块数
            for(j=1; j <= J; j += 2)    //
            {
                F_2k_1 += f(a + j*h);    //f(2k-1)累加和
            }     
            T[i+1] = T[i]/2 + h * F_2k_1;        //递推梯形公式
            res_T =  T[i+1];
            
            if(fabs(T[i+1] - T[i]) < EPS)
                if(i < 3) continue;
                else break;
        }
        
        for(i=1; i<ITERATE; i++)
        {
            S[i] = (4 * T[i] - T[i-1]) / 3;        //递推辛普森公式
            res_S = S[i];
            if(fabs(S[i] - S[i-1]) < EPS)
                if(i < 3) continue;
                else break;
        }
        
        printf("Simpson_Formula : the value of function is %.10f\n", res_S);
        //return res_S;
    }

    //布尔公式
    void Bool_Formula(double f(double), double a, double b)        //布尔公式
    {
        unsigned int i, j, M, J=1;
        double h;
        double F_2k_1 = 0;
        double T[32], S[32], B[32];
        double res_T = 0, res_S = 0, res_B;
        
        T[0] = (f(a) + f(b)) * (b-a)/2;
        for(i=0; i<ITERATE; i++)
        {
            F_2k_1 = 0;
            J = 1;
            
            for(j=0; j<=i; j++)    //区间划分
                J <<=  1;    //2^J
            h = (b - a) /J;    //步长
            //M = J/2;            //2M表示将积分区域划分的块数
            for(j=1; j <= J; j += 2)    //
            {
                F_2k_1 += f(a + j*h);    //f(2k-1)累加和
            }     
            T[i+1] = T[i]/2 + h * F_2k_1;        //递推公式
            res_T =  T[i+1];
            
            if(fabs(T[i+1] - T[i]) < EPS)
                if(i < 3) continue;
                else break;
        }
        
        for(i=1; i<ITERATE; i++)
        {
            S[i] = (4 * T[i] - T[i-1]) / 3;                //递推辛普森公式
            res_S = S[i];
            if(fabs(S[i] - S[i-1]) < EPS)
                if(i < 3) continue;
                else break;
        }
        
        for(i=1; i<ITERATE; i++)
        {
            B[i] = (16 * S[i] - S[i-1]) / 15;            //递推布尔公式
            res_B = B[i];
            if(fabs(B[i] - B[i-1]) < EPS)
                if(i < 3) continue;
                else break;
        }
        
        printf("Bool_Formula : the value of function is %.10f\n", res_B);
        //return res_B;
    }

    //采用二分区间的方法迭代,实际运行速度很慢
    void Formula(int method, double f(double), double a, double b)//递推梯形公式
    {
        unsigned int i, j, M, J=1;
        double h;
        double F_2k_1 = 0;
        double T[32], S[32], B[32];
        double res_B = 0, res_S = 0, res_T = 0;
        
        T[0] = (f(a) + f(b)) * (b-a)/2;
        for(i=0; i<ITERATE; i++)
        {
            F_2k_1 = 0;
            J = 1;
            
            for(j=0; j<=i; j++)    //区间划分
                J <<=  1;    //2^J
            h = (b - a) /J;    //步长
            //M = J/2;            //2M表示将积分区域划分的块数
            for(j=1; j <= J; j += 2)    //
            {
                F_2k_1 += f(a + j*h);    //f(2k-1)累加和
            }     
            T[i+1] = T[i]/2 + h * F_2k_1;        //递推公式
            res_T =  T[i+1];
            
            if(fabs(T[i+1] - T[i]) < EPS)
                if(i < 3) continue;
                else break;
        }
        
        switch(method)
        {
            default :
            case TRAPEZOID_FORMULA :
                    printf("Trapezoid_Formula : the value of function is %.10f\n", res_T);
                    //return res_T;
                    break;
            case SIMPSON_FORMULA :
                    for(i=1; i<ITERATE; i++)
                        {
                            S[i] = (4 * T[i] - T[i-1]) / 3;
                            res_S = S[i];
                            if(fabs(S[i] - S[i-1]) < EPS)
                                if(i < 3) continue;
                                else break;
                        }
                    printf("Simpson_Formula : the value of function is %.10f\n", res_S);
                    //return res_S;
                    break;
            case BOOL_FORMULA :
                    for(i=1; i<ITERATE; i++)
                    {
                        S[i] = (4 * T[i] - T[i-1]) / 3;
                        res_S = S[i];
                        if(fabs(S[i] - S[i-1]) < EPS)
                            if(i < 3) continue;
                            else break;
                    }
                    for(i=1; i<ITERATE; i++)
                    {
                        B[i] = (16 * S[i] - S[i-1]) / 15;
                        res_B = B[i];
                        if(fabs(B[i] - B[i-1]) < EPS)
                            if(i < 3) continue;
                            else break;
                    }
                    
                    printf("Bool_Formula : the value of function is %.10f\n", res_B);
                    //return res_B;
                break;
        }
        
    }





    展开全文
  • 数值积分,无约束优化问题最速下降法,无约束优化问题的0.618法,共轭梯度法,n_单纯形法,数值法求微分方程,最小二乘算法,多项式乘除,
  • 1. 常用数值积分方法(特别是梯形法、Simpson方法、 Cotes公式、Romberg算法以及Gauss 求积公式)的原理;2. 复合梯形公式,复合Simpson公式,Romberg算法的程序
  • 常用数值计算方法c++源代码实现

    热门讨论 2011-10-11 10:33:28
    二分法求解;.牛顿法求解;高斯消去法求解;雅可比迭代法求解;拉格朗日插值;牛顿插值;最小二乘法拟合;龙贝格方法计算积分;欧拉方法求解初值问题;
  • 数值积分

    2019-06-17 22:47:26
    数值积分法也是计算机仿真模拟中常用的一种方法。在已知函数的微分方程时,求解函数下一时刻的值,我们主要有欧拉法、梯形法和龙格库塔法。 欧拉积分 欧拉积分法是这些方法中精度最低的,但也是最容易编程实现的一...
  • 业 论 文 题 目 数值积分常用算法设计与实现 学 生 * 指导教师 * 教授 年 级 2005?级 专 业 计算机科学与技术 系 别 计算机系 学 院 计算机科学技术与信息工程学院 * 2019?年?06?月 论文提要 本文应用插值积分法和...
  • 积分基本概念2.欧拉积分3.中值积分4.RK4积分(4阶龙格库塔法) 1.积分基本概念 非线性微分方程: 在有限的时间间隔Δt积分: 连续时间内: 2.欧拉积分 欧拉方法假设导数f(·)在区间内是恒定的,有...
  • MATLAB进行数值积分的主要函数: 1.trapz 梯形法求解积分 x=0:pi/10:pi; y=sin(x); trapz(x,y) 2.quad 基于变步长simpso法求积分 q = quad(fun,a,b,tol) 其中fun是被积函数文件名或函数句柄,a, b是积分下限和...
  • 1. 积分基本概念 设F(x)为函数f(x)的一个原函数,我们把函数f(x)的所有原函数F(x)+C(C为任意常数)叫做函数f(x)的不定积分(indefiniteintegral)...他们的区别就是如何用数值方法模拟一个斜率。这里简单总结如下: ...
  • 常微分方程的解法 (一): 常微分方程的离散化 :差商近似导数、数值积分方法、Taylor 多项式近似 常微分方程的解法 (二): 欧拉(Euler)方法 常微分方程的解法 (三): 龙格—库塔(Runge—Kutta)方法 、线性多步法 ...
  • 高斯勒让德数值积分公式

    千次阅读 2019-12-15 23:06:00
    高斯勒让德数值积分公式1 引言2 高斯积分公式2.1 一维区间上高斯数值积分公式2.2 二维三角形上的高斯数值积分...在插值型积分公式中,高斯积分公式具有最高的代数精度,所以是常用数值积分公式。 2 高斯积分公式...
  • MATLAB数值积分

    2019-09-25 07:01:24
    许多工程技术和数学研究中要用到定积分,如果无法直接算不出精确值(如含在积分方程中的积分)或计算困难但可用近似值近似时,就用数值积分方法加以解决。常用的算法有:复化梯形、辛甫生(Simpson)、柯特斯...
  • Matlab数值微分与数值积分

    千次阅读 2020-06-11 22:40:13
    matlab中数值微分与数值积分的实现
  • 数值积分之初步介绍

    千次阅读 2018-07-04 22:13:25
    为什么要数值积分呢?也就是数值积分的必要性问题?数值积分的必要性源自计算函数的原函数的困难性。利用原函数计算定积分的方法建立在牛顿—莱布尼兹之上。然而,原函数可以用初等函数表示的函数为数不多,大部分的...
  • 通过MATLAB实现微积分数值解法
  • Java 常用数值算法集

    2009-04-17 20:43:57
    本书共有数值计算中常用的Java方法近200个,内容包括:解线性代数方程组、插值、数值积分、特殊函数、函数逼近、随机数、排序、特征值问题、数据拟合、方程求根和非线性方程组求解、函数的极值和最优化、傅里叶变换谱...
  • Visual C++ 常用数值算法集

    热门讨论 2012-03-19 11:57:59
    Visual C++ 常用数值算法集 作者:何光渝编 出版社:科学出版社 出版日期:2002年7月 ISBN:703010498 序 前言 第1章 线性代数方程组的解法 1.1全主元高斯-约当(Gauss-Jordan)消去法 1.2LU分解法 1.3追赶法 ...
  • 数值积分法原理及matlab程序实现

    万次阅读 2016-04-16 10:14:48
    数值积分法是求定积分的近似值的数值方法。即用被积函数的有限个抽样值的离散或加权平均近似值代替定积分的值。  求某函数的定积分时,在多数情况下,被积函数的原函数很难用初等函数表达出来,另外,许多实际问题...
  • 在科学研究和工程设计中经常需要做大量的数值...主要内容有:线性代数方程组的数值解法,非线性方程和方程组的迭代解法,矩阵特征值和特征向量的计算,函数的插值与逼近,数值积分,常微分方程和偏微分方程的数值解法.
  • MATLAB学习笔记:数值积分

    万次阅读 多人点赞 2018-01-11 20:37:41
    当 (1)被积函数的原函数不能用...就应该建立定积分的近似计算方法:数值积分方法。 梯形法: z=trapz(x,y) >> x=0:0.5:1; >> y=exp(-x.^2); >> z=trapz(x,y)z = 0.7314 >> x=0:0.05:1;
  • 布料仿真中常用积分方法

    千次阅读 2016-11-28 22:36:41
    1. 简介布料仿真中,我们通常将布料剖分为三角形网格(或四边形网格),并用弹簧-质点模型构造动力学系统:质点即三角形的顶点...常用方法有:显式/隐式欧拉法,Symplectic Euler,Midpoint method,Leapfrog integrati
  • Visual C++常用数值算法集 本收共不数值计算中常用的Visual C++子过程近200个, 内容包括:解线性代数议程组、插值、数值积分、特殊函数、函数逼近、随机数、排序、 特征值问题、数据拟合、方程求根和非线性方程组...
  • Visual FORTRAN 常用数值算法集是一本十分经典的Fortran数值算法介绍书籍。在介绍Fortran基本语法的基础上,集中讲述利用Fortran语言实现常用的数值算法,包括:插值、积分、线性方程组解法、逼近、拟合、特征值计算...
  • 内容包括:解线性代数方程组、插值、数值积分、特殊函数、函数逼近、随机数、排序、特征值问题、数据拟合、方程求根和非线性方程组求解、函数的极值和最优化、傅里叶变换谱方法、数据的统计描述、解常微分方程组、两...

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 12,252
精华内容 4,900
关键字:

常用的数值积分方法