数值分析法求正弦余弦积分函数之欧阳化创编.docx
《数值分析法求正弦余弦积分函数之欧阳化创编.docx》由会员分享，可在线阅读，更多相关《数值分析法求正弦余弦积分函数之欧阳化创编.docx（11页珍藏版）》请在冰豆网上搜索。
数值分析法求正弦余弦积分函数之欧阳化创编
天津职业技术师范大学
时间：
2021.02.12
创作人：
欧阳化
课程设计任务书
理学院数学1403班学生张群
课程设计课题：
用数值积分法计算正弦积分函数和余弦积分函数
一、课程设计工作日自2016年7月4日至2016年7月5日
二、同组学生：
无
三、课程设计任务要求（包括课题来源、类型、目的和意义、基本要求、完成时间、主要参考资料等）：
课题来源：
教师自拟
类型：
理论研究
目的和意义：
培养学生对数值分析中主要算法的应用能力，探索相关算法之间的内在联系。
基本要求：
根据数值分析课程所学的知识，应用Matlab软件编写程序，完成对算法及其内在原理的实验研究。
完成时间：
参考资料：
《数值分析》李庆扬等清华大学出版社
指导教师签字：
教研室主任签字：
一、问题叙述
用数值积分法计算正弦积分函数和余弦积分函数
提示：
正弦积分，余弦
函数
要求：
（1）编写函数，对任意给定的x，可求出值。
（2）使用尽可能少的正余弦函数的调用，计算更精确的值。
（用多种方法，创新方法）
二、问题分析
1、数值积分基本原理：
用数值分析求解积分的数值方法有很多，如简单的梯形法、矩形法、辛普森（Simpson）法、牛顿-科斯特（Newton-Cotes）法等都是常用的方法。
它们的基本思想都是将整个积分区间[a，b]分成n个子区间[xi,xi+1]，i=1，2，…，n，其中x1=a，xn+1=b。
这样求定积分问题就分解为求和问题。
2、本题要求用数值积分法计算正弦积分函数和余弦函数积分，那么应该从编写函数的入手，建立function的m文件，通过对函数的调用，对任意跟定的x值，求出积分函数的值。
三、数值积分法求解问题
1、梯形公式、矩形公式
首先，积分中值定理告诉我们，在积分区间[a，b]内存在一点
，成立
dx=（b-a）f（
），就是说，底为b-a而高为f（
）的矩形面积恰等于所求区边梯形的面积。
如果我们用两端点“高度”f（a）与f（b）的算术平均值作为平均高度f（
）的近似值，这样导出的求积公式
dx≈
[f（a）+f（b）]便是我们熟悉的梯形公式。
将积分区间[a,b]n等分，每个小区间宽度均为h=
，h称为积布步长。
记a=x0＜x1＜…＜xk…＜xo=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'）
holdon
由图可知这种方法精度太低，应选择其他方法。
2、quad函数、quan1函数
正弦：
functiony=si（t）
a=1e-8;%函数在0点无界，去掉0点
y=quad（'sin（x）./x',a,t）
y=quadl（'sin（x）./x',a,t）
余弦：
functiony=ci（t）
a=-1e1;%函数在0点无界，去掉0点
y=quad（'cos（x）./x',a,t）
y=quadl（'cos（x）./x',a,t）
图像:
x=1:
100;
fori=1:
100
y2（i）=si（x（i））;
end
plot（x,y2,'r'）
title（'辛普森'）
x=1:
100;
fori=1:
100
y2（i）=ci（x（i））;
end
plot（x,y2,'b'）
title（'辛普森'）
给定任意x值，均可计算出对应的正弦、余弦函数积分。
但从结果可以看出精度不是很高。
3、复合求积公式
由于牛顿-科特斯公式在n≥8时不具有稳定性，故不可能通过提高阶的方法来提高求积精度。
为了提高精度通常可把积分区间分成若干子区间（通常是等分），再在每个子区间上用低级求积公式。
这种方法为复合求积法。
3.3.1复合梯形公式
将区间
划分为n等分，分点
在每个子区间
上采用梯形公式，则得
记
，
称为复合梯形公式。
复合梯形公式的余项
由于
且
所以
使
于是复合梯形公式的余项为
事实上只要设
，则可得收敛性，只要把
改写成为
程序如下：
正弦：
functionT_n=fhtxs（a,b,n）
h=（b-a）/n;
fork=0:
n
x（k+1）=a+k*h;
ifx（k+1）==0
x（k+1）=10^（-10）;
end
end
T_1=h/2*（SS（x
（1））+SS（x（n+1）））;
fori=2:
n
F（i）=h*SS（x（i））;
end
T_2=sum（F）;
T_n=T_1+T_2;
余弦：
functionT_n=fhtxc（a,b,n）
h=（b-a）/n;
fork=0:
n
x（k+1）=a+k*h;
ifx（k+1）==0
x（k+1）=10^（-10）;
end
end
T_1=h/2*（CC（x
（1））+CC（x（n+1）））;
fori=2:
n
F（i）=h*CC（x（i））;
end
T_2=sum（F）;
T_n=T_1+T_2;
图像：
正弦余弦
3.3.2复合新普斯求积公式
将区间
划分为n等分，在每个子区间
上采用辛普森公式，若记
则得
称为复合辛普森求积公式。
程序如下：
正弦
functionS_n=fhxpss（a,b,n）
h=（b-a）/n;
fork=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）=10^（-10）;
x_k（k+1）=10^（-10）;
end
end
S_1=h/6*（SS（x
（1））+SS（x（n+1）））;
fori=2:
n
F_1（i）=h/3*SS（x（i））;
end
forj=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;
余弦：
functionS_n=fhxpsc（a,b,n）
h=（b-a）/n;
fork=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）=10^（-10）;
x_k（k+1）=10^（-10）;
end
end
S_1=h/6*（CC（x
（1））+CC（x（n+1）））;
fori=2:
n
F_1（i）=h/3*CC（x（i））;
end
forj=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;
图像与复合梯形所得图像基本相同，深入分析两只复合函数的优劣，
对于积分函数
假设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=[Inf1.75520.54030.0472-0.2081-0.3205-0.3300-0.2676-0.1634-0.04680.0567
];%所求积分函数的数值
pp=csape（x0,y0）;%默认的边界条件，Lagrange边界条件
formatlongg
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;
holdon
plot（x,z,x,y,'k--',x,y1,'r'）
plot（x0,y0,'*'）
holdoff
clear
x0=0:
0.5:
5;
y0=[Inf1.75520.54030.0472-0.2081-0.3205-0.3300-0.2676-0.1634-0.04680.0567
];%所求积分函数的数值
pp=csape（x0,y0）;%默认的边界条件，Lagrange边界条件
formatlongg
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;
holdon
plot（x,z,x,y,'k--',x,y1,'r'）
plot（x0,y0,'*'）
holdoff
如图所示：
5、高斯求积公式
function[ql,Ak,xk]=gsqj（fun,a,b,n,tol）
ifnargin==1
a=-1;b=1;n=7;tol=1e-8;
elseifnargin==3
n=7;tol=1e-8;
elseifnargin==4
tol=1e-8;
elseifnargin==2||nargin>5
error（'TheNumberofInputArgumentsIsWrong!
'）;
end
%计算求积节点
symsx
p=sym2poly（diff（（x^2-1）^（n+1）,n+1））/（2^n*factorial（n））;
tk=roots（p）;%求积节点
%计算求积系数
Ak=zeros（n+1,1）;
fori=1:
n+1
xkt=tk;
xkt（i）=[];
pn=poly（xkt）;
fp=@（x）polyval（pn,x）/polyval（pn,tk（i））;
Ak（i）=quadl（fp,-1,1,tol）;%求积系数
end
%积分变量代换，将[a,b]变换到[-1,1]
xk=（b-a）/2*tk+（b+a）/2;
%检验积分函数fun有效性
fun=fcnchk（fun,'vectorize'）;
%计算变量代换之后积分函数的值
SS=fun（xk）*（b-a）/2;
%计算积分值
ql=sum（Ak.*SS）;
计划表
7月3号
熟悉问题，准备工作，借阅相应的资料，搞清楚题目的用意
题目要求多种方法计算，并尽量减少函数的调用。
7.4
归纳总结多种数值求积的方法，找到各种方法对应的matlab程序。
梯形
辛普森公式
复化辛普森、复化梯形④高斯勒让德
⑤插值
7.5
编写程序，对所找的方法逐一处理编程，思考跟好的方法。
运行编程结果，进行检查改进
7.6
编写报告总结，对程序进行系统性总结，完成课程设计
7.7
修改论文，修改程序，检查修正出现的错误
时间：
2021.02.12
创作人：
欧阳化
升级会员 