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

1、 課 程 設(shè) 計(jì) 任 務(wù) 書(shū) 課程設(shè)計(jì)課題： 用數(shù)值積分法計(jì)算正弦積分函數(shù)和余弦積分函數(shù)一、課程設(shè)計(jì)工作日自 2016 年 7 月 4 日至 2016 年 7 月 5日二、同組學(xué)生： 無(wú) 三、課程設(shè)計(jì)任務(wù)要求（包括課題來(lái)源、類型、目的和意義、基本要求、完成時(shí)間、主要參考資料等）：課題來(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)題敘述用

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

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

4、函數(shù)值為小矩形的高,則分別得到兩個(gè)曲邊梯形面積的近似計(jì)算公式。具體程序如下：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) hold on 由圖可知這種方法精度太低，應(yīng)選擇其他方法。2、quad函數(shù)、quan1函數(shù)正弦：function y=si(t) a=1e-8; %函數(shù)在0點(diǎn)無(wú)界，去掉0點(diǎn) y=quad(sin(x)./x,a,t) y=quadl(si

5、n(x)./x,a,t)余弦：function y=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;for i=1:100 y2(i)=si(x(i);endplot(x,y2,r)title(辛普森)x=1:100;for i=1:100 y2(i)=ci(x(i);endplot(x,y2,b)title(辛普森) 給定任意x值，均可計(jì)算出對(duì)應(yīng)的正弦、余弦函數(shù)積分。但從結(jié)果可以看出精度不是很高。3、 復(fù)合求積公式由于牛頓-科特斯公式在n8時(shí)不具有穩(wěn)定性，故不可能通過(guò)提高階的

6、方法來(lái)提高求積精度。為了提高精度通常可把積分區(qū)間分成若干子區(qū)間（通常是等分），再在每個(gè)子區(qū)間上用低級(jí)求積公式。這種方法為復(fù)合求積法。3.3.1 復(fù)合梯形公式將區(qū)間劃分為n等分，分點(diǎn)在每個(gè)子區(qū)間上采用梯形公式，則得記， 稱為復(fù)合梯形公式。復(fù)合梯形公式的余項(xiàng)由于且所以使 于是復(fù)合梯形公式的余項(xiàng)為 事實(shí)上只要設(shè)，則可得收斂性，只要把改寫成為程序如下：正弦：function T_n=fhtxs(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)=10(-10); endendT_1=h/2*(SS(x(1)+SS(x(n+1);for

7、i=2:n F(i)=h*SS(x(i);endT_2=sum(F);T_n=T_1+T_2;余弦：function 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)=10(-10); endendT_1=h/2*(CC(x(1)+CC(x(n+1);for i=2:n F(i)=h*CC(x(i);endT_2=sum(F);T_n=T_1+T_2;圖像： 正弦 余弦3.3.2 復(fù)合新普斯求積公式將區(qū)間劃分為n等分，在每個(gè)子區(qū)間上采用辛普森公式，若記則得l 稱為復(fù)合辛普森求積公式。程序如下：正弦funct

8、ion 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)=10(-10); x_k(k+1)=10(-10); endendS_1=h/6*(SS(x(1)+SS(x(n+1);for i=2:n F_1(i)=h/3*SS(x(i);endfor j=1:n F_2(j)=2*h/3*SS(x_k(j);endS_2=sum(F_1)+sum(F_2);S_n=S_1+S_2;余弦：function S_n=fhxpsc(a,b,

9、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)=10(-10); x_k(k+1)=10(-10); endendS_1=h/6*(CC(x(1)+CC(x(n+1);for i=2:n F_1(i)=h/3*CC(x(i);endfor j=1:n F_2(j)=2*h/3*CC(x_k(j);endS_2=sum(F_1)+sum(F_2);S_n=S_1+S_2;圖像與復(fù)合梯形所得圖像基本相同，深入分析兩只復(fù)合函數(shù)的優(yōu)劣，對(duì)于積分函數(shù) 假設(shè)x=1，則將區(qū)

10、間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é)果T8只有兩位有效數(shù)字，而復(fù)合辛普森的結(jié)果卻有六位有效數(shù)字，所以復(fù)合辛普森公式計(jì)算出的結(jié)果更加的精確。4、 插值型的求積公式clc, clearx0=0:0.5:5; y0= Inf 1.7552 0.5403 0.0472 -0.2081 -0.3205 -0.3300 -0.2676

11、 -0.1634 -0.0468 0.0567;%所求積分函數(shù)的數(shù)值pp=csape(x0,y0) ; %默認(rèn)的邊界條件，Lagrange邊界條件format long gchazhi=pp.coefs %顯示每個(gè)區(qū)間上三次多項(xiàng)式的系數(shù)s=quadl(t)ppval(pp,t),0,5) %求積分format %恢復(fù)短小數(shù)的顯示格式x=0:0.1:5;y=cos(x)/x;y1=spline(x0,y0,x);z=0*x;hold onplot(x,z,x,y,k-,x,y1,r)plot(x0,y0,*)hold offclearx0=0:0.5:5; y0= Inf 1.7552 0.54

12、03 0.0472 -0.2081 -0.3205 -0.3300 -0.2676 -0.1634 -0.0468 0.0567; %所求積分函數(shù)的數(shù)值pp=csape(x0,y0) ; %默認(rèn)的邊界條件，Lagrange邊界條件format long gchazhi=pp.coefs %顯示每個(gè)區(qū)間上三次多項(xiàng)式的系數(shù)s=quadl(t)ppval(pp,t),0,5) %求積分format %恢復(fù)短小數(shù)的顯示格式x=0:0.1:5;y=cos(x)/x;y1=spline(x0,y0,x);z=0*x;hold onplot(x,z,x,y,k-,x,y1,r)plot(x0,y0,*)ho

13、ld off如圖所示：5、 高斯求積公式function ql,Ak,xk=gsqj(fun,a,b,n,tol)if nargin=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|nargin5 error(The Number of Input Arguments Is Wrong!);end% 計(jì)算求積節(jié)點(diǎn)syms xp=sym2poly(diff(x2-1)(n+1),n+1)/(2n*factorial(n);tk=roots(p); % 求積節(jié)點(diǎn)% 計(jì)算求積系數(shù)Ak=zeros(n+1,1);for i=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); % 求積系數(shù)end% 積分變量代換，將a,b變換到-1,1xk=(b-a)/2*tk+(b+a)/2;% 檢驗(yàn)積分函數(shù)fun有效性fun=fcnchk(fun,vectorize);% 計(jì)算變量代換之后積分函數(shù)的值SS=fun(xk)*(b-a)/2;% 計(jì)算積分值ql=su