数值积分用matlab实现.doc

.

东北大学秦皇岛分校
数值计算课程设计报告

数值积分及Matlab实现

学 院
数学与统计学院
专 业
信息与计算科学
学 号
5133117
姓 名
楚文玉
指导教师
张建波 姜玉山
成 绩

教师评语：

指导教师签字：

2015年07月14日

.
.
1 绪论
在科研计算中，经常会碰到一些很难用公式定理直接求出精确解的积分问题，对于这类问题，我们一般转化为数值积分问题，用计算机来实现求解问题．
1.1 课题的背景
对于定积分在求某函数的定积分时，在一定条件下，虽然有牛顿-莱布里茨公式可以计算定积分的值，但在很多情况下的原函数不易求出或非常复杂．被积函数的原函数很难用初等函数表达出来，例如等；有的函数的原函数存在，但其表达式太复杂，计算量太大，有的甚至无法有解析表达式．因此能够借助牛顿-莱布尼兹公式计算定积分的情形是不多的．另外，许多实际问题中的被积函数往往是列表函数或其他形式的非连续函数，对这类函数的定积分，也不能用不定积分方法求解，只能设法求其近似值．因此，探讨近似计算的数值积分方法是有明显的实际意义的，即有必要研究定积分的数值计算方法，以解决定积分的近似计算．而数值积分就是解决此类问题的一种有效的方法，它的特点是利用被积函数在一些节点上的信息求出定积分的近似值．微积分的发明是人类科学史上一项伟大的成就，在科学技术中，积分是经常遇到的一个重要计算环节数值积分是数学上重要的课题之一，是数值分析中重要的内容之一．随着计算机的出现，近几十年来，对于数值积分问题的研究已经成为一个很活跃的研究领域．现在，数值积分在计算机图形学，积分方程，工程计算，金融数学等应用科学领域都有着相当重要的应用，所以研究数值积分问题有着很重要的意义．国内外众多学者在数值积分应用领域也提出了许多新方法．在很多实际应用中，只能知道积分函数在某些特定点的取值，比如天气测量中的气温、湿度、气压等，医学测量中的血压、浓度等等．通过这个课题的研究，我们将会更好地掌握运用数值积分算法求出特殊积分函数的定积分的一些基本方法、理论基础；并且通过Matlab软件编程的实现，应用于实际生活中．
1.2 课题的主要内容框架
1.2.1 数值积分各求积公式简介
简介牛顿-柯特斯求积公式及其辛普森求积公式，龙贝格求积公式，高斯求积公式的基本理论基础和方法．
1.2.2 求积公式的代码实现
通过理解各种数值积分求积公式的原理方法，通过Matlab软件编程，实现以上求积公式．
1.2.3 应用举例
通过简单举例，自建一个相对简单和复杂的函数，用上面编写的Matlab源程序来解决实际问题，体会数值积分和Matlab的优势．
2 牛顿－柯特斯公式及Matlab实现
2.1 牛顿－柯特斯公式的基本原理方法
设将积分区间[a, b]划分为n等分，步长为，选取等距节点构造出的差值型求积公式
， （2.1）
称为牛顿-柯特斯公式，式中称为柯特斯系数．根据
， （2.2）
引进变量代换，则有
（2.3）
当n = 2时，此时柯特斯系数为，，，相应的求积公式就是辛普森求积公式：
（2.4）
2.2 牛顿－柯特斯公式的Matlab实现
function[C, g] = NCotes(a, b, n, m)
% a，b分别为积分的上下限；
% n是子区间的个数；
% m是被调用第几个被积函数；
% 当n=1时计算梯形公式；当n=2时计算辛普森公式，以此类推；
I = n;
h = (b - a) / i;
z = 0;
for j = 0 : i
x(j + 1) = a + j * h;
s = 1;
if j == 0
s = s;
else
for k=1 : j
s =s * k;
end
end
r = 1;
if i - j == 0
r = r;
else
for k = 1 : (I - j)
r = r * k;
end
end
if mod((I - j), 2) == 1
q = -(I * s * r);
else
q = i * s * r;
end
y = 1;
for k = 0 : i
if k ~= j
y = y * (sym('t') - k);
end
end
l = int(y, 0 , i);
C(j + 1)= l / q;
z = z + C(j + 1)*f1(m, x(j + 1));
end
g=(b - a)*z
3 复合求积公式及Matlab实现
3.1 复合梯形公式的基本原理
将区间[a, b]划分成n等分，分点，,在每个子区间[]（）上采用梯形公式得：
（3.1）
记
（3.2）
称式(3.2)为复合梯形公式．
3.2 复合梯形公式的Matlab实现
function s = trapr1(f, a, b, n)
% f表示被积函数；
% a，b表示积分上下限；
% n是子区间的个数；
h = (b - a) / n;
s = 0;
for k = 1 : (n - 1)
x = a + h * k;
s = s + feval('f', x);
end
format long
s = h*(feval('f', a) + feval('f', b)) / 2 + h * s;
3.3 复合辛普森求积公式的基本原理
将区间[a,b]分等分，在每一个子区间]上采用辛普森公式，若记，则得
（3.3）
记
（3.4）
称式(3.4)为复合辛普森求积公式．
3.4 复合辛普森求积公式的Matlab实现
function s = simpr1(f, a, b, n)
% f表示被积函数；
% a，b表示积分上下限；
% n是子区间的个数；
h = (b - a) / (2 * n);
s1 = 0;
s2 = 0;
for k = 1 : n
x = a + h * (2*k - 1);
s1 = s1 + feval('f', x);
end
for k = 1 : (n - 1)
x = a + h * 2 * k;
s2 = s2 + feval('f', x);
end
s = h*(feval('f', a) + feval('f', b) + 4 * s1 + 2 * s2) / 3;
4 龙贝格求积公式及Matlab实现
4.1 龙贝格算法的基本原理
由梯形的递推法可以看出，将积分区间等分时，用复化梯形公式计算的结果作为积分的近似值，其误差近似值为
（4.1）
可以设想，如果用这个误差作为的一种补偿，即将
（4.2）
作为积分的近似值，可望提高其精确度．直接根据复化求积公式，不难验证
（4.3）
这说明，将区间对分前后两次复化梯形公式的值，按式
（4.4）
作线性组合恰好等于复合辛普森公式的值，它比更接近于近似值． 同样，根据
（4.5）
用于作线性组合会得到比更精确的值，且通过直接验证可得
（4.6）
再由
（4.7）
用与作线性组合，又可得到比更精确的值，通常记为，即
（4.8）
此式（4.8）就称为龙贝格求积公式．
上述用若干积分近似值推算出更为精确的积分近似值得方法，称为外推法．我们将
序列分别称为梯形序列，辛普森序列，柯特斯序列和龙贝格序列．由龙贝格序列求积的算法称为龙贝格算法．具体步骤为：
第一步：算出和的值，根据公式
（4.9）
求出;
第二步：将区间[a,b]分半，算出的值，并根据(4.3)和(4.9)式计算和;
第三步：再将区间分半，算出及的值，并根据(4.3)和(4.9)式计算和，再有公式(4.4)求出；
第四步：再将区间分半，计算，, ,并根据公式(4.5)计算．
第五步：再将区间分半，类似上述过程计算，, ，．
重复以上步骤即可得到， ,一直到龙贝格序列中前后两项的绝对值差不超过给定的误差险为止．
4.2 龙贝格算法的Matlab实现
function [R, quad, err, h]=romber(f, a, b, n, delta)
% f表示被积函数；
% a，b表示积分上下限；
% n是子区间的个数
% delta是误差限
M = 1;
h = b - a;
err = 1
J = 0;
R = zeros(4, 4);
R(1, 1) = h * (feval('f', a) + feval('f', b)) / 2
while ((err > delta) & (J < n) ) | (J < 4)
J = J + 1;
h = h / 2;
s = 0;
for p = 1 : M
x = a + h * (2*p - 1);
s = s + feval('f', x);
end
R(J + 1, 1) = R(J, 1) / 2 + h*s;
M = 2*M;
for K = 1 : J
R(J + 1,K + 1) = R(J + 1, K) + (R(J + 1,K) - R(J, K)) / (4^K - 1);
end
err = abs(R(J, J) - R(J + 1,K + 1));
end
quad = R(J + 1, J + 1)；
5 高斯－勒让德求积公式及Matlab实现
5.1 高斯－勒让德求积公式的基本原理
在高斯求积公式
（5.1）
中，若取权函数=1，区间[-1,1]，则得公式
（5.2）
我们知道勒让德多项式是区间[-1,1]上的正交多项式，因此，勒让德多项式的零点就是求积公式(5.2)的高斯点，形如(5.1)式的高斯公式特别的称为高斯-勒让德求积公式．如下表5.1所示为高斯-勒让德求积公式的节点数和系数．
5.2 高斯－勒让德求积公式的Matlab实现
function quad = gauss8(f ,a, b, x, A)
N = length(x);
T = zeros(1, N);
T = (a + b) / 2 + ((b - a) / 2) * x;
quad = ((b - a) / 2) * sum(A. * feval('f', T));

表5.1 高斯-勒让德求积公式的节点数和系数
n

n

0
0.0000000
2.0000000
3
0.8611363
0.3399810
0.3478548
0.6521452
1
0.5773503
1.0000000
4
0.9061978
0.5384693
0.0000000
0.2369269
0.4786287
0.5688889
2
0.7745967
0.0000000
0.5555556
0.8888889
5
0.9324695
0.6612094
0.2386192
0.1713245
0.3607616
0.4679139

6 各个求积公式的应用举例与比较分析
6.1 简单数值积分的解(精确值0.946083070367183)
6.1.1 牛顿－柯特斯当n=1时的梯形算法和n=2时的辛普森算法的结果
解：先用M文件定义一个f1.m的函数
function f = f1(i, x)
g(1) = sqrt(x);
if x == 0
g(2) = 1;
else
g(2) = sin(x) / x;
end
f = g(i);
输入>>Ncotes(0, 1, 1, 2)
回车得到；
ans =
0.9270
输入>>Ncotes（0, 1, 2, 2）
回车得到：
ans =
0.9461
6.1.2 复合梯形公式和复合辛普森求积公式的计算结果
解：建立一个M文件定义一个f.m函数
function y = f(x)
if x == 0
y = 1;
else
y = sin(x) / x
end
输入>>trapr1(‘f’, 0, 1, 10)
回车得到：
ans =
0.9458
输入>>simpr1(‘f’, 0, 1, 10)
ans =
0.9461
6.1.3 龙贝格求积公式的计算结果（取误差不超过）
解：建立一个M文件定义一个f.m函数
function y = f(x)
if x == 0
y = 1;
else
y = sin(x) / x
end
输入>>romber（‘f’, 0, 1, 10，0.5*10^(- 10)）
回车得到：
quad =
0.9461
ans =
0.9207 0 0 0 0
0.9398 0.9461 0 0 0
0.9445 0.9461 0.9461 0 0
0.9457 0.9461 0.9461 0.9461 0
0.9460 0.9461 0.9461 0.9461 0.9461
6.1.4 高斯－勒让德求积公式的计算结果(给定节点3)
解：建立一个M文件定义一个f.m函数
function y = f(x) ；
y = sin(x) / x;
输入
>>gauss8（‘f’,0,1,[-0.7745966692,0.7745966692,0],[0.5555555556,0.5555555556,0.8888888888]）
回车得到：
ans =
0.8956
6.2 复杂数值积分的解(精确值-1.8785)
6.2.1 牛顿－柯特斯当n = 1时的梯形算法和n = 2时的辛普森算法的结果
解：先用M文件定义一个f1.m的函数
function f = f1(i, x)
g(1) = sqrt(x);
if x == 0
g(2) = 1;
else
g(2) = cos(x) – 1 / (1 + x.^2) – 1 / (4.*sqrt(4 + x.^2));
end
f = g(i);
输入>>Ncotes(0, 2*pi, 1, 2)
回车得到
ans =
-1.8692
输入>>Ncotes(0, 2*pi, 2, 2)
回车得到
ans =
-1.8704
6.2.2 复合梯形公式和复合辛普森求积公式的计算结果
解：先用M文件定义一个f1.m的函数
function y = f(x)
if x == 0
y = 1;
else
y = cos(x) – 1 / (1 + x.^2) – 1 / (4.*sqrt(4 + x.^2));
end
输入>>trapr1(‘f’, 0, 2*pi, 10)
回车得到：
ans =
-1.8748
输入>>simpr1(‘f’, 0, 2*pi, 10)
ans =
-1.8694
6.2.3 龙贝格求积公式的计算结果（取误差不超过）
解：先用M文件定义一个f1.m的函数
function y = f(x)
if x == 0
y = 1;
else
y = cos(x) – 1 / (1 + x.^2) – 1 / (4.*sqrt(4 + x.^2));
end
输入>>romber（‘f’, 0, 2*pi ,10，0.5*10^(-10)）
回车得到：
quad =
-1.8764
6.3 各个求积公式的比较分析（以的各个积分结果为例）
表6.1各个数值积分的比较
积分方法
牛顿柯特斯（梯形）
牛顿-柯特斯（辛普森）
复合梯形
复合辛普森
龙贝格
高斯-勒让德
精确值
0.9461
0.9461
0.9461
0.9461
0.9461
0.9461
实际值
0.9270
0.9460
0.9458
0.9460
0.9460
0.8956
误差
0.189
0.0001
0.0003
0.0001
0.0001
0.0505
代数精度
1
2
3
3
10
5

牛顿-柯特斯方法是一种利用插值多项式来构造数值积分的常用方法，这其中梯形
积分方法的误差最大，近似效果最差，辛普森方法的精度比梯形积分高了一个数量级，它的代数精度比梯形积分的代数精度高，能更好地近似积分值；牛顿-柯特斯积分方法的误差比辛普森积分精度高两个数量级．复合梯形积分方法比单独的梯形积分精度高，龙贝格方法收敛速度快、计算精度较高，但是计算量较大．高斯求积方法积分精度高、数值稳定、收敛速度较快，但是节点与系数的计算较麻烦，而且要求已知积分函数．
结 论
本文主要讨论了数值积分的计算方法并通过MATLAB软件编程实现，通过前面的研究我们知道求数值积分近似值的计算方法很多，包括牛顿-柯特斯求积公式、复合求积公式、龙贝格求积公式、高斯求积公式等等．其中牛顿-柯特斯方法是一种利用插值多项式来构造数值积分的常用方法，这其中梯形积分方法的误差最大，近似效果最差，辛普森方法的精度比梯形积分高了一个数量级，它的代数精度比梯形积分的代数精度高，能更好地近似积分值；牛顿-柯特斯积分方法的误差比辛普森积分精度高两个数量级．因此，一般情况下，代数精度越高，积分公式计算精度也越高．但是高阶的牛顿-柯特斯方法的收敛性没有保证，因此，在实际计算中很少使用高阶的牛顿-柯特斯公式．复合梯形积分方法比单独的梯形积分精度高，它的积分精度和被积函数有关，还和复合积分时的步长有关．复合辛普森积分公式比单独的辛普森积分公式高近7个数量级，效果明显．龙贝格方法收敛速度快、计算精度较高，但是计算量较大．高斯求积方法积分精度高、数值稳定、收敛速度较快，但是节点与系数的计算较麻烦、而且要求已知积分函数．一般来说，牛顿-柯特斯方法的代数精度越高，数值积分的效果越好、越精确．当积分区间比较大的时候，可以采用复化积分方法可以得到更好的效果；变步长积分方法不仅可以很好地控制计算误差，并且可以寻找到适当的积分步长；龙贝格积分方法可以更好地利用变步长复合积分公式得到的积分序列从而得到更为精确的数值结果，是一个较好的数值积分方法．高斯求积方法精确度高，收敛性快也是一种很优秀的数值积分方法．
参考文献
[1] 张德丰. Matlab数值分析与应用[M]. 北京: 国防工业出版社, 2007. 1-1
[2] 胡良剑, 孙晓君.MATLAB数学实验[M].  北京: 高等教育出版社, 2006
[3] 李庆扬, 王能超, 易大义. 数值分析[M]. 北京: 清华大学出版社, 2008