數(shù)值分析法求正弦余弦積分函數(shù)

1、天津職業(yè)技術(shù)師范大學(xué)課程設(shè)計(jì)任務(wù)書理學(xué)院數(shù)學(xué)1403班學(xué)生.張群課程設(shè)計(jì)課題：用數(shù)值積分法計(jì)算正弦積分函數(shù)和余弦積分函數(shù)一、課程設(shè)計(jì)工作日自2016年7月4日至2016年7月5日二、同組學(xué)生：無(wú)三、課程設(shè)計(jì)任務(wù)要求（包括課題來(lái)源、類型、目的和意義、基本要求、完成時(shí)問(wèn)、主要參考資料等）：課題來(lái)源：教師自擬類型：理論研究目的和意義：培養(yǎng)學(xué)生對(duì)數(shù)值分析中主要算法的應(yīng)用能力，探索相關(guān)算法之間的內(nèi)在聯(lián)系。基本要求：根據(jù)數(shù)值分析課程所學(xué)的知識(shí)，應(yīng)用Matlab軟件編寫程序，完成對(duì)算法及其內(nèi)在原理的實(shí)驗(yàn)研究。完成時(shí)間：參考資料：數(shù)值分析李慶揚(yáng)等清華大學(xué)出版社指導(dǎo)教師簽字:教研室主任簽字:一、問(wèn)題敘述用數(shù)值積

2、分法計(jì)算正弦積分函數(shù)和余弦積分函數(shù)提示：正弦積分，余弦si(x)=/半dt函數(shù)ci(x)=i等dt要求：(1)編寫函數(shù)，對(duì)任意給定的x,可求出值。(2)使用盡可能少的正余弦函數(shù)的調(diào)用，計(jì)算更精確的值。(用多種方法，創(chuàng)新方法)二、問(wèn)題分析1、數(shù)值積分基本原理：用數(shù)值分析求解積分的數(shù)值方法有很多，如簡(jiǎn)單的梯形法、矩形法、辛普森(Simpson)法、牛頓-科斯特(Newton-Cotes)法等都是常用的方法。它們的基本思想都是將整個(gè)積分區(qū)間a,b分成n個(gè)子區(qū)間xi,xi+i,i=1,2,n,其中xi=a,xn+i=b。這樣求定積分問(wèn)題就分解為求和問(wèn)題。2、本題要求用數(shù)值積分法計(jì)算正弦積分函數(shù)和余弦函

3、數(shù)積分，那么應(yīng)該從編寫函數(shù)的入手，建立function的m文件，通過(guò)對(duì)函數(shù)的調(diào)用，對(duì)任意跟定的x值，求出積分函數(shù)的值。三、數(shù)值積分法求解問(wèn)題1、梯形公式、矩形公式首先，積分中值定理告訴我們，在積分區(qū)間a,b內(nèi)存在一點(diǎn)%成立ff(x)dx=(b-a)f(，)，就是說(shuō)，底為b-a而高為f(，)的矩a形面積恰等于所求區(qū)邊梯形的面積。如果我們用兩端點(diǎn)“高度”f(a)與f(b)的算術(shù)平均值作為平均高度f(wàn)(D的近似值，這樣導(dǎo)出的求積公式j(luò)f(x)dxb2af(a)+f(b)便是我們熟悉的梯形公式。將積分區(qū)間a,bn等分，每個(gè)小區(qū)間寬度均為h=*吃，h稱為積布n步長(zhǎng)。記a=X0<Xi<,<

4、Xk,<Xo=b,在小區(qū)間上用小矩形面積近似小曲邊梯形的面積，若分別取左端點(diǎn)和右端點(diǎn)的函數(shù)值為小矩形的高，則分別得到兩個(gè)曲邊梯形面積的近似計(jì)算公式。具體程序如下：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由圖可知這種方法精度太低，應(yīng)選擇其他方法。2、quad函數(shù)、quanl函數(shù)正弦：functiony=si(t)a

5、=1e-8;%函數(shù)在0點(diǎn)無(wú)界，去掉0點(diǎn)y=quad('sin(x)./x',a,t)y=quadl('sin(x)./x',a,t)余弓玄:functiony=ci(t)a=-1e1;%函數(shù)在0點(diǎn)無(wú)界，去掉0點(diǎn)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值，均可計(jì)算出對(duì)應(yīng)的正弦、余弦函數(shù)積分。但從結(jié)果可以看出精度不是很高。3、復(fù)合求積公式由于牛頓-科特斯公式在nA8時(shí)不具有穩(wěn)定性，故不可能通過(guò)提高階的方法來(lái)提高求積精度。為了提高精度通?？砂逊e分區(qū)間分成若干子區(qū)間（通常是等分），再在每個(gè)子區(qū)間上用低級(jí)求積公式。這種方法為復(fù)合求積法

7、。3.3.1復(fù)合梯形公式，n,在每個(gè)子將區(qū)間a,b劃分為n等分，分點(diǎn)Xk=a+kh,h=b-a,k=n區(qū)間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稱為復(fù)合梯形公式。復(fù)合梯形公式的余項(xiàng)nfFk)J1h3Rfkz0_12由于f(x)WC2a,b,且n1"I一"minfkfkmax0-：k-：n1-10:kn-1nk=0所以非w(a,

8、b便1nf'1f''knk=0于是復(fù)合梯形公式的余項(xiàng)為事實(shí)上只要設(shè)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;圖像：余弦正弦復(fù)化梯形2i1-1.81.6八.-!X/V1.41.21，0.8''''010203040503.3.2復(fù)合新普斯求積公式將區(qū)間a,b劃分為n等分，在每個(gè)子區(qū)間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稱為復(fù)合辛普森求積公式。程序如下：正弦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;圖像與復(fù)合梯形所得圖像基本相同，深入分析兩只

12、復(fù)合函數(shù)的優(yōu)劣，對(duì)于積分函數(shù)si(x)=廣；tdt假設(shè)x=1,則將區(qū)間0,1劃分為8等份,應(yīng)用復(fù)合梯形求得T8=0.9456909而如果將0,1分為4等份，應(yīng)用復(fù)合辛普森有S4=0.9460832通過(guò)參考數(shù)值分析(李慶陽(yáng))的結(jié)論，發(fā)現(xiàn)無(wú)論是復(fù)合梯形公式還是復(fù)合辛普森公式，最終結(jié)果都會(huì)隨著h值的減小而更加精確。對(duì)復(fù)合梯形公式和復(fù)合辛普森公式計(jì)算出的結(jié)果進(jìn)行比較，發(fā)現(xiàn)復(fù)合梯形法的結(jié)果丁8只有兩位有效數(shù)字，而復(fù)合辛普森的結(jié)果卻有六位有效數(shù)字，所以復(fù)合辛普森公式計(jì)算出的結(jié)果更加的精確。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;%所求積分函數(shù)的數(shù)值pp=csape(x0,y0);%默認(rèn)的邊界條件，Lagrange邊界條件formatlonggchazhi=pp.coefs%顯示每個(gè)區(qū)間上三次多項(xiàng)式的系數(shù)s=quadl(t)ppval(pp,t),0,5)%t積分format%恢復(fù)短小數(shù)的顯示格式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%所求積分函數(shù)的數(shù)值%默認(rèn)的邊界條件，Lagrange邊界條件%顯示每個(gè)區(qū)間上三次多項(xiàng)式的系數(shù)%求積分%恢復(fù)短小數(shù)的顯示格式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%計(jì)算求積節(jié)點(diǎn)symsxp=sym2poly(diff(xA2-1)A(n+1),n+1)/(2An*factorial(n);tk=roots(p);%求積節(jié)點(diǎn)%計(jì)算求積系數(shù)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);%求積系數(shù)end%積分變量代換，將a,b變換到-1,1xk=(b-a)/2*tk+(b+a)/2;%檢驗(yàn)積分函數(shù)fun有效性fun=fcnchk(fun,'vectorize');%計(jì)算變量代換之后