复化求积公式的算法及其应用

摘要

在数值计算中，低阶牛顿柯特斯求积方法存在很多缺陷，从余项公式可以看出其要求提高求积公式的代数精度，必须增加结点个数，会导致插值多项式出现龙格现象，且数值稳定性不能保证.基于以上原因，我们往往采用复化求积方法，此方法不仅可以克服以上缺点而且便于在计算机上实现，值得研究和学习.

在本课程设计中，我们首先从复化求积公式的思想引入，然后详细介绍复化梯形求积公式、复化辛普森求积公式和复化柯特斯求积公式的推导过程和相关性质，再对三种求积公式进行比较和总结，其次画出三种求积公式的流程图，最后通过求解例题写出三种求积算法的程序设计.

关键词复化求积算法;流程图;程序设计

目录

引言 (1)

第一章复化求积算法 (2)

§1.1复化求积公式 (2)

§1.1复化求积公式的思想 (3)

§1.2复化求积公式的构造 (3)

§1.2复化梯形求积公式 (3)

§1.2.1复化梯形求积公式的推导过程 (3)

§1.2.2复化梯形求积公式的性质 (3)

§1.3复化辛普森求积公式 (4)

§1.3.1复化辛普森求积公式的推导过程 (4)

§1.3.2复化辛普森求积公式的性质 (4)

§1.4复化柯特斯求积公式 (5)

§1.4.1复化柯特斯求积公式的推导过程 (5)

§1.4.2复化柯特斯求积公式的性质 (5)

§1.5三种复化求积公式的比较及总结 (6)

第二章复化求积公式算法的流程图及其应用 (9)

§2.1 流程图 (9)

§2.2 应用 (12)

参考文献 (15)

附录A (16)

附录B (17)

附录C (18)

引言

积分计算在分析数学领域里是个古老的问题，在数值分析中已被广泛应用.但在计算机上却不能像在分析数学中那样，用原函数[满足)()('x f x F =的函数)(x F 就是函数)(x f 的原函数]计算积分.这是因为在实际问题中，函数关系往往是用列表数据或曲线给出的.即使知道了函数的表达式，求其一个原函数并非一个简单问题.许多函数难以用初等函数表示（如2

,/sin x e x x -等）.在计算机上，通常利用函数的若干个离散值，以代数运算近似计算积分值，这类近似计算法称为数值积分法.

设给定区间],[b a 上的函数)(x f .需要建立计算积分dx x f f I b

a ?=)()(的近似方法.

数值积分的基本思想是试图用一个简单又易于积分的函数逼近)(x f ，以计算积分

)(f I .显然插值多项式是一个很好的选择，因为插值多项式可由)(x f 的若干值构造出

来，其积分很容易计算.为此，需将],[b a 分为n 等分n i x x i i ,,2,1],,[1 =+，其中

b x x x x a n =<<<<=+1321 .分割步长h ，因此,1,3,2,/)1(1+=-+=n i h i x x i 对应的函数值)()(,),(),()(121b f x f x f x f a f n ==+ .显然)(f I 可以表示为所有小区间上各函数的积分的和，即

)()(1f I f I n

i i ∑==

其中 dx x f I i i

x x i ?

+=1

)(

通常把为每个)(f I i 建立的计算公式简称为求积公式，而把)(f I 建立的求积公式

称为复化求积公式.

由于在实际计算时，不宜使用高阶的牛顿——柯特斯公式，但若积分区间较大，单独用一个低阶的牛顿——柯特斯公式来计算积分的近似值，显然精度不好，为了提高数值求积的精确度，可利用积分对区间的可加性来解决这个问题，这就是通常采用

的复合求积法.而且使用这种方法之后，求积公式的收敛性和稳定性也得到了改善.

第一章 复化求积算法

牛顿—柯特斯公式的求积余项表明，求积节点n 越大，对应的求积公式精度越高，但由于牛顿—柯特斯公式在8>n 时数值不稳定，因此不能用增加求积节点数的方法来提高计算精度.实用中常将求积区间],[b a 分成若干个小区间，然后在每个小区间上采用数值稳定的牛顿—柯特斯公式求小区间上的定积分，最后把所有小区间上的计算结果相加来作为原定积分的近似值.采用这种方法构造的求积公式就称为复合求积公式.复合求积公式具有计算简单且可以任意逼近所求定积分值的特点，这是牛顿—柯特斯公式一般做不到的.常用的复合求积公式有复合梯形求积公式和复合辛普森求积公式以及复合柯特斯求积公式.以下我们将从三种复化求积算法的构造、余项、稳定性、收敛性等几方面进行讨论，并写出相应的流程图以及应用中所涉及到的算法的程序设计.

§1.1复化求积公式

§1.1.1 复化求积公式的思想

n 很大时，牛顿——柯特斯求积公式出现了不稳定、不收敛现象，往往使用低阶牛顿——柯特斯求积公式，误差比较大，故将],[b a 若干等分，在每个子区间上反复使用低阶牛顿——柯特斯公式，进行累加.而构造出来的新的求积公式，称之为复化求积公式.在构造求积公式的过程中，我们将求积区间],[b a 进行等距细分：

n i n

a

b i

a x i ,,1,0, =-+=，在每个小区间],[1i i x x -上用相同的“基本”求积公式(如梯形公式；中矩形公式；左(右)矩形公式或辛普森公式)计算出dx x f i i x x ?

-1

)(的近似值i S .

§1.1.2 复化求积公式的的构造

将定积分?b

a dx x f )(的区间],[

b a 划分为n 等分，各节点为kh a x k +=，

n k ,,1,0 =，n

a

b h -=

，在子区间)1,,1,0](,[1-=+n k x x k k 上使用牛顿——柯特公式，将],[1+k k x x 分割为l 等份，步长为l h

，节点为

1,,2,,+=+++k k k k k x l

lh

x l h x l h x x

记

121,,,,++

+

+

=k l

l k l

k l

k k x x

x

x

x

为在],[1+k k x x 上作)(x f 的l 阶牛顿——柯特斯求积公式.

∑∑?

=+

+=+=-=≈+l

i l

i k l i l

i k l

i l i k k k i x x x

f C h x

f C x x I dx x f k

)(0

)(1)

()

()()()(1

由积分区间的可加性，可得

n

l

i k n k l

i l i n k k l n k k k b

a

I x

f C h I dx

x f dx x f ==≈=+-==-=-=+∑∑∑∑?

?

)()()(1

00

)(1

)(1

1

§1.2 复化梯形求积公式

§1.2.1 复化梯形求积公式的的推导过程

将积分区间],[b a 划分等分，步长n

a

b h -=

，求积节点kh a x k +=，n k ,,1,0 =在每个小区间)1,,1,0](,[1-=+n k x x k k 上应用梯形公式

)]()([2

)(11

++≈

?+k k x x x f x f h

dx x f k k

然后将它们累加求和，作为所求积分I 的近似值.

])()(2)([2)]())()()((2)([2

)]

()([2

)()(1

1

121011

01

1

∑∑∑?

?---+-=-=++=+++++=+≈==+n i k n n k k n k n k x x b

a

b f x f a f h

x f x f x f x f x f h

x f x f h

dx x f dx x f I k k

记

n T )]()(2)([21

1

b f x f a f h

n i k ++=∑-=

式为复化梯形求积公式，下标n 表示将区间n 等分，若把区间n 2等分，在每个小区间上仍用梯形求积公式，则可得到n n T T ,2和n H 间的关系为：

)(2

1

2n n n H T T +=

其中

∑=--+=n

k n n

a

b k a f h H 1]2)

12([ §1.2.2复化梯形求积公式的性质

性质1.1复化梯形求积公式余项

当)(x f 在],[b a 上有连续的二阶导数，则复化梯形公式的余项：

)(12

)()('

'2ηf h a b T dx x f R n b

a T --

=-=? ],[b a ∈η 性质1.2稳定性

若],[,)(''b a x M x f ∈≤，则有估计式

M n

a b R n

T 2

312)(-≤ 复化梯形求积公式的系数均大于零，且满足

a b nh n h

A n

i i -==+-+=∑

=]1)1(21[2

因此，复化梯形求积公式的计算过程是数值稳定的.

性质1.3收敛性

可证复化梯形求积公式是收敛的. 性质1.4代数精度

定义1.1 若积分?b a

dx x f )(的数值积分公式?b

a

dx x f )()(0

k n

k k x f A ∑=≈对于任意一个

次数不高于m 次的多项式都精确成立，且存在一个1+m 次多项式使之不精确成立，则称该数值积分公式的代数精度为m .

可证复化梯形求积公式的代数精度为2.

§1.3 复化辛普森求积公式

§1.3.1 复化辛普森求积公式的的推导过程

将积分区间],[b a 划分等分，记子区间],[1+k k x x 的中点为h x x k k 2

1

2

1+=+

在每个

小区间上应用辛普森公式，则有

))

()(2)(4)((6)

444(6)]

()(4)([6)()(101

12

12

11223112101211

01

1

b f x f x f a f h

f f f f f f f f f h

x f x f x f h

dx

x f dx x f I n k n k k k n n n k k k n k n k x x b

a

k k

+++=+++++++++=

++≈==∑∑∑∑?

?-=-=+--++

-=-=+

其中

h x x

k k 2

12

1

+=+

记 )]()(2)(4)([61

1102

1b f x f x f a f h

S n k k n k k n +++=∑∑-=-=+

式为复化辛普森求积公式

§1.3.2复化辛普森求积公式的性质

性质1.5复化辛普森求积公式余项

当)(x f 在],[b a 上有连续的四阶导数，复化辛普森公式的求积余项为：

)(2880

)()2(180)

4(4)4(4ηηf h a b f h a b R S --=--

= ],[b a ∈η 性质1.6稳定性

同复化梯形求积公式，复化辛普森求积公式的系数均大于零，且满足总和为a b - 因此，复化辛普森求积公式的计算过程是数值稳定的.

性质1.7收敛性

可证复化辛普森求积公式是收敛的. 性质1.8代数精度

可证复化辛普森求积公式的代数精度为4.

§1.4 复化柯特斯求积公式

§1.4.1 复化柯特斯求积公式的的推导过程

将积分区间],[b a 划分等分，若把每个子区间],[1+k k x x 四等份，内点依次记为

4

32

14

1,,+

+

+

k k k x

x

x

，同理可得复化柯特斯求积公式

)](7)(14)(32)(12)(32)(7[9010101

14

3211041b f x f x f x f x f a f h

C n k n k n k k k k n k k n +++++=∑∑∑∑-=-=-=++-=+

(1-1)

其中

h x x

h x x h x x k k k k k k 4

3;21

;41

4

32

14

1+

=+=+=+

++ 记(1-1)为复化柯特斯求积公式

§1.4.2复化柯特斯求积公式的性质

性质1.9复化柯特斯求积公式余项

当)(x f 在],[b a 上有连续的四阶导数，复化柯特斯公式的求积余项为：

)()4

(945)(2)

6(6

ηf h a b R c --

= ],[b a ∈η

性质1.10稳定性

同复化梯形求积公式，复化柯特斯求积公式的系数均大于零，且满足总和为a b - 因此，复化柯特斯求积公式的计算过程是数值稳定的.

性质1.11收敛性

可证复化柯特斯求积公式是收敛的. 性质1.12代数精度

可证复化柯特斯求积公式的代数精度为6.

§1.5 三种复化求积公式的比较及总结

为了更形象的表述三种复化求积公式之间的关系，我们通过一个例子来进行比较

例1.1使用各种复化求积公式计算定积分dx x

x

I ?=10sin 为简单起见，依次使用8阶复化梯形公式、4阶复化辛普森公式和2阶复化柯特斯公式,可得各节点的值如下表

表1-1节点值

94569086

.0)]1()(2)0([161

7

1

8=++=∑=f x f f T k k 94608331

.0)]1()(2)(4)0([241

30312

14=+++=∑∑==+f x f x f f S k k k k 94608307

.0)]1(7)(14)](32)(12)(32[)0(7[18011

1104

342412=+++++=∑∑==+++f x f x f x f x f f C k k k k k k 比较三个公式的结果：

精度最低 94569086.08=T 精度次高 94608331.04=S 精度最高 94608307.02=C

原积分的精确值为

6719460830703.0sin 1

0==?

dx x

x

I . 我们知道，三种求积公式的余项分别如表1-2

表1-2 复化梯形、辛普森、柯特斯求积公式的余项

定义1.2对于复化求积公式n I 若存在0>p 及0≠c ，使其余项n I I -满足

c h I I p

n

h =-→0

lim

则称复化求积公式n I 是p 阶收敛的 P 阶收敛性的意义：

对于一个数值求积公式来说，收敛阶越高，近似值n I 收敛到真值dx x f b

a ?)(的速度

就越快.

由于三种求积公式的余项分别是h 的2,4,6阶无穷小量 所以n n n C S T ,,趋于定积分I 的速度依次更快.

从这三种求积公式的构造过程中可以看出，它们都属于机械求积公式，但不属于插值行和牛顿柯特斯公式.都具有稳定性和收敛性，且收敛速度一个比一个快，一个比一准确.在使用函数值个数相等的情况下，248,,C S T 的精度逐渐升高.

第二章 复化求积公式算法的流程图及其应用

§2.1 流程图

1. 复化梯形求积公式

图2.1 复化梯形求积公式算法的流程图

Step1给出被积函数)(x f 、区间],[b a 端点b a ,和等分数n ； Step2求出,kh x k =n

a

b h -=

； Step3计算∑-=1

)(),(),(n k k x f b f a f ；

Step4得)]()()([21

1

b f x f a f h T n k k n ++=∑-=

2. 复化辛普森求积公式

图2.2 复化辛普森求积公式算法的流程图

Step1 给出被积函数)(x f 、区间],[b a 端点b a ,和等分数n ； Step2求出,kh x k =n

a

b h -=

； Step3计算∑∑-=+-=1

2

11

)(,)(),(),(n k k n k k x

f x f b f a f ；

Step4得)]()(2)(4)([61

1102

1b f x f x f a f h

S n k k n k k n +++=∑∑-=-=+

3. 复化柯特斯求积公式

图2.3 复化柯特斯求积公式算法的流程图

Step1给出被积函数)(x f 、区间],[b a 端点b a ,和等分数n ；

Step2求出,kh x k =n

a

b h -=

； Step3计算∑∑∑∑-=-=+

-=+

-=+1

1

1

4

31

2

11

4

1)(,)(,)(,)(),(),(n k k n k k n k k n k k x f x

f x

f x

f b f a f ；

Step4

得)](7)(14)(32)(12)(32)(7[9010101

14

3211041b f x f x f x f x f a f h

C n k n k n k k k k n k k n +++++=∑∑∑∑-=-=-=++-=+

§2.2 应用

例2.1.分别用复化梯形，复化辛普森，复化柯特斯公式计算函数32)(x x x f -=在区间]1,0[上的弧长S .(要求写出源程序和运行结果) *注 在],[b a 上的弧长

dx x f S b

a

?

+=2'))((1

1.用复化梯形公式计算S 的过程：

(1).写出变量说明

表2-1 复化梯形求积公式程序设计的变量说明

Step1 输入n ，n

a

b h -=

，被积函数0),(1=s x f ； Step2 for 1=k to 1-n ；

{计算11)(s kh a f s →++} ))(2)((2

1b f s a f h

s ++=

；

Step3 输出近似值s .

(3) 写出源程序和运行结果(见附录A) 2.用复化辛普森公式计算S 的过程: (1).写出变量说明

表2-2 复化辛普森求积公式程序设计的变量说明

Step1:输入n ，n

a

b h -=

，被积函数0),(1=s x f 0,2=s ； Step2:for 1=i to 1-n ，2+=i i ；

{计算11)2/*(s h i a f s →++} Step3:for 2=j to 1-n ，2+=j j ； {计算22)2/*(s h j a f s →++}

))(24)((6

21b f s s a f h

s +++=

； Step4:输出近似值s .

(3).写出源程序和运行结果(见附录B) 3.用复化柯特斯公式计算S 的过程: (1).写出变量说明

表2-3 复化柯特斯求积公式程序设计的变量说明

Step1输入n ，n

a

b h -=

，被积函数0),(1=s x f 0,2=s 0,3=s ； Step2 for 1=i to 1-n ，2+=i i ；

{计算11)4/*(s h i a f s →++} Step3:for 2=j to 1-n ，4+=j j ； {计算22)4/*(s h j a f s →++} Step4: for 4=k to 2-n ，2+=k k ； {计算33)4/*(s h k a f s →++}

))(141232)((90

321b f s s s a f h

s ++++=

； Step5:输出近似值s .

(3).写出源程序和运行结果(见附录C)

根据运行结果可知，由三种复化求积公式求得的S 的值分别为064837.1、061199.1、061189.1，精度逐渐升高.

参考文献

[1] 薛毅，耿美英.数值分析[M]. 北京:北京工业大学出版社.2003年.

[2] 刘长安.数值分析教程[M].西安:西北工业大学出版社.2005年.

[3] 朝伦巴根，贾德彬.数值计算方法[M].北京:中国水利水电出版社.2007年.

[4] 韩旭里，万中.数值分析与实验[M].北京: 科学出版社.2006年.

[5] 林成森.数值分析[M].北京: 科学出版社.2007年.

[6] 封建湖，车刚明，聂玉峰.数值分析原理. 北京: 科学出版社.2001年.

附录A 1.复化梯形求积公式的程序设计：

(1).源程序：

#include

#include

double f(double x)

{double z;

z=sqrt(1+pow((2*x-3*pow(x,2)),2));

return z;

}

main()

{ int n,k;

float h;

float a;

float b;

double s=0.0;

double s1=0.0;

double t;

printf("Please input the deng fen ;"); scanf("%d",&n);

printf("Please input qujian a ;");

scanf("%f",&a);

printf("Please input qujian b ;");

scanf("%f",&b);

h=(b-a)/n;

for (k=1;k

{ t=a+k*h;

s1=s1+f(t);}

s=(h/2)*(f(a)+2*s1+f(b));

printf("%f\n",s);

}

(2).运行结果：

图1 复化梯形求积公式计算弧长结果

附录B

2.复化辛普森求积公式的程序设计：

(1).源程序：

#include

#include

double f(double x)

{double z;

z=sqrt(1+pow((2*x-3*pow(x,2)),2));

return z;}

main()

{ int n,i,j;

float h;

float a;

float b;

double s=0.0;

double s1=0.0,s2=0.0;

double t,l;

printf("Please input the deng fen ;");

scanf("%d",&n);

printf("Please input qujian a ;");

scanf("%f",&a);

printf("Please input qujian b ;");

scanf("%f",&b);

h=(b-a)/n;

for(i=1;i<8;i=i+2)

{t=a+i*h/2;

s1=s1+4*f(t);}

for(j=2;j<8;j=j+2)

{l=a+j*h/2;

s2=s2+2*f(l);}

s=(h/6)*(f(a)+s1+s2+f(b));

printf("%f\n",s);

}

(2).运行结果：

图2 复化辛普森求积公式计算弧长结果

附录C

3.复化柯特斯求积公式的程序设计：

(1).源程序：

#include

#include

double f(double x)

{double z;

z=sqrt(1+pow((2*x-3*pow(x,2)),2));

return z;

}

main()

{

int n,i,j,k;

float h;

float a;

float b;

double s=0.0;

double s1=0.0,s2=0.0,s3=0.0;

double t,l,m;

printf("Please input the deng fen ;");

scanf("%d",&n);

printf("Please input qujian a ;");

scanf("%f",&a);

printf("Please input qujian b ;");

scanf("%f",&b);

h=(b-a)/n;

for(i=1;i<8;i=i+2)

{

t=a+i*h/4;

s1=s1+32*f(t);

}

for(j=2;j<7;j=j+4)

{l=a+j*h/4;

三种复化求积分算法的精度分析

【摘要】分别利用复化梯形公式、复化simpson公式和复化gauss-legendre i型公式对定积分进行运算，得到近似数值解，并对各算法的精度和计算复杂度进行了比较与分析。数值举例结果表明，三种复化求积分算法的运算结果均在绝对误差限ε=5e-8内，并且在相同的精度下，复化gauss-legendre i型公式的步长和计算量最小。 【关键词】复化梯形公式；复化simpson公式；gauss-legendre公式 1 引言 数值积分是计算数学的基本内容，在工程技术和科学计算中起着十分重要的作用，当积分的精确值不能不能求出时，数值积分就变得越来越重要。通常数值积分的计算常利用机械积分来实现，其基本思想为： （1） 2 理论模型 复化梯形求积公式 将区间[a，b]划分成n等分，分点xk=a+kh（，k=1，2，3…n），在每个子区间[xk，xk+1] （k=1，2，3 …n-1）上采用梯形式，则得到 （2） 记 （3） 上式（3）为复化梯形公式，其余项可由式 ，（a≤η≤b）（4） 得 ，ηk∈[xk，xk+1] （5） 由于 f（x）∈c2[a，b] 且 ，（0≤k≤n-1）（6） 所以∈（a，b），使 （7） 于是复化梯形公式余项为 （8） 复化simpson求积公式 将区间[a，b]划分为n等分，在每个子区间[xk，xk+1]上采用simpson式，若记，则得 （9） 记 （10） 上式（10）为复化simpson求积公式，其余项可由式 ，（a≤η≤b）（11） 得 ，ηk∈[xk，xk+1] （12） 于是当f（x）∈c4[a，b]时，与复化梯形公式相似有 ，η∈[a，b] （13） 复化gauss-legendre i型求积公式 gauss型求积公式是具有最高代数精度的插值求积公式。通过适当选取求积公式（1）的节点ε=5e-8和求积系数ak≥0和xk∈[a，b] （k=1，2，3…n），可使其代数精度达到最高的2n+1次。利用特殊区间[-1，1]上n+1次legendre正交多项式的根作为节点，我们可以建立gauss-legendre型求积公式。将区间[a，b]划分成n等分，分点xk=a+kh（，

编程MATLAB程序实现复化梯形和辛普森数值积分

数值分析实验报告—— 实验目的[1] 掌握复化梯形和辛普森数值积分法的基本原理和方法； [2] 编程MA TLAB程序实现复化梯形和辛普森数值积分 实验内容与步骤1.编程序实现复化梯形数值积分求积公式 function y=f(x) y=sqrt(x).*log(x); function T_n=F_H_T(a,b,n) h=(b-a)/n; for k=0:n x(k+1)=a+k*h; if x(k+1)==0 x(k+1)=10^(-10); end end T_1=h/2*(f(x(1))+f(x(n+1))); for i=2:n F(i)=h*f(x(i)); end T_2=sum(F); T_n=T_1+T_2;

实验内容与步骤运行结果： >> T_n=F_H_T(0,1,20) T_n = -0.4336 2.编程序实现复化辛普森数值积分求积公式 function y=f(x) y=sqrt(x).*log(x); function S_n=S_P_S(a,b,n) h=(b-a)/n; for k=0:n x(k+1)=a+k*h; x_k(k+1)=x(k+1)+1/2*h; if (x(k+1)==0)|(x_k(k+1)==0) x(k+1)=10^(-10); x_k(k+1)=10^(-10); end

S_1=h/6*(f(x(1))+f(x(n+1))); for i=2:n F_1(i)=h/3*f(x(i)); end for j=1:n F_2(j)=2*h/3*f(x_k(j)); end S_2=sum(F_1)+sum(F_2); S_n=S_1+S_2; 运行结果： >> S_n=S_P_S(0,1,20) S_n = -0.4423 实验心得 通过此次实验的操作，我掌握了复合梯形公式和复合辛普森公式，对编程又有了新的突破！

数值分析作业复化求积公式

数值计算方法上机题目3 计算定积分的近似值： 2 2 1x e xe dx =? 要求： （1）若用复化梯形公式和复化Simpson 公式计算，要求误差限7102 1-?=ε，分别利用他们的余项估计对每种算法做出步长的事前估计； （2）分别利用复化梯形公式和复化Simpson 公式计算定积分； （3）将计算结果与精确解比较，并比较两种算法的计算量。 解： （1) x xe x f =)(,所以x x k xe ke x f +=)()(，x x xe e x f +=2)(''，x x xe e x f +=4)()4( x x xe e x f +=6)()6( 对于复化梯形公式： )(12)(''2ηf h a b f R n --=，2max ''4)(e f =η，n h 1= 代入数据可知 722102 1124-?≤n e ，57.7018≥n 取7019=n 对于复化Simpson 公式 )()2(180)()4(4ηf h a b f R n --=，2max )4(6)(e f =η，n h 1= 代入数据可知 742102 128806-?≤n e ，56.23≥n 取24=n （2）复化梯形公式： 函数 function y=fun(x) y=x*exp(x); 程序： clc Clear % 复化梯形计算 format long

a=1;b=2; n=7019;% 区间划分为m等份 h=(b-a)/n;% 步长，根据误差限由该算法的余项作事前估计得到 ty1=fun(a)+fun(b); ty2=0; for i=1:n-1 x=a+i*h1; ty2=ty2+fun(x); end T=h*(ty1+2*ty2)/2; T 对于复化Simpson公式 clc Clear % 复化Simpon计算 format long a=1;b=2; n=24;% 区间划分为n等份 h=(b-a)/(2*n);% 步长，根据误差限由该算法的余项作事前估计得到sy1=fun(a)+fun(b); sy2=0;sy3=0; for j=1:2*n-1 x=a+j*h2; if rem(j,2)==0 sy3=sy3+fun(x); else sy2=sy2+fun(x); end end S=h*(sy1+4*sy2+2*sy3)/3; S % 精确值 Exactanswer=exp(2) 运算结果： T = 7.389056127230221 S = 7.389056126214707 Exactanswer = 7.389056098930650 （3）比较可知复化梯形公式的计算量较大

辛普森求积公式

摘要 在工程实验及研究中，实际工作中，变量间未必都有线性关系，如服药后血药浓度与时间的关系；疾病疗效与疗程长短的关系；毒物剂量与致死率的关系等常呈曲线关系。曲线拟合是指选择适当的曲线类型来拟合观测数据，并用拟合的曲线方程分析两变量间的关系.可以说，曲线拟合模型与我们的生活生产密切相关. 本课题着重介绍曲线拟合模型及其应用，其中包括它的基本思想、模型的建立、以及具体应用.为了更好的了解曲线拟合模型，可以将它分为线性与非线性模型，在模型建立的基础上我们可以用最小二乘法来解决一些我们日常所应用的问题. 关键词曲线拟合；线性与非线性模型；最小二乘发

目录 引言 (1) 第一章曲线拟合 (2) §1.1 基本思想及基本概念 (2) §1.1.1 方法思想 (2) §1.1.2几个基本概念 (2) §1.2辛普森算法基本定义及其应用 (4) §1.2.1辛普森求积公式的定义 (4) §1.2.2辛普森求积公式的几何意义 (5) §1.2.3辛普森求积公式的代数精度及其余项 (5) §1.2.4辛普森公式的应用 (6) 第二章辛普森求积公式的拓展及其应用 (7) §2.1 复化辛普森求积公式 (7) §2.1.1问题的提出 (7) §2.1.2复化辛普森公式及其分析 (7) §2.1.3复化辛普森公式计算流程图 (8) §2.1.4复化辛普森公式的应用 (9) §2.2 变步长辛普森求积公式 (10) §2.2.1变步长辛普森求积公式的导出过程 (10) §2.2.2变步长辛普森求积公式的加速过程 (12) §2.2.3变步长辛普森求积公式的算法流程图 (13) §2.2.4变步长辛普森公式算法程序代码 (14) §2.2.5变步长辛普森求积公式的应用 (14) §2.2.6小结 (14) §2.2.7数值求积公式在实际工程中的应用 (14) 参考文献 (16) 附录A (17)

复化积分法(复化梯形求积-复化Simpson公式-变步长求积法)MATLAB编程实验报告 (1)

复化积分法（复化梯形求积，复化Simpson 公式，变步长求积法） MATLAB 编程实验报告 一、 问题描述： 编写函数实现复化积分法。 二、 实验步骤（过程）： （一）复化求积法 （1）复化梯形求积：用复化梯形求积公式求解 dx x x ?10sin function [f]=Tn(a,b,n,y) syms t; h=(b-a)/n; f=0; for k=1:n+1 x(k)=a+(k-1)*h z(k)=subs(y,t,x(k)); end for i=2:n f=f+z(i); end q=subs(y,t,a); if y=='sin(t)/t'&&a==0 q=1; end p=subs(y,t,b); T=h/2*(q+p+2*f); T=vpa(T,7) clc,clear; syms t; a=0;b=1; y=sin(t)/t; n=8; Tn(a,b,n,y); （2）复化Simpson 公式：用复化Simpson 公式求解?211dx e x function [f]=simpson(a,b,n,y)

syms t; h=(b-a)/n; f=0;l=0; for k=1:n+1 x(k)=a+(k-1)*h w(k)=0.5*h+x(k) z(k)=subs(y,t,x(k)); end for i=2:n f=f+z(i); end for i=1:n l=l+w(i); end q=subs(y,t,a); if y=='sin(t)/t'&&a==0 q=1; end p=subs(y,t,b); T=h/2*(q+p+2*f); T=vpa(T,7) clc,clear; syms t; a=1;b=2; y=exp(1/t); n=5; simpson(a,b,n,y); （3）变步长求积法：以书本例4.5为例function [f]=TN(a,b,y,R0) syms t; T=[]; f=0; q=subs(y,t,a); if y=='sin(t)/t'&&a==0 q=1; end p=subs(y,t,b); T(1)=(b-a)/2*(q+p); i=2; n=i-1; h=(b-a)/n; z1=a+h/2; z2=subs(y,t,z1);

数值分析与实验复化辛卜生公式龙贝格算法

数值分析与实验课程设计 班级: 姓名: 学号:

08级应用数学《数值分析与实验（实践）》任务书 一、设计目的 通过《数值分析与实验（实践）》实践环节，掌握本门课程的众多数值解法和原理，并通过编写C语言或matlab程序，掌握各种基本算法在计算机中的具体表达方法，并逐一了解它们的优劣、稳定性以及收敛性。在熟练掌握C语言或matlab语言编程的基础上，编写算法和稳定性均佳、通用性强、可读性好，输入输出方便的程序，以解决实际中的一些科学计算问题。 二、设计教学内容 1、数值方法的稳定性； 2、禾U用牛顿法和割线法程序求出非线性方程的解，并比较它们之间的优 劣； 3、高斯消去法和列主元高斯消去法求解线性方程组； 雅克比法和高斯-赛德尔迭代法解方程组； 4、利用Lagrange插值多项式求未知点的近似值； 5、利用所给数据进行数据的多项式和可转化成多项式形式的函数拟合； 6、编写复化辛卜生公式和龙贝格算法，通过实际计算体会各种方法的精确 \ 度； 7、利用改进Euler方法和四阶Runge-Kutta方法求解初值问题的微分方程组； &利用幕法求矩阵按模最大的特征值及对应特征向量； \ （8个中选取1个） 二、设计时间 2011 —2012学年第1学期：第16周共计一周 教师签名: 2011年12月12日

、八 刖 数值计算方法是一种利用计算机解决数学 .言 问题的数值近似解方法， 特别是无法用人工过计算器计算的数学问题。数值计算方法常用于矩阵高次代数方程矩阵特征值与特征向量的数值解法，插值法，线性方程组迭代法，函数逼近，数值积分与微分，常微分方程初值问题数值解等。 作为数学与计算机之间的一条通道，数值计算的应用范围已十分广泛，作为用计算机解决实际问题的纽带，数值算法在求解线性方程组，曲线拟合、数值积分、数值微分，迭代方法、插值法、拟合法、最小二乘法等应用广泛。 数值计算方法是和计算机紧密相连的，现代计算机的出现为大规模的数值计 算创造了条件，集中而系统的研究适用于计算机的数值方法是十分必要的。数值计算方法是在数值计算实践和理论分析的基础上发展起来的。 通过数值计算方法与实验将有助于我们理解和掌握数值计算方法基本理论和相关软件的掌握，熟练求解一些数学模和运算，并提高我们的编程能力来解决实际问题。

计算方法期中测试(二)答案

期中测试（下） 班级： 姓名： 学号： 分数： 一、填空题（20分） 1、计算积分1 ? 2、5个节点的牛顿- 3、求积公式 0()()b n k k k a f x dx A f x =≈∑? 4 、数值积分公式 ()1 1 ()29[(1)8(0)(1)]f x dx f f f -'≈-++? 5、求解一阶常微分方程初值问题00(,),()y f x y y x y '==的改进欧拉公式为 二、选择题（6分） 1、舍入误差是( A )产生的误差。 A. 只取有限位数 B ．模型准确值与用数值方法求得的准确值 C ． 观察与测量 D ．数学模型准确值与实际值 2、用 1+x 近似表示e x 所产生的误差是( C )误差。 A ． 模型 B ． 观测 C ． 截断 D ． 舍入 3、解线性方程组的主元素消去法中选择主元的目的是( A )。 A ．控制舍入误差 B ． 减小方法误差 C ．防止计算时溢出 D ． 简化计算 4、求解初值问题00(,),()y f x y y x y '==欧拉法的局部截断误差是( A )；中心欧拉法的局部截断误差是( B )； 改进欧拉法的局部截断误差是( B )；四阶龙格－库塔法的局部截断误差是( D ) A. 2()O h B. 3()O h C. 4()O h D. 5()O h 三、计算题（64分） 1、（10分）试分别推导复化梯形和复化辛普森求积公式。 证明：以积分 ()b f x dx 为例。将积分区间[,]a b 做n 等分，步长()/h b a n =-。

2、(10分) 求A 、B 使求积公式1 1 ()[(1)(1)][(0.5)(0.5)]f x dx A f f B f f -≈-++-+? 的代数精度尽量高, 并求其代数精度； 利用此公式求2 1 1 I dx x = ? (保留4位小数)。 解：2 ,,1)(x x x f =是精确成立，即 3、(12分) 取5个等距节点 ，分别用复化梯形公式和复化辛普森公式计算积分2 201 I dx = ?的近似值（保留4位小数）。

复化抛物线积分公式

摘要 求函数在给定区间上的定积分，在微积分学中已给出了许多计算方法，但是，在实际问题计算中，往往仅给出函数在一些离散点的值，它的解析表达式没有明显的给出，或者，虽然给出解析表达式，但却很难求得其原函数。这时我们可以通过数值方法求出函数积分的近似值。 在用近似值代替真实值时，遇到的问题就是近似值的代数精度是否足够。当代数精度不足够时，很显然提高插值函数的次数是一种方法，但是考虑到数值计算的稳定性，当次数过高时，会出现龙格现象，用增大n的方法来提高数值积代数精度是不可取的。因此，提出类似于分段插值，为了减少数值积分的误差，可以把积分区间分成若干个小区间，在每个小区间上采用低阶数值积分公式，然后把这些小区间上的数值积分结果加起来作为函数在整个区间上的近似值，这个就是复化数值积分的思想。 本实验针对在每个小区间上利用抛物线积分公式，即阶数为2，进行实验。 关键词：龙格现象复化数值积分代数精度复化抛物线积分公式

1、实验目的 1)通过本次实验体会并学习复化抛物线积分公式的优点。 2)通过对复化抛物线积分公式进行编程实现，提高自己的编程能力。 3)用实验报告的形式展现，提高自己在写论文方面的能力。 2、算法流程 已知定积分的抛物线积分公式及其误差为 根据数学知识，我们知道积分区间可划分，且不改变积分值，即如下所示：针对上式，在每一个小区间上利用抛物线积分公式有 得到 其中,令

当作为积分的近似值时，其误差为 若,则由介值定理推得 设,得到误差限 由上式可以进行计算精度控制。这样就给出了n+1点复化抛物线积分公式及其误差 3、算法实例 用复化抛物线积分公式计算积分 解：具体程序如下：

教案二 复化求积公式、变步长求积方法

教案二 复化求积公式、变步长求积方法 基本内容提要 1 复化求积公式的基本思想 2 基本复化求积公式及其余项公式 3 复化求积公式收敛阶的概念 4 变步长求积方法的基本思想与计算过程 教学目的和要求 1 理解复化求积思想 2 掌握低阶复化牛顿－科特斯求积公式及其余项公式 3 理解复化牛顿－科特斯求积公式的收敛阶概念 4 理解变步长求积方法的基本思想方法、掌握变步长求积方法及其优越性 教学重点 1 复化求积公式的基本思想与基本公式 2 变步长求积方法的基本思想与计算过程 教学难点 1 低阶复化牛顿－科特斯求积公式及其余项公式的应用 2 利用变步长求积方法的基本思想进行数值积分计算的过程及其优越性 课程类型 新知识理论课 教学方法 结合提问，以讲授法为主 教学过程 问题引入 通过分析各低阶Newton-Cotes 求积公式的余项公式，揭示其与积分区间长度的关系，引出复化求积的思想。 §3.2复化求积公式 典型的复化求积公式有复化梯形求积公式和复化辛普森求积公式。 复化梯形求积公式如下推导：选定.n a b h ?= ∫∫ ∫∫?+++==2 1 1 1 )()()()(][x x x x x x b a n n dx x f dx x f dx x f dx x f f I L

≈ ))]()(())()(())()([(2 12110n n x f x f x f x f x f x f h ++++++?L =.)()](21 )()(21[0 h T b f x f a f h n i i ∧++∑= )(h T 即为复化梯形求积公式。其截断误差 ∑=?+?=?=n i i i x f x f h T f I f R 1 1)]()([)()(][ =),()(12 2 ηf a b h ′′?? 其中，],[b a ∈η满足条件 .) ()()()(21n f f f f n ηηηη′′++′′+′′= ′′L 类似于上述过程，可以推导复化辛普森求积公式：选定.2n a b h ?= ∫∫ ∫∫?+++==4 21 22 22 )()()()(][x x x x x x b a n n dx x f dx x f dx x f dx x f f I L = ))]()(4)(())()(4)(())()(4)([(3 21222432210n n n x f x f x f x f x f x f x f x f x f h +++++++++??L =.)()]()()(4)([311 11 212h S b f x f x f a f h n i n i i i ∧+++∑∑=?=? .)(h S 即为辛普森求积公式。其截断误差 ),()(180)(][][44 ηf a b h h S f I f R ???= 其中，].,[b a ∈η 例3.2.1 利用梯形求积公式计算： dx x x ∫1 0sin 使截断误差不超过3105.0?×.

3.2复化求积公式习题及解答

3.2-3.5习题 一、填空题 1. 梯形求积公式和复化梯形公式都是插值型求积公式_____（对或错）。 (答案：错) 2.已知(1)1.2,(2)1.4,f f f = ==，则用复合梯形公式计算求得 3 1 ()f x dx ≈? ， （答案：2.75） 3. 已知，在[0, 1] 内 ，有一位整数，用复合 梯形求积公式计算要保证有3位有效数字，至少应将[0, 1]（ ）等分。 A. 4 B. 5 C. 6 D. 7 4、(1)1,(2)2,(3)2f f f -===-，则[1,2,3]f -=_________，三点高斯求积公式 2 ()f x dx ≈? ______________. 答案： )531(95)1(98)531(95,1213+++-- f f f 二、计算题 1．建立Gauss 型求积公式111220 ()()A f x A f x ≈+? 答案 12120.0455363610.6421930581.035301293 0.964698706x x A A ==== 2． 试确定常数A ，B ，C 和，使得数值积分公式 有尽可能高的代数精度。试问所得的数值积分公式代数精度是多少？它是否为Gauss 型的？ 答案： ， 该数值求积公式具有5次代数精确度，它是Gauss 型的。 3、用Romberg 算法计算积分3 0?（只作两次外推）。 解：取2,0,3t a b === (0)11 (()())14.230249472 T f a f b =+=， 1t =，0 2(1)(0)1 10 1 111 (()11.17136992 222i T T f a b a =??=++ -= ? ?? ? ∑ ，

辛普森公式

Simpson算法及其推广形式 摘要：本文研究了辛普森公式的数值积分的计算方法问题，并且更进一步研究了变步长复化的辛普森公式和二重积分的辛普森公式的问题。首先是对 一维辛普森公式和变步长复化辛普森公式以及二维辛普森公式的推导及 其算法，进行误差分析，并且列举了实例。然后，对辛普森公式进行改 进，这里的改进最主要是对辛普森公式的代数精度进行提高，从而使辛 普森公式对积分的计算更加精确。另外，还研究了辛普森公式的推广形 式。最后，在结论的当中列举了一个例子。 关键词：辛普森公式算法改进推广形式二重积分的辛普森公式

Abstract：This paper first studies the calculation methods of the numerical integration in simpson formula, and then study of the long-simpson formula and the double integral simpson formula problem. First, study the algorithm and derived of one-dimensional simpson formula and step-change in simpson formula, as well as two-dimensional simpson formula, and then analysis the error. Finally , list the example. In this , improve the simpson formula. This improved the most important is to incre ase the simpson formula’s accuracy of algebra. Besides, we study the simpson formula’s promotion of forms. At the last, we list a example in the conclusion. Key word:The simpson formula, Algorithm, Improve, Promotion of forms, The simpson formula of the two-dimensional integral.

复式求积法(计算方法实验报告)

实验名称：实验6 复化求积法 实验题目：测得飞机在高度h 时上升速度v （h ）的数据如表所列： 飞机从地面上升到H km 高度所需时间可用公式 dh h v 1t H 0?=）（ 分别用复化梯形公式和复化Simpson 公式计算飞机上升到8km 高空所需时间。 实验目的：用复化求积法解决实际问题 基础理论：复化求积公式：复化梯形公式和复化Simpson 公式 实验环境：操作系统：windows 7 实验平台：Visual C++ 实验过程：根据教材书本P-86-92复化求积公式的概念和复化梯形公式和复化Simpson 公式以及算法，在VS2010中编写C++代码，调试并运行程序。 结果分析： 附 录： 程序清单： #include "stdafx.h"

int _tmain(int argc, _TCHAR* argv[]) { return 0; } #include #include #define n 4; using namespace std; void main() { int i; float h,temp,xk,xm,xr,yl,ym,yr; float B[6]={50.0,46.0,40.0,32.2,22.5,10.0}; cout<

复化梯形法复化矩形法变步长梯形变步长辛普森

陕西科技大学 机械教改班 用C++的积分 其实积分的思想就是,微分—>求和—>取极限，如果是用纯手工法那就是先对一个函数微分，再求出它的面积，在取极限，因为我们的计算速度和计算量有限，现在有了计算机这个速度很快的机器，我们可以把微分后的每个小的面积加起来，为了满足精度，我们可以加大分区，即使实现不了微分出无限小的极限情况，我们也至少可以用有限次去接近他，下面我分析了四种不同的积分方法，和一个综合通用程序。 一．积分的基本思想 1、思路：微分—>求和—>取极限。 2、Newton —Leibniz 公式 ?-=b a a F b F dx x f ) ()()( 其中，)(x F 被积函数 )(x f 的原函数。 3、用计算机积分的思路 在积分区间内“微分—>求和—>控制精度”。因为计算机求和不可以取极限，也就是不可以无限次的加下去，所以要控制精度。 二．现有的理论 1、一阶求积公式---梯形公式 ?=+-=b a T b f a f a b dx x f )]()([2 )( 他只能精确计算被积函数为0、1次多项式时的积分。 2、二阶求积分公式——牛顿、科特斯公式 ?=+++-=b a S b f a b f a f a b dx x f )]()2(4)([6)( 他只能精确计算被积函数为0、1、2、3次多项式时的积分。 三．四种实现方法 1.复化矩形法 将积分区间[a,b]等分成n 个子区间: ],[],[],[],[],[112322110n n n n x x x x x x x x x x ---、、、 则h=(b-a)/n,区间端点值k x =a+kh

6.3 复化求积公式

§3 复化求积公式 ● 复化求积法的基本思想： 将积分区间],[b a n 等分，可得到1+n 个求积节 点：kh a x k +=，),,1,0(n k Λ=，其中n a b h -=，对 积分1 1 1 ()()k k n n b x k a x k k I f x dx f x dx I +--==== =∑∑?? 在每一个小区间1[,]k k x x +上利用n 阶牛顿-柯特斯公式计算，然后对每个区间的近似积分值求和，用所得的值近似代替原积分值。如此得到的求积公式称为复化求积公式。 ● 复化梯形公式：（每个小区间上利用梯形公式求积） 1 1 1 110 ()()(()()) 2k k n b x a x k n k k k k k I f x dx f x dx x x f x f x +-=-++===-≈+∑?? ∑ 求和展开得：

011211 1 (()())(()()) 2 (()()) (()2()())2n n n n k k h T f x f x f x f x f x f x h f a f x f b --==++++++=++∑L 其中，n a b h -= 复化辛甫生公式: （每个小区间上用辛甫生公式求积） 1、公式： 1 12 10 1 110 ()()(()4()())6k k n b x a x k n k k k k k k I f x dx f x dx x x f x f x f x +-=-+++===-≈++∑?? ∑ 12 k x +表示为区间1[,]k k x x +的中点。 求和展开得： 1322 12 01121((()4()())(()4()6 ())(()4()())n n n n h S f x f x f x f x f x f x f x f x f x --=+++++++++L

复化梯形求积公式

第二章 1.1 复合梯形求积公式 复合梯形求积公式是复合求积法的一种，在本章中，将从其原理、概念等方面对它做一个详细介绍。在本章的最后，会对复合梯形求积法进行程序设计，使得可以从不同的方面对这种方法有更深的理解。 1.1.1 复合梯形求积公式的理论 当积分区间[a ，b]的长度较大，而节点个数1+n 固定时，直接使用Newton-Cotes 公式的余项将会较大。但是如果增加节点个数，即1+n 增加时，公式的舍入误差又很难得到控制。为了提高公式的精度,又使算法简单易行,往往使用复化方法。 即将积分区间][b a ， 分成若干子区间，然后在每个小区间上使用低阶Newton-Cotes 公，最后将每个小区间上的积分的近似值相加，这就叫做复合求 积法。而复合梯形求积公式就是复合求积法的一种。 1.1.2复合求积公式的原理 将区间[]b a ,划分为n 等分，分点,,,1,0,,n k n a b h kh a x k =-= += 在每个子区间[](),1,,1,0,1-=+n k x x k k 上采用梯形公式，则得 []) ()()(2)()(11 1 1 f R x f x f h dx x f dx x f I n k n k k b a n k x x k k ++===+-=-=∑?∑? + 记 ()[()]()[()()]∑∑-=+-=++=+=1 1 110222n k b k k n k k n x f x f a f h x f x f h T ， （1.1） 称为复合梯形公式，其余项可由 )(). ,(),(12 ][''3 b a f a b f R ∈-- =ηη 得 ()()() 11 0''3,12+-=∈?? ? ???-=-=∑k k k n k k n n x x f h T I f R ηη 由于[],,)(2b a C x f ∈ 且 ()(),max 1min 1010 ' '' '10-≤≤-=-≤≤≤≤∑n k k n k k n k f n f ηη

数值积分的辛普森方法

实习七 数值积分的辛普森方法 一、实习目的 1．掌握计算定积分近似值的辛普森方法； 2．理解复化辛普森求积公式。 二、相关知识 抛物线公式（辛普森公式） 将积分区间],[b a 作2n 等分：n i ih a x n a b h i 2,,2,1,0,,2 =+=-=，现在考察由分点22-k x 和k x 2形成的一个小区间],[222k k x x -，（12-k x 为中点），n k ,,2,1 =，在每一个 小区间],[222k k x x -上，作一条抛物线k k k x x y γβα++=2通过三点))(,(2222--k k x f x ， ))(,(1212--k k x f x 和))(,(22k k x f x ，这样就产生关于未知系数k α，k β和k γ的线性方程组 ?????=++=++=++------)() ()(222212122122222222k k k k k k k k k k k k k k k k k k x f x x x f x x x f x x γβαγβαγβα （7-1） 显然上述方程组有唯一解（由高等代数知识知）。 现在，以)(x f y =为顶的曲边梯形用以抛物线k k k x x y γβα++=2为顶的曲边梯形来 代替，其面积 dx x x dx x f k k n k x x k b a k k )()(1 2222γβα++≈∑??=-∑=--=n k k k x x 12226]4)(2)2([){(2222222222222222k k k k k k k k k k k k k x x k x x x x x x γβαγβα++++++++----- }222k k k k k x x γβα+++)]()(4)([6212221222k k k n k k k x f x f x f x x ++-=--=-∑ （7-2） 得抛物线公式，记为n S 2，化简后： {})()(4)(2)(4)(2)(4)(6212432102n n n x f x f x f x f x f x f x f n a b S +++++++-=- 在实际求解数值积分时，我们总是采用成倍加密节点的方法，就抛物线公式而言，若n S 2被认为精度不够，则接着计算n S 4，而精度是否达到要求，又以n n S S 24-是否足够小作为判

复化求积法

复化求积法 1.实验题目 H1/v(h)dh t=∫ 分别用复化梯形公式和复化Simpson公式计算飞机上升到8km高空所需的时间. 2.实验目的 验证复化梯形公式和复化Simpson公式对求积分的数值结果。用C++程序设计语言编程实现这一算法，加深对其的优缺点的理解。 3.基础理论 4.实验环境 Visual C++ 语言 5.实验过程 #include #define f(v) 1/v double g(double x) { double h[6]={0,2,4,6,8,10}; double v[6]={50,46,40,32.2,22.5,10}; for(int i=0;i<6;i++) { if(x==h[i]) return v[i]; if(x

double a=0,b=8; int n=4; int h=2; double t=0; double x=a; double xx=g(x); t=t+f(xx); x=a+h/2; xx=g(x); t=t+4*f(xx); for(int i=0;i