数值分析法求正弦余弦积分函数

天津职业技术师范大学

课程设计任务书

理学院数学1403 班学生张群

课程设计课题：

用数值积分法计算正弦积分函数和余弦积分函数

一、课程设计工作日自 2016 年 7 月 4 日至 2016 年 7 月 5日

二、同组学生：无

三、课程设计任务要求（包括课题来源、类型、目的和意义、基本要求、完成时

间、主要参考资料等）：

课题来源：教师自拟

类型：理论研究

目的和意义：培养学生对数值分析中主要算法的应用能力，探索相关算法之间的内在联系。

基本要求：根据数值分析课程所学的知识，应用Matlab软件编写程序，完成对算法及其内在原理的实验研究。

完成时间：

参考资料：《数值分析》李庆扬等清华大学出版社

指导教师签字：教研室主任签字：

一、问题叙述

用数值积分法计算正弦积分函数和余弦积分函数

提示：正弦积分，余弦0sin ()x

t si x dt t =⎰函数cos ()x

t ci x dt t

-∞=⎰

要求：（1）编写函数，对任意给定的x ，可求出值。

（2）使用尽可能少的正余弦函数的调用，计算更精确的值。（用多种方法，创新方法）

二、问题分析

1 、数值积分基本原理：用数值分析求解积分的数值方法有很多，如简单的梯形法、矩形法、辛普森（Simpson ）法、牛顿-科斯特（Newton-Cotes ）法等都是常用的方法。它们的基本思想都是将整个积分区间[a ，b]分成n 个子区间[x i ,x i+1]，i=1，2，…，n ，其中x 1=a ，x n+1=b 。这样求定积分问题就分解为求和问题。

2、本题要求用数值积分法计算正弦积分函数和余弦函数积分，那么应该从编写函数的入手，建立function 的m 文件，通过对函数的调用，对任意跟定的x 值，求出积分函数的值。 三、数值积分法求解问题 1、 梯形公式、矩形公式

首先，积分中值定理告诉我们，在积分区间[a ，b]内存在一点ξ，成立⎰b

a x f ）（dx=（b-a ）f （ζ），就是说，底为b-a 而高为f （ζ）的矩形面积恰等于所求区边梯形的面积。如果我们用两端点“高度”f （a ）与f （

b ）的算术平均值作为平均高度f （ξ）的近似值，这样导出的求积公式⎰b

a x f ）（dx ≈

2

a

-b [f （a ）+f （b ）]便是我们熟悉的梯形公式。

将积分区间[a,b]n 等分，每个小区间宽度均为h=

n

a -

b ）

（，h 称为积布步长。记a=x 0＜x 1＜…＜x k …＜x o =b ，在小区间上用小矩形面积近似小曲边梯形的面积,若分别取左端点和右端点的函数值为小矩形的高,则分别得到两个曲边梯形面积的近似计算公式。具体程序如下：

clear x=linspace(0,pi); dx=x(2); y=sin(x); s1=sum(y)*dx s2=trapz(y)*dx

sc1=cumsum(y)*dx; sc2=cumtrapz(y)*dx;

plot(x,-cos(x)+1,x,sc1,'.',x,sc2,'o')

hold on

0.5

1

1.5

2

2.5

3

3.5

00.20.40.60.811.21.41.61.82

由图可知这种方法精度太低，应选择其他方法。

2、quad 函数、quan1函数

正弦：function y=si(t)

a=1e-8; %函数在0点无界，去掉0点 y=quad('sin(x)./x',a,t) y=quadl('sin(x)./x',a,t) 余弦：function y=ci(t)

a=-1e1; %函数在0点无界，去掉0点 y=quad('cos(x)./x',a,t) y=quadl('cos(x)./x',a,t) 图像:

x=1:100; for i=1:100 y2(i)=si(x(i)); end

plot(x,y2,'r') title('辛普森')

10

20

30

40

50

60

70

80

90

100

0.911.11.21.31.41.51.61.71.8

1.9辛普森

x=1:100; for i=1:100

y2(i)=ci(x(i)); end

plot(x,y2,'b') title('辛普森')

0102030405060708090100

-400

-200020040060080010001200

1400辛普森

给定任意x 值，均可计算出对应的正弦、余弦函数积分。但从结果可以看出精度不是很高。 3、复合求积公式

由于牛顿-科特斯公式在n ≥8时不具有稳定性，故不可能通过提高阶的方法来提高求积精度。为了提高精度通常可把积分区间分成若干子区间（通常是等分），再在每个子区间上用低级求积公式。这种方法为复合求积法。

3.3.1 复合梯形公式

将区间[]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

1f 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 ，

称为复合梯形公式。 复合梯形公式的余项

()()()

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 ηη 所以()b a ,∈∃η使 ()()k n k f n f ηη∑-==10

'

''

'1

于是复合梯形公式的余项为

()()η'

'212

f h a b f R n --

=

事实上只要设()[]b a C x f ,∈，则可得收敛性，只要把n T 改写成为

()()]∑∑=-=-+⎢⎣⎡-=n

k k n k k n x f n a b x f n a b T 1

1021 程序如下： 正弦：

function T_n=fhtxs(a,b,n)

h=(b-a)/n; for k=0:n