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

1、天津职业技术师范大学课程设计任务书理学院数学1403班学生.张群课程设计课题：用数值积分法计算正弦积分函数和余弦积分函数一、课程设计工作日自2016年7月4日至2016年7月5日二、同组学生：无三、课程设计任务要求（包括课题来源、类型、目的和意义、基本要求、完成时问、主要参考资料等）：课题来源：教师自拟类型：理论研究目的和意义：培养学生对数值分析中主要算法的应用能力，探索相关算法之间的内在联系。基本要求：根据数值分析课程所学的知识，应用Matlab软件编写程序，完成对算法及其内在原理的实验研究。完成时间：参考资料：数值分析李庆扬等清华大学出版社指导教师签字:教研室主任签字:一、问题叙述用数值积

2、分法计算正弦积分函数和余弦积分函数提示：正弦积分，余弦si(x)=/半dt函数ci(x)=i等dt要求：(1)编写函数，对任意给定的x,可求出值。(2)使用尽可能少的正余弦函数的调用，计算更精确的值。(用多种方法，创新方法)二、问题分析1、数值积分基本原理：用数值分析求解积分的数值方法有很多，如简单的梯形法、矩形法、辛普森(Simpson)法、牛顿-科斯特(Newton-Cotes)法等都是常用的方法。它们的基本思想都是将整个积分区间a,b分成n个子区间xi,xi+i,i=1,2,n,其中xi=a,xn+i=b。这样求定积分问题就分解为求和问题。2、本题要求用数值积分法计算正弦积分函数和余弦函

3、数积分，那么应该从编写函数的入手，建立function的m文件，通过对函数的调用，对任意跟定的x值，求出积分函数的值。三、数值积分法求解问题1、梯形公式、矩形公式首先，积分中值定理告诉我们，在积分区间a,b内存在一点%成立ff(x)dx=(b-a)f(，)，就是说，底为b-a而高为f(，)的矩a形面积恰等于所求区边梯形的面积。如果我们用两端点“高度”f(a)与f(b)的算术平均值作为平均高度f(D的近似值，这样导出的求积公式jf(x)dxb2af(a)+f(b)便是我们熟悉的梯形公式。将积分区间a,bn等分，每个小区间宽度均为h=*吃，h称为积布n步长。记a=X0<Xi<,<

4、Xk,<Xo=b,在小区间上用小矩形面积近似小曲边梯形的面积，若分别取左端点和右端点的函数值为小矩形的高，则分别得到两个曲边梯形面积的近似计算公式。具体程序如下：clearx=linspace(0,pi);dx=x(2);y=sin(x);s1=sum(y)*dxs2=trapz(y)*dxsc1=cumsum(y)*dx;sc2=cumtrapz(y)*dx;plot(x,-cos(x)+1,x,sc1,'.',x,sc2,'o')holdon由图可知这种方法精度太低，应选择其他方法。2、quad函数、quanl函数正弦：functiony=si(t)a

5、=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:100y2(i)=si(x(i);endplot(x,y2,'r')title('辛普森')1.91.81.71.61.51.41.31.21.110.901020

6、30405060708090100辛普森x=1:100;fori=1:100y2(i)=ci(x(i);endplot(x,y2,'b')title('辛普森')辛普森1400120010008006004002000-200-4000102030405060708090100给定任意X值，均可计算出对应的正弦、余弦函数积分。但从结果可以看出精度不是很高。3、复合求积公式由于牛顿-科特斯公式在nA8时不具有稳定性，故不可能通过提高阶的方法来提高求积精度。为了提高精度通常可把积分区间分成若干子区间（通常是等分），再在每个子区间上用低级求积公式。这种方法为复合求积法

7、。3.3.1复合梯形公式，n,在每个子将区间a,b划分为n等分，分点Xk=a+kh,h=b-a,k=n区间Xk，Xk¥tk=0,1;n.1,)上采用梯形公式，则得bnJx卜nI=f(x)dx=£1+f(x)dx=c£If(xk)+f(xkQ】+R(f)a*k0k2kz0n1n1Tn=-Hf(xk)+f(xk+)=【f(a)+2£f(xk)+”),2kz0称为复合梯形公式。复合梯形公式的余项nfFk)J1h3Rfkz0_12由于f(x)WC2a,b,且n1"I一"minfkfkmax0-：k-：n1-10:kn-1nk=0所以非w(a,

8、b便1nf'1f''knk=0于是复合梯形公式的余项为事实上只要设fxca,bi,则可得收敛性，只要把Tn改写成为程序如下:1ba,.b_a/viTn二2一二k/xk；kJxk正弦：functionT_n=fhtxs(a,b,n)h=(b-a)/n;fork=0:nx(k+1)=a+k*h;ifx(k+1)=0x(k+1)=10A(-10);endendT_1=h/2*(SS(x(1)+SS(x(n+1);fori=2:nF(i)=h*SS(x(i);endT_2=sum(F);T_n=T_1+T_2;余弦：functionT_n=fhtxc(a,b,n)h=(b-a)

9、/n;fork=0:nx(k+1)=a+k*h;ifx(k+1)=0x(k+1)=10A(-10);endendT_1=h/2*(CC(x(1)+CC(x(n+1);fori=2:nF(i)=h*CC(x(i);endT_2=sum(F);T_n=T_1+T_2;图像：余弦正弦复化梯形2i1-1.81.6八.-!X/V1.41.21，0.8''''010203040503.3.2复合新普斯求积公式将区间a,b划分为n等分，在每个子区间kk,xk41上采用辛普森公式，若记1 .Xk+2=Xk+h,则信2bn-1I=f(x)dx=1f(x)dxak=0n.1f(xk

10、)4f(4.12)f(xk)R(f).6称为复合辛普森求积公式。程序如下：正弦functionS_n=fhxpss(a,b,n)h=(b-a)/n;fork=0:nx(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);endendS_1=h/6*(SS(x(1)+SS(x(n+1);fori=2:nF_1(i)=h/3*SS(x(i);endforj=1:nF_2(j)=2*h/3*SS(x_k(j);endS_2=sum(F_1)+sum(F_2);S_n=S

11、_1+S_2;余弦：functionS_n=fhxpsc(a,b,n)h=(b-a)/n;fork=0:nx(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);endendS_1=h/6*(CC(x(1)+CC(x(n+1);fori=2:nF_1(i)=h/3*CC(x(i);endforj=1:nF_2(j)=2*h/3*CC(x_k(j);endS_2=sum(F_1)+sum(F_2);S_n=S_1+S_2;图像与复合梯形所得图像基本相同，深入分析两只

12、复合函数的优劣，对于积分函数si(x)=广；tdt假设x=1,则将区间0,1划分为8等份,应用复合梯形求得T8=0.9456909而如果将0,1分为4等份，应用复合辛普森有S4=0.9460832通过参考数值分析(李庆阳)的结论，发现无论是复合梯形公式还是复合辛普森公式，最终结果都会随着h值的减小而更加精确。对复合梯形公式和复合辛普森公式计算出的结果进行比较，发现复合梯形法的结果丁8只有两位有效数字，而复合辛普森的结果却有六位有效数字，所以复合辛普森公式计算出的结果更加的精确。4、插值型的求积公式clc,clearx0=0:0.5:5;y0=Inf1.75520.54030.0472-0.20

13、81-0.3205-0.3300-0.2676-0.1634-0.04680.0567;%所求积分函数的数值pp=csape(x0,y0);%默认的边界条件，Lagrange边界条件formatlonggchazhi=pp.coefs%显示每个区间上三次多项式的系数s=quadl(t)ppval(pp,t),0,5)%t积分format%恢复短小数的显示格式x=0:0.1:5;y=cos(x)/x;y1=spline(x0,y0,x);z=0*x;holdonplot(x,z,x,y,'k-',x,y1,'r')plot(x0,y0,'*')ho

14、ldoffclear0.0472-0.2081-0.3205-0.3300-0.2676%所求积分函数的数值%默认的边界条件，Lagrange边界条件%显示每个区间上三次多项式的系数%求积分%恢复短小数的显示格式x0=0:0.5:5;y0=Inf1.75520.5403-0.1634-0.04680.0567;pp=csape(x0,y0);formatlonggchazhi=pp.coefss=quadl(t)ppval(pp,t),0,5)formatx=0:0.1:5;y=cos(x)/x;y1=spline(x0,y0,x);z=0*x;holdonplot(x,z,x,y,'

15、k-',x,y1,'r')plot(x0,y0,'*')holdoff如图所示:5、高斯求积公式functionql,Ak,xk=gsqj(fun,a,b,n,tol)ifnargin=1a=-1;b=1;n=7;tol=1e-8;elseifnargin=3n=7;tol=1e-8;elseifnargin=4tol=1e-8;elseifnargin=2|nargin>5error('TheNumberofInputArgumentsIsWrong!');end%计算求积节点symsxp=sym2poly(diff(xA2-1)A(n+1),n+1)/(2An*factorial(n);tk=roots(p);%求积节点%计算求积系数Ak=zeros(n+1,1);fori=1:n+1xkt=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,1xk=(b-a)/2*tk+(b+a)/2;%检验积分函数fun有效性fun=fcnchk(fun,'vectorize');%计算变量代换之后