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

天津职业技术师大学
课程设计任务书
____________ 理_________ 学院___________ 数学1403 __________ 班学生___________ 课程设计课题：
用数值积分法计算正弦积分函数和余弦积分函数
一、课程设计工作日自2016 年_1_—月_4________ 日至2016 年7
月5 日
二、同组学生：无_______________________________________________________
、课程设计任务要求（包括课题来源、类型、目的和意义、基本要求、完成时间、主要参考资料等）：
课题来源：教师自拟
类型：理论研究
目的和意义：培养学生对数值分析中主要算法的应用能力，探索相关算法之
间的在联系。

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

完成时间：
参考资料：《数值分析》庆扬等清华大学
教研室主任签字:
指导教师签字:
一、问题叙述
用数值积分法计算正弦积分函数和余弦积分函数
提示：正弦积分，余弦si(x)函数ci(x) x嶼dt
0 t t
要求：(1)编写函数，对任意给定的X,可求出值。

(2)使用尽可能少的正余弦函数的调用，计算更精确的值。

(用多种方法，创新方法)
二、问题分析
1、数值积分基本原理：用数值分析求解积分的数值方法有很多，
如简单的梯形法、矩形法、辛普森( Simpson)法、牛顿-科斯特(Newton-Cotes)法等都是常用的方法。

它们的基本思想都是将整个积分区间［a , b］分成n个子区间［x i,x i+i］ , i=1 , 2,…，n,其中X i=a, X n+i二b。

这样求定积分问题就分解为求和问题。

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

三、数值积分法求解问题
1、梯形公式、矩形公式
首先，积分中值定理告诉我们，在积分区间［a, b］存在一点，成立
b f (x)dx= (b-a) f (),就是说，底为b-a而高为f ()的矩形a
面积恰等于所求区边梯形的面积。

如果我们用两端点“高度” f (a) 与f (b)的算术平均值作为平均高度f ()的近似值，这样导出的求积公式^f (x)dx 〜号卩(a) +f (b)］便是我们熟悉的梯形公式。

将积分区间[a,b]n等分，每个小区间宽度均为山=出),h称为积布
n
步长。

记a=x o< x i v…v X k…v x o=b，在小区间上用小矩形面积近似小曲边梯形的面积，若分别取左端点和右端点的函数值为小矩形的高贝S分别得到两个曲边梯形面积的近似计算公式。

具体程序如下：
clear
x=li nspace(0,pi);
dx=x (2);
y=s in (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
2
1.8
1.6
1.4
1.2
1
0.8
0.6
0.4
0.2
0 0.5 1 1.5 2 2.5 3 3.5 由图可知这种方法精度太低，应选择其他方法。

2、quad 函数、quan1函数 正弦：function y=si(t) a=1e-8;
y=quad('si n( x)./x',a,t) y=quadl('si n( x)./x',a,t) 余弦：function y=ci(t) a=-1e1;
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( '辛普森')
x=1:100; for i=1:100
%函数在0点无界，去掉0点
%函数在0点无界，去掉0点
1.9 1
c
1
1
1.8
J
-
1.7 ■ 1
1.6 J 1
1 i
\ 1 1 /
1 / 1 1 J
\J\T\
1.5 1
1: ¥
1.4 —
—
1.3 =
-
1.2 一
-
1.1
一
-
1 --
■-
1
r
1
1
r r \ r
1
辛普森
10
30
40 50 60 70 80 90 100
20
y2(i)=ci(x(i));
end plot(x,y2, 'b') title( '辛
普森')
1400 1200
1000
800 600
400 200
0 -200 -400 辛普森
0 10 20 30 40 50 60 70 80 90 100
给定任意X值，均可计算出对应的正弦、余弦函数积分。

但从结果可以看出精度不是很高。

3、复合求积公式
由于牛顿-科特斯公式在n》8时不具有稳定性，故不可能通过提高阶的方法来提高求积精度。

为了提高精度通常可把积分区间分成若干子区间(通常是等分)，再在每个子区间上用低级求积公式。

这种方法为复合求积法。

331 复合梯形公式
将区间a,b 划分为n 等分，分点x k a kh,h -_ , k 0,1, n
区间
x k
, x k i
k 0,1, ,n 1,
上米用梯形公式，则得
称为复合梯形公式。

复合梯形公式的余项
事实上只要设f X Ca,b ，则可得收敛性，只要把
n 1
n
十 1 b a 「
b a ‘
T n
f X k
f X k 2 n k 0
n k 1
程序如下： 正弦：
fun cti on T_n=fhtxs(a,b, n) h=(b-a)/n; for k=0: n
f(x)dx
X k x k
f(x)dx
h n 1
-
f(x
k
)
f(x k 1
)
2 k 0
R n (f )
T n
1
f X k
X k 1 X k f x -
,n,在每个子
R n f I
T n
12
k x k x k 1
由于 f(x) C 2 a,b ,
min f
0 k n 1
丄n1
k _
n k 0
max ,
0 k n 1
所以 a,b 使
于是复合梯形公式的余项为
R n f
b 12
a h 2
T n 改写成为
x(k+1)=a+k*h;
if x(k+1)==0
x(k+1)=10A(-10);
end
end
T_1= h/2*(SS(x(1))+SS(x( n+1)));
for i=2: n
F(i)=h*SS(x(i));
end
T_2=sum(F);
T_n二T_1+T_2;
余弦：
fun cti on T_n=fhtxc(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)=10A(-l0);
end
end
T_1=h/2*(CC(x(1))+CC(x( n+1)));
for i=2: n
F(i)=h*CC(x(i));
end
T_2=sum(F);
T_n=T_1+T_2;
图像：
正弦余弦
332复合新普斯求积公式
h n 1
6k0
[f(x k ) 4f(x
k12) f(x k )] R n (f)
-
称为复合辛普森求积公式。

程序如下：
正弦
fun cti on S_n=fhxpss(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)=10A (-10); x_k(k+1)=10A (-10); end end
S_1=h/6*(SS(x(1))+SS(x( n+1))); for i=2: n
F_1(i)=h/3*SS(x(i)); end for j=1: n
F_2(j)=2*h/3*SS(x_k(j)); end
S_2=sum(F_1)+sum(F_2); S_n=S_1+S_2;
余弦:
fun cti on S_n=fhxpsc(a,b ,n)
h=(b-a)/n; for k=0:n x(k+1)=a+k*h;
将区间［a,b ］划分为n 等分，
在每个子区间X k ，X ki 上采用辛普森公式，若记
1 X k1
2 X k 2
h,
则得
f(x)dx
f(x)dx
k 0
x_k(k+1)=x(k+1)+1/2*h;
if (x(k+1)==0)||(x_k(k+1)==0)
x(k+1)=10A(-10);
x_k(k+1)=10A(-10);
end
end
S_1=h/6*(CC(x(1))+CC(x( n+1)));
for i=2: n
F_1(i)=h/3*CC(x(i));
end
for j=1: n
F_2(j)=2*h/3*CC(x_k(j));
end
S_2=sum(F_1)+sum(F_2);
S_n=S_1+S_2;
图像与复合梯形所得图像基本相同，深入分析两只复合函数的优劣，对于积分函数si(x):平dt假设x=1，则将区间［0,1］划分为8等份, 应用复合梯形求得
T8=0.9456909
而如果将［0，1］分为4等份，应用复合辛普森有
S4=0.9460832
通过参考数值分析（庆阳）的结论，发现无论是复合梯形公式还是复合辛普森公式，最终结果都会随着h值的减小而更加精确。

对复合梯形公式和复合辛普森公式计算出的结果进行比较，发现复合梯形法的结果T8只有两位有效数字，而复合辛普森的结果却有六位有效数字，所以复合辛普森公式计算出的结果更加的精确。

4、插值型的求积公式
clc, clear
x0=0:0.5:5;
y0=［ Inf 1.7552 0.5403 0.0472 -0.2081 -0.3205 -0.3300
-0.2676 -0.1634 -0.0468 0.0567
]; %所求积分函数的数值
pp=csape(x0,y0) ; %默认的边界条件，Lagrange 边界条件
format long g
chazhi=pp.coefs %显示每个区间上三次多项式的系数
s=quadl((t)ppval(pp,t),0,5)%求积
分
format %恢复短小数的显示格式
x=0:0.1:5;
y=cos(x)/x;
y1=spline(x0,y0,x);
z=0*x;
hold on
plot(x,z,x,y, 'k--' ,x,y1, 'r')
plot(x0,y0, '*' )
hold off
clear x0=0:0.5:5; y0=[ Inf 1.7552 0.5403 -0.1634 -0.0468 0.0567 ];
pp=csape(x0,y0) ; format long g chazhi=pp.coefs s=quadl((t)ppval(pp,t),0,5) format x=0:0.1:5;
y=cos(x)/x;
y1=spline(x0,y0,x); z=0*x;
hold on
plot(x,z,x,y,'k--',x,y1,'r') plot(x0,y0,'*') hold off
如图所示：
0.0472 -0.2081 -0.3205 -0.3300
%所求积分函数的数值
%默认的边界条件，Lagrange 边界条件
%显示每个区间上三次多项式的系数
%求积分
%恢复短小数的显示格式
-0.2676
-a ma s; ^ra 53 n -CL
B3[arQCPUn2Ff 0.033653 LL<64.^9£ -fl.a8LW3«fFlllM]
■ OLMQBMAIF 禅 1
附 -A
-C HJDL34I
E
■H 仲12吧骷林肪e?
-LDMOJIiaBlZket ■0机
展斷2i2»<$時】 -
Q.1^B946M«17V9 -
CL10rT943F3l ：W72
fl
朋制対创J 駁
CL 17抄1鮒拥附ID
血f
2加CBIWT D 恥JD 眩 L 5IJ!3rBD!：Sfil f?9' fl-25T9SF94K3TO122
山.JIESMMHS 却調 1
Uui^otdariSTirM CL
CL ： 3 S3] UI 354 2
5、高斯求积公式
fun cti on
[ql,Ak,xk]=gsqj(fu n,a,b, n,tol)
if n arg in==1
a=-1;b=1; n=7;tol=1e-8; elseif nargin==3 n=7;tol=1e-8; elseif nargin==4
tol=1e-8;
elseif nargin==2|| nargin>5
error( 'The Number of In put Argume nts Is Wrong!'
);
end
%计算求积节点 syms x
p=sym2poly(diff((x A 2-1)A ( n+1), n+1))/(2A n*factorial( n)); tk=roots(p); % 求积节点
%计算求积系数 Ak=zeros( n+1,1);
£31« I LJ t. D Q 血14 血却 iLfedsr jfa 血
|D 上x
■為"--
W b
強 ■: ifH 哄 恋咗■ClVAILiB^MSCh
*■
yirrlpJ! [£jl Hp-fc *1.1 也划 _l|、-M' l b*
iTffl'jrari 1 throuih 2
CB I IJ1F F 1 1hr DIL|}| 4
iL022asr^s meres
・九龙初P 爲
-4L O32»52225^21QL4
-叽OMSK 題開・詐帕】 l.r55£
0.5405 0.04^3 -0. ~ n0| -L 32D6 -0.訶兀
0.1634
CL^WLZIiBa^VMB
for i=1: n+1
xkt=tk;
xkt(i)=[];
pn=poly(xkt);
fp=(x)polyval(p n, x)/polyval(p n,tk(i));
Ak(i)=quadl(fp,-1,1,tol); % 求积系数end
%积分变量代换，将[a,b] 变换到[-1,1]
xk=(b-a)/2*tk+(b+a)/2;
%检验积分函数fun有效性
fun=fcn chk(fu n, 'vectorize' );
%计算变量代换之后积分函数的值
SS=fu n(xk)*(b-a)/2;
%计算积分值
ql=sum(Ak.*SS);
计划表。