数值积分用matlab实现

东北大学秦皇岛分校

数值计算课程设计报告

数值积分及Matlab实现

学院数学与统计学院

专业信息与计算科学

学号*******

姓名楚文玉

指导教师张建波姜玉山

成绩

教师评语：

指导教师签字：

2015年07月14日

1 绪论

在科研计算中，经常会碰到一些很难用公式定理直接求出精确解的积分问题，对于这类问题，我们一般转化为数值积分问题，用计算机来实现求解问题． 1.1 课题的背景

对于定积分()b

a f x dx ⎰在求某函数的定积分时，在一定条件下，虽然有牛顿-莱布里

茨公式()()()b

a

I f x dx F b F a ==-⎰可以计算定积分的值，但在很多情况下的原函数()

f x 不易求出或非常复杂．被积函数的原函数很难用初等函数表达出来，例如

2

sin (),x x f x e x

-=

等；有的函数()f x 的原函数()F x 存在，但其表达式太复杂，计算量太大，有的甚至无法有解析表达式．因此能够借助牛顿-莱布尼兹公式计算定积分的情形是不多的．另外，许多实际问题中的被积函数()f x 往往是列表函数或其他形式的非连续函数，对这类函数的定积分，也不能用不定积分方法求解，只能设法求其近似值．因此，探讨近似计算的数值积分方法是有明显的实际意义的，即有必要研究定积分的数值计算方法，以解决定积分的近似计算．而数值积分就是解决此类问题的一种有效的方法，它的特点是利用被积函数在一些节点上的信息求出定积分的近似值．微积分的发明是人类科学史上一项伟大的成就，在科学技术中，积分是经常遇到的一个重要计算环节数值积分是数学上重要的课题之一，是数值分析中重要的内容之一．随着计算机的出现，近几十年来，对于数值积分问题的研究已经成为一个很活跃的研究领域．现在，数值积分在计算机图形学，积分方程，工程计算，金融数学等应用科学领域都有着相当重要的应用，所以研究数值积分问题有着很重要的意义．国内外众多学者在数值积分应用领域也提出了许多新方法．在很多实际应用中，只能知道积分函数在某些特定点的取值，比如天气测量中的气温、湿度、气压等，医学测量中的血压、浓度等等．通过这个课题的研究，我们将会更好地掌握运用数值积分算法求出特殊积分函数的定积分的一些基本方法、理论基础；并且通过Matlab 软件编程的实现，应用于实际生活中． 1.2 课题的主要内容框架 1.2.1 数值积分各求积公式简介

简介牛顿-柯特斯求积公式及其辛普森求积公式，龙贝格求积公式，高斯求积公式的基本理论基础和方法． 1.2.2 求积公式的代码实现

通过理解各种数值积分求积公式的原理方法，通过Matlab 软件编程，实现以上求积公式． 1.2.3 应用举例

通过简单举例，自建一个相对简单和复杂的函数，用上面编写的Matlab 源程序来解决实际问题，体会数值积分和Matlab 的优势．

2 牛顿－柯特斯公式及Matlab 实现

2.1 牛顿－柯特斯公式的基本原理方法

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

h n

-=，选取等距节点k x a kh =+构造出的差值型求积公式

()0()()n

n n k k k I b a C f x ==-∑， （2.1）

称为牛顿-柯特斯公式，式中()n k C 称为柯特斯系数．根据

()b

k a A l x dx =⎰， 0,1,2,

,.k n = （2.2）

引进变量代换x a th =+，则有

()0000

(1)()!()!n k

n n n n n k

j j j k

j k

h t j C

dt t j dt b a k j nk n k -==≠≠--==----∏∏⎰⎰ （2.3） 当n = 2时，此时柯特斯系数为(2)

16C =，(2)146C =，(2)

2

16

C =，相应的求积公式就是辛普森求积公式：

()4()()62b a a b S f a f f b -+⎡⎤=

++⎢⎥⎣⎦

（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 等分，分点k x a kh =+，,0,1,

b a

k n n

-=,在每个子区间[1,k k x x +]

（0,1,1k n =-）上采用梯形公式得： 1

1

1

10

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

n n b

x k k n a

x k k h I f x dx f x dx f x f x R f +--+=====++∑∑⎰⎰

（3.1）

记

11

100

[()()][()2()()]22n n n k k k k k h h

T f x f x f a f x f b --+===+=++∑∑ （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);