数值分析法求正弦余弦积分函数.docx
《数值分析法求正弦余弦积分函数.docx》由会员分享，可在线阅读，更多相关《数值分析法求正弦余弦积分函数.docx（16页珍藏版）》请在冰豆网上搜索。
数值分析法求正弦余弦积分函数
天津职业技术师范大学
课程设计 任 务 书
理 学院 数学1４03 班 学生张群
课程设计课题:
用数值积分法计算正弦积分函数和余弦积分函数
一、课程设计工作日自 ２01６ 年 7 月4 日至 ２016年 ７月 5日
二、同组学生：
无
三、课程设计任务要求（包括课题来源、类型、目的和意义、基本要求、完成时间、主要参考资料等）：
课题来源：
教师自拟
类型：
理论研究
目的和意义:
培养学生对数值分析中主要算法的应用能力,探索相关算法之间的内在联系。
基本要求：
根据数值分析课程所学的知识，应用Matｌab软件编写程序,完成对算法及其内在原理的实验研究。
完成时间:
参考资料：
《数值分析》李庆扬等清华大学出版社
指导教师签字：
教研室主任签字:
一、问题叙述
用数值积分法计算正弦积分函数和余弦积分函数
提示：
正弦积分,余弦
函数
要求:
（1）编写函数，对任意给定的ｘ，可求出值.
（2）使用尽可能少的正余弦函数的调用,计算更精确的值.（用多种方法,创新方法）
2、问题分析
1、数值积分基本原理：
用数值分析求解积分的数值方法有很多,如简单的梯形法、矩形法、辛普森（Sｉmｐsoｎ）法、牛顿-科斯特（Newton—Cotes）法等都是常用的方法。
它们的基本思想都是将整个积分区间［a，ｂ]分成n个子区间[xｉ,xｉ＋１],i=1,2，…，n，其中x１=a,xｎ+1＝b。
这样求定积分问题就分解为求和问题。
2、本题要求用数值积分法计算正弦积分函数和余弦函数积分，那么应该从编写函数的入手，建立fｕnction的m文件，通过对函数的调用，对任意跟定的x值，求出积分函数的值。
3、数值积分法求解问题
1、 梯形公式、矩形公式
首先,积分中值定理告诉我们，在积分区间［a，b］内存在一点
成立
dx=（b—ａ）ｆ（
），就是说,底为b－a而高为f（
）的矩形面积恰等于所求区边梯形的面积。
如果我们用两端点“高度”ｆ（a）与f（b）的算术平均值作为平均高度f（
）的近似值,这样导出的求积公式
dｘ≈
［f（a）＋f（ｂ）]便是我们熟悉的梯形公式。
将积分区间[a,b]n等分,每个小区间宽度均为h=
，h称为积布步长.记a=x０＜x1＜…＜xk…＜ｘｏ=b,在小区间上用小矩形面积近似小曲边梯形的面积，若分别取左端点和右端点的函数值为小矩形的高，则分别得到两个曲边梯形面积的近似计算公式.具体程序如下：
cｌear
x＝lｉｎｓpaｃe（０,pi）；
ｄx＝ｘ
（2）;
y=sin（x）;
s1=sｕｍ（y）＊ｄx
s2＝ｔrapｚ（ｙ）＊dx
sｃ1=ｃumsuｍ（ｙ）＊dx;
sc2＝ｃｕmｔraｐz（y）*dx;
pｌoｔ（x,—cｏs（x）+１，x,sc1，＇。
'，x，sc2,'o＇）
ｈoldon
由图可知这种方法精度太低，应选择其他方法.
2、quad函数、quaｎ1函数
正弦：
fｕnctiｏn ｙ=si（t）
a=1e—8; ％函数在0点无界，去掉０点
y=ｑuad（'sin（ｘ）．／x’,a,ｔ）
y=ｑuadl（＇sin（x）．/x'，a，t）
余弦：
fuｎctｉｏｎ y=ci（t）
a=-１e1； ％函数在0点无界，去掉0点
y＝quad（＇cos（x）．/x＇,ａ，t）
y=ｑuadl（＇ｃoｓ（x）。
/x'，a,ｔ）
图像：
ｘ＝1:
1００；
fori＝1：
10０
y2（i）=ｓi（x（i））；
end
ploｔ（x，y2，＇r’）
titｌe（’辛普森'）
x=1：
100；
for i=1:
100
y２（i）=cｉ（x（i））;
eｎd
plot（x,y２,’b’）
titｌe（'辛普森＇）
给定任意x值,均可计算出对应的正弦、余弦函数积分.但从结果可以看出精度不是很高.
3、复合求积公式
由于牛顿—科特斯公式在n≥8时不具有稳定性，故不可能通过提高阶的方法来提高求积精度。
为了提高精度通常可把积分区间分成若干子区间（通常是等分），再在每个子区间上用低级求积公式。
这种方法为复合求积法.
3。
3.１ 复合梯形公式
将区间
划分为n等分，分点
在每个子区间
上采用梯形公式，则得
记
称为复合梯形公式。
复合梯形公式的余项
由于
ﻩ且
所以
使
于是复合梯形公式的余项为
事实上只要设
，则可得收敛性,只要把
改写成为
程序如下:
正弦:
functioｎT_ｎ=ｆhtｘs（a,b,n）
h=（b—ａ）／ｎ；
ｆork=0：
n
x（ｋ+1）=ａ＋k＊h;
ifx（k+1）=＝0
x（k+１）=10^（-10）；
end
eｎd
Ｔ_1＝h/2*（SS（x
（1））+SS（ｘ（n＋1）））；
ｆoｒi=2：
ｎ
F（i）=h*SS（x（i））;
end
T_2＝sum（F）；
T_n=T_1+Ｔ＿２;
余弦：
ｆuｎcｔiｏｎT_n=fhtxc（ａ,b，n）
ｈ=（b-ａ）/ｎ；
ｆork=0:
n
x（k＋1）=ａ＋ｋ*h；
ifx（k+1）==0
x（k＋１）=1０＾（-10）；
eｎd
ｅnd
T_1＝ｈ／２＊（ＣＣ（x
（1））+CC（x（n+1）））；
foｒi=2：
n
F（i）＝h*CＣ（x（i））;
end
T＿２=suｍ（F）；
T＿n=Ｔ_１+T＿2;
图像:
正弦 余弦
3.3．2复合新普斯求积公式
将区间
划分为n等分，在每个子区间
上采用辛普森公式，若记
则得
●
称为复合辛普森求积公式。
程序如下：
正弦
ｆunctiｏnS＿n=fhxｐｓｓ（a，b，n）
h=（b－ａ）／n;
ｆｏｒk=0:
ｎ
ｘ（k+１）=a+ｋ＊ｈ;
x_k（k+1）=ｘ（k+1）＋1/2＊h;
iｆ（x（ｋ+1）＝=0）|｜（x_k（k+１）=＝０）
x（k+1）＝1０^（—10）；
x_ｋ（k+1）＝1０＾（-１0）；
end
ｅnd
S＿1＝h/6*（ＳS（x（１））+ＳS（ｘ（ｎ+１）））；
fori=2：
n
F_1（i）=h/３*SS（x（i））；
enｄ
fｏrj=1:
n
F_2（j）＝2*h／3*SS（x＿ｋ（ｊ））;
ｅnd
Ｓ_２=sｕｍ（Ｆ_1）+sum（Ｆ_2）；
S_ｎ=S＿1+S＿2;
余弦：
fｕｎctionS_n＝fhｘpsｃ（ａ,ｂ,n）
ｈ＝（b-a）/n;
fork=0：
n
x（ｋ＋１）=a+ｋ＊h；
x＿ｋ（ｋ＋1）＝x（k＋1）+1/２＊ｈ；
if（ｘ（k＋１）=＝0）｜|（x_ｋ（ｋ+１）＝=0）
ｘ（k＋1）=10^（-10）;
ｘ_k（k＋1）=１0^（-10）;
ｅnｄ
ｅnd
S＿1＝h／6＊（CC（ｘ
（1））+CC（x（n+1）））；
fori=２:
n
F_1（ｉ）=h/3*CC（x（i））；
ｅnｄ
ｆorj=１：
ｎ
F_2（j）＝2＊h／3*CＣ（x_k（j））;
ｅnd
S_2=sum（F＿1）+sｕm（F＿２）；
Ｓ_n=S_１+S_2；
图像与复合梯形所得图像基本相同，深入分析两只复合函数的优劣，
对于积分函数
假设ｘ=1,则将区间[0，１］划分为8等份，应用复合梯形求得
T8＝0。
94569０９
而如果将[0，1］分为4等份，应用复合辛普森有
S4＝0.9460８32
通过参考数值分析（李庆阳）的结论,发现无论是复合梯形公式还是复合辛普森公式,最终结果都会随着h值的减小而更加精确.对复合梯形公式和复合辛普森公式计算出的结果进行比较，发现复合梯形法的结果Ｔ８只有两位有效数字，而复合辛普森的结果却有六位有效数字,所以复合辛普森公式计算出的结果更加的精确。
4、插值型的求积公式
clｃ， cｌear
ｘ0=0:
0．5：
5；
ｙ０=［ Ｉnｆ 1.７５５２ 0。
5403 0。
0472－0.20８1—0。
３20５ -0。
3３00 －0。
２67６ -0。
16３4 －0.０468 0。
0567
]；％所求积分函数的数值
pp=ｃsaｐe（x0,ｙ0） ;％默认的边界条件,Ｌａgｒangｅ边界条件
ｆorｍat ｌｏｎｇg
cｈazｈi＝ｐp。
coefｓ ％显示每个区间上三次多项式的系数
s=quaｄl（@（t）ppvａl（pp，t），０，5） %求积分
foｒmat ％恢复短小数的显示格式
x=0:
0。
1:
5；
ｙ=coｓ（x）／x；
y1=splｉne（ｘ0，y０，x）;
z=0＊ｘ；
hold on
ｐlot（x,z，ｘ，y，＇k—-’，x,y1，'ｒ'）
plot（x0，y０，’＊’）
holdoff
ｃｌeaｒ
x0=0：
0.5：
５;
ｙ0=［Ｉｎf 1．7552 0.5403 0．04７2—０。
2081 -０.3205 —0.3３00 —0．2６７6 -0.1634 -0。
04680.05６7
］； %所求积分函数的数值
pｐ=ｃsape（x0,y0）； %默认的边界条件，Lagrangｅ边界条件
ｆoｒmatｌoｎｇｇ
ｃｈazhi=pp。
coefs ％显示每个区间上三次多项式的系数
ｓ＝quａdｌ（@（ｔ）ppvａl（pp,t）,0,5） ％求积分
foｒｍat ％恢复短小数的显示格式
x＝0:
0。
1：
５；
ｙ=cos（ｘ）/x;
y1=ｓｐline（ｘ0,y0，x）;
ｚ＝0＊x；
hold ｏn
plｏt（x,z，x，ｙ，’k——’，x,y1，’r'）
pｌot（x0，y0，'*'）
holｄｏｆf
如图所示:
5、高斯求积公式
funｃtｉｏn[ql，Aｋ,xk］=gｓｑj（ｆun,a,b，n，tol）
ifnａrｇｉn＝＝１
a=—1；b=１；n=7;tｏｌ=1ｅ-8；
elseｉfnａrgiｎ==３
n＝7;tol=1e-8；
elｓｅifnargin==４
tol=1e－8;
ｅlseiｆnarｇin==2｜｜naｒgin>5
erｒor（＇ＴｈｅNumbｅrofInpｕｔAｒgumｅntsIsWrｏnｇ！
＇）；
end
%计算求积节点
symsx
p＝sym２ｐoｌy（difｆ（（ｘ^２-1）^（n+1）,n＋1））/（2^ｎ*ｆａctｏｒｉａｌ（ｎ））;
tk=rｏots（p）；％求积节点
％计算求积系数
Aｋ＝ｚｅros（ｎ+1,1）;
fｏｒ i=1:
ｎ+1
xkｔ=tk;
xkt（i）=［］;
pｎ=poly（xkｔ）；
fp=＠（x）ｐoｌyvａl（pn，ｘ）/ｐolyvaｌ（ｐn，tｋ（i））；
Aｋ（i）=qｕａdl（ｆp,—1，1,tｏl）;%求积系数
end
％ 积分变量代换，将［a，ｂ]变换到[-1,1]
xk=（b－a）/2＊tｋ+（ｂ+a）/2;
%检验积分函数fun有效性
ｆun=ｆｃｎｃhｋ（fun,'ｖectorｉze’）；
％ 计算变量代换之后积分函数的值
ＳS=fun（xｋ）＊（b-a）/２；
％ 计算积分值
ｑｌ=sum（Ak．＊SＳ）；
计划表
7月3号
熟悉问题，准备工作，借阅相应的资料,搞清楚题目的用意
题目要求多种方法计算，并尽量减少函数的调用。
7.4
归纳总结多种数值求积的方法,找到各种方法对应的ｍａtlab程序。
梯形
②辛普森公式
③复化辛普森、复化梯形④高斯勒让德
⑤插值
7。
5
编写程序,对所找的方法逐一处理编程，思考跟好的方法.
运行编程结果，进行检查改进
7。
６
编写报告总结,对程序进行系统性总结,完成课程设计
7。
７
修改论文，修改程序,检查修正出现的错误
升级会员 