就能推荐高次多项式求积分的难点,常用的算法有

小编:凯鲁嘎吉 – 天涯论坛
http://www.cnblogs.com/kailugaji/

 

  用程序来求积分的格局有为数非常多,那篇小说首即使有关Newton-科特斯公式。

一、实验指标

机械求积法

  学过插值算法的校友最轻巧想到的正是用插值函数代替被积分函数来求积分,但骨子里在大比很多场合下那是低效的。

洋洋工程技艺和数学切磋中要用到定积分,要是不能够直接算不出准确值(如含在积分方程中的积分)或总括困难但可用近似值近似时,就用数值积分法方法加以化解。常用的算法有:复化梯形、辛甫生(Simpson)、柯特斯(Cotes)求积法; 龙贝格(Romberg)算法;高斯(Gauss)算法。

转发请注解出处!

  插值函数一般是五个不当先n次的多项式,假若用插值函数来求积分的话,就能够推荐高次多项式求积分的难题。那样会将原本的求积分难点带到另贰个求积分难点:怎么样求n次多项式的积分,何况当次数变高时,会冒出龙悲歌现象,零值误差反而恐怕会附加,并且高次的插值求积公式有望会变得不牢固:详细原因不赘述。

二、实验原理

一、引言

  随着人工智能的兴起,在微型Computer领域又二遍吸引了数学热,不管是守旧的机械学习,依旧明天的纵深学习,都离不开积分的帮忙,那Computer在尾巴部分到底是怎么求积分的呢?作者同大家齐声研商。

  Newton-科特斯公式消除这一题指标方式是将大的插值区间分为一批小的插值区间,使得多项式的次数不会太高。然后经过引进参数函数

图片 1

二、理论推导

    大家明白,在大家所学的微积分中我们是经过牛顿-莱布尼兹公式实行求解,但是在其实使用中大家反复会蒙受比较复杂的函数,他们的原函数大家再三是找不到的,那个时候大家理应怎么求解呢?

  大家不难想的措施是定义法,也等于把区间举办剪切,当分点非常多的时候我们就可以用矩形面积代替曲线所围成的面积,然则大家为了获得精度非常高的结果往往必要划分等多区间,那样测算的次数将大大扩大。

  那应该怎么优化呢?这里大家介绍一种求积分的议程:机械求积法。

  机械求积分法前戏:

    
在微积分中我们求定积分时不唯有有Newton-莱布尼兹公式,同有的时候候还应该有积分中值定理:

     若函数图片 2闭区间图片 3
上三回九转,,则在积分区间图片 4上至少存在三个点图片 5,使下式成立

    图片 6

    其中,a、b、图片 7满足:a≤图片 8≤b。\[

    我们从上式子出发,将求积分的标题转化为找某四个图片 9的标题,那怎么替代呢?

    即便大家的函数是二个0次函数也正是一个F(x)
= C
是贰个常函数时候,在[a,b]间隔的定积分正是(a-b)*f(a),也正是矩形面积,那样我们的f(图片 10)=f(a),我们后续加盟是叁次的啊?F(x)
= ax + b

    大家很轻易想到这几个图形是二个梯形。

                                              图片 11

    因此其积分=(b-a)*[f(a)+f(b)]/2,那就是我们的梯形公式,这里的f(图片 12)=[f(a)+f(b)]/2,那难点来了,假诺函数是一个畸形的曲线呢?那大家该怎么办吗?

 

图片 13

    很轻松想到的就是把区间细分成相当多少个小区间,然后在各样区间找叁个老少咸宜的f(图片 14),然后再求积分,再求和就好了,是的那正是们的思路,那么些艺术正是机械求积法。大家上边给出定义:

    机械求积法:

      

                      图片 15

图片 16

三、实验程序

三、代数精度概念

  
在知晓机械求积公式之后,这大家怎么核准贰个求积公式的好坏的?这里我们引入代数精度的定义。

  定义:

     如若某些求积公式对于次数不当先m的多项式均能纯粹创设,而对此m+1次的的多项式不标准成立,则称该公式具备m次代数精度。

     (代数精度越高,越标准)

  我们来探视梯形公式的代数精度:(在证实时候倘诺取1,x,x^2,x^3…….等就行了,其余都足以由那一个组合而成)

    F(x) = (b-a)*[f(a)+f(b)]/2,当 f(x) =
1时候,设a=0,b=1,明显创设,而f(x) = x
时,用Newton-莱布尼兹公式算得:x^2/2|10
= 0.5,而用梯形公式获得的也是0.5,那样对于f(x) = x时也规范创制,而当f(x) =
x^2时候则不成了,那么大家就说

  梯形公式具备叁遍代数精度,其余的代数精度的规定方法也同上。

 

将带有幂的项的取值范围固定在一个恒定范围内,那样一来就将多项式带有幂的一部分的求积变为叁个一定的常数,只需手工业算出来就能够。那么些常数能够一直带入多项式求积函数。

上面给出复化Simpson求积法程序(梯形及柯特斯复化求积分程序可根据编制):

四、插值型求积公式

  大家在其次有个别时候知道,大家的目标就是每每的分开区间直到有主意准确的求出积分(找到f(图片 17)),其余呢,大家还是能够透过插值法拟合曲线,把难点转化为数值难题。

  (未有接触插值法的能够活动)见鬼吧拉格朗日插值法

  这里我们用插值获得的插值多项式子代替原函

    图片 18

   

 图片 19

  

那正是插值型积分公式,那样我们就因此插值法能够减流年算,然则那并非大家最终要的,那只是个起来,更牛的还在后面,请持续关注杰克总结方式层层博客。

转发请注脚出处!

  上式中x的求积分区间为[a, b],h = (b – a)/n,
那样一来积分区间变为[0,
n],必要留神的是从这么些公式能够阅览叁个大的间距被分成n个等长的小区间。
这一部分具体请参见任意一本有关数值计算的书!

图片 20

 点击看看主页别的剧情呗! 

 

     

   n是一个刚开始阶段明确好的值。

四、实验内容

  又因为一个大的插值区间须要被分为等长的四个小区间,并在那几个小区间上分别张开插值和积分,由此此时的Newton-科特斯公式被叫做:复化Newton-科特斯公式。

选一可准确算值的定积分,用复化的梯形法及复化Simpson求积法作近似总结,并比较结实。

   况且对于n的例外取值Newton-科特斯有两样的称呼:
当n=1时,叫做复化梯形公式,复化梯形公式也便是将每叁个小区间都看为贰个梯形(高为h,上底为f(t),
下底为f(t+1))。那与积分的本色:最为分隔  相同。

 

  当n=2时,复化Newton-科特斯公式被叫作复化Simpson公式(非United States法律界盛名的特别Simpson)。

五、解答(按如下顺序提交电子版)

 

1.(程序)

  小编那篇小说完毕的是复化梯形公式:

xps.m:

    图片 21

function y=xps(x)
y=x^(3/2);

 

复化梯形公式:

  首先写八个函数求节点函数值求和那部分:

   trap.m:

    

function [T,Y,esp]=trap(a,b,n)
h=(b-a)/n;
T=0;
for i=1:(n-1)
    x=a+h*i;
    T=T+xps(x);
end
T=h*(xps(a)+xps(b))/2+h*T;
syms x
Y=vpa(int(xps(x),x,a,b),8);
esp=abs(Y-T);
"""
@brief: 求和 ∑f(xk) : xk表示等距节点的第k个节点,不包括端点
        xk = a + kh (k = 0, 1, 2, ...)
        积分区间为[a, b]

@param: xk      积分区间的等分点x坐标集合(不包括端点)
@param: func    求积函数
@return: 返回值为集合的和
"""
def sum_fun_xk(xk, func):
    return sum([func(each) for each in xk])

复化辛甫生(Simpson)公式:

  

   simpson.m:

  然后就能够写整个求积分函数了:

function [SI,Y,esp]=simpson(a,b,m)
%a,b为区间左右端点,xps(x)为求积公式,m*2等分区间长度
h=(b-a)/(2*m);
SI0=xps(a)+xps(b);
SI1=0;
SI2=0;
for i=1:((2*m)-1)
    x=a+i*h;
    if mod(i,2)==0
        SI2=SI2+xps(x);
    else 
        SI1=SI1+xps(x);
    end
end
SI=h*(SI0+4*SI1+2*SI2)/3;
syms x
Y=vpa(int(xps(x),x,a,b),8);
esp=abs(Y-SI);
"""
@brief: 求func积分 :

@param: a  积分区间左端点
@param: b  积分区间右端点
@param: n  积分分为n等份(复化梯形求积分要求)
@param: func  求积函数
@return: 积分值
"""   
def integral(a, b, n, func):
    h = (b - a)/float(n)
    xk = [a + i*h for i in range(1, n)]
    return h/2 * (func(a) + 2 * sum_fun_xk(xk, func) + func(b))

2.(运算结果)

  非常的差不离

>> [T,Y,esp]=trap(1,2,8)

T =

    1.8636


Y =

1.8627417


esp =

0.0008089288247354886607354274019599
>> [SI,Y,esp]=simpson(1,2,8)

SI =

    1.8627


Y =

1.8627417


esp =

0.000000020499792974248975951923057436943

 

从总括结果看:复化Simpson公式更确切。

  

3.(拓展(方法创新、体会等))

  试验:

MATLAB中有局地停放函数,用于施行自适应求积分,都以基于Gander和Gautschi构造的算法编写的。

  当把大间隔分为五个小区间时:

quad:使用Simpson求积,对于低精度恐怕不光滑函数功效更加高

    图片 22

quadl:该函数使用了可以称作洛巴托求积(Lobatto
Quadrature)的算法,对于高精度和光滑函数效用越来越高

  分为二十一个小区间时:

使用方法:

  图片 23

I=quad(func,a,b,tol);

  求的积分值正是那几个多彩的梯形面积之和。

func是被积函数,a,b是积分限,tol是可望的相对抽样误差(若是不提供,默以为1e-6)

 

例如说对于函数f=xe^x在[0,3]上求积分,分明能够通过解析解知道结果是2e^3+1=41.171073846375336

  测量试验代码:

先成立一个M文件xex.m

if __name__ == "__main__":

    func = lambda x: x**2
    a, b = 2, 8
    n = 20
    print integral(a, b, n, func)

    ''' 画图 '''
    import matplotlib.pyplot as plt
    plt.figure("play")
    ax1 = plt.subplot(111)
    plt.sca(ax1)

    tmpx = [2 + float(8-2) /50 * each for each in range(50+1)] 
    plt.plot(tmpx, [func(each) for each in tmpx], linestyle = '-', color='black')

    for rang in range(n):
        tmpx = [a + float(8-2)/n * rang, a + float(8-2)/n * rang, a + float(8-2)/n * (rang+1), a + float(8-2)/n * (rang+1)]
        tmpy = [0, func(tmpx[1]), func(tmpx[2]), 0] 
        c = ['r', 'y', 'b', 'g']
        plt.fill(tmpx, tmpy, color=c[rang%4])
    plt.grid(True)
    plt.show()

内容如下:

  注意上边代码中的n并不是上文开篇提到的公式中的n,开篇提到的n是指将每叁个切实的插值区间(也正是小间隔)等距插n个节点,复化梯形公式的n是一定的为1.

function f=xex(x)

  而代码中的n指将大间隔分为n个小间隔。

f=x.*exp(x);

 

下一场调用:

>> format long

>> format compact

>> quad(@xex,0,3)

ans =

  41.171073850902332

足见有9位有效数字,精度依然蛮高的。

借使利用quadl:

>> quadl(@xex,0,3)

ans =

  41.171074668001779

相反只有7位有效数字

相关文章