常用的算法有,怎么着求n次多项式的积分

笔者:凯鲁嘎吉 – 网易
http://www.cnblogs.com/kailugaji/

 

  用程序来求积分的办法有广大,这篇作品紧如若关于牛顿(牛顿)-科特斯公式。

一、实验目标

机械求积法

  学过插值算法的同学最容易想到的就是用插值函数代替被积分函数来求积分,但其实在多数情景下这是无用的。

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

转载请讲明出处!

  插值函数一般是一个不领先n次的多项式,假如用插值函数来求积分的话,就会推荐高次多项式求积分的题材。这样会将本来的求积分问题带到另一个求积分问题:怎么着求n次多项式的积分,而且当次数变高时,会并发龙悲歌现象,误差反而可能会附加,并且高次的插值求积公式有可能会变得不安宁:详细原因不赘述。

二、实验原理

一、引言

  随着人工智能的勃兴,在处理器世界又一回吸引了数学热,不管是观念的机器学习,仍旧前些天的吃水学习,都离不开积分的扶助,这总结机在底部到底是何等求积分的啊?小编同我们一齐探索。

  Newton-科特斯公式解决这一问题的措施是将大的插值区间分为一堆小的插值区间,使得多项式的次数不会太高。然后经过引入参数函数

图片 1

二、理论推导

    大家了然,在我们所学的微积分中我们是经过牛顿-莱布尼兹公式举办求解,可是在其实使用中咱们反复会遭遇相比较复杂的函数,他们的原函数我们往往是找不到的,这么些时候咱们相应怎么求解呢?

  我们不难想的点子是定义法,也就是把区间举行剪切,当分点相当多的时候我们就可以用矩形面积代替曲线所围成的面积,不过我们为了拿到精度很高的结果往往需要划分等多区间,那样统计的次数将大大扩大。

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

  机械求积分法前戏:

    
在微积分中我们求定积分时不仅有牛顿(牛顿)-莱布尼兹公式,同时还有积分中值定理:

     若函数图片 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)))。

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

 

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.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:使用辛普森求积,对于低精度或者不光滑函数成效更高

    图片 22

quadl:该函数使用了名为洛巴托求积(Lobatto
Quadrature)的算法,对于高精度和光滑函数效用更高

  分为20个小区间时:

应用方法:

  图片 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位有效数字

相关文章