數(shù)值分析課程設(shè)計積分方法的實際應(yīng)用

1、數(shù)值分析課程設(shè)計報告求積公式的實際應(yīng)用學(xué) 院數(shù)學(xué)與統(tǒng)計學(xué)院專 業(yè)信息與計算科學(xué)學(xué) 號姓 名指導(dǎo)教師成 績教師評語：指導(dǎo)教師簽字： 2018年 1 月 8 日1 緒論數(shù)值分析是計算數(shù)學(xué)的一個主要部分，計算數(shù)學(xué)是數(shù)學(xué)科學(xué)的一個分支，它研究用計算機求解各種數(shù)學(xué)問題的數(shù)值檢索方其理論與軟件的實現(xiàn)。隨著計算機和計算方法的飛速發(fā)展，幾乎所有學(xué)科都走向定量化和精確化，從而產(chǎn)生了一系列計算性的學(xué)科分支，如計算物理、計算化學(xué)、計算生物學(xué)、計算地質(zhì)學(xué)、計算氣象學(xué)和計算材料學(xué)等，計算數(shù)學(xué)中的數(shù)值計算方法則是解決“計算”問題的橋梁和工具。我們知道，計算能力是計算工具和計算方法的效率的乘積，提高計算方法的效率與提高計算

2、機硬件的效率同樣重要?？茖W(xué)計算已用到科學(xué)技術(shù)和社會生活的各個領(lǐng)域中。數(shù)值計算方法，是一種研究并解決數(shù)學(xué)問題的數(shù)值近似解方法， 是在計算機上使用的解數(shù)學(xué)問題的方法，簡稱計算方法。在科學(xué)研究和工程技術(shù)中都要用到各種計算方法。 例如，在航天航空、地質(zhì)勘探、汽車制造、橋梁設(shè)計、天氣預(yù)報和漢字字樣設(shè)計中都有計算方法的蹤影。計算方法既有數(shù)學(xué)類課程中理論上的抽象性和嚴(yán)謹(jǐn)性，又有實用性和亞丁實驗性的技術(shù)特征，計算方法是一門理論性和實踐性都很強的學(xué)科。在70年代，大多數(shù)學(xué)校僅在數(shù)學(xué)系的計算數(shù)學(xué)專業(yè)和計算機系開設(shè)計算方法這門課程。隨著計算機技術(shù)的迅速發(fā)展和普及，現(xiàn)在計算方法課程幾乎已成為所有理工科學(xué)生的必修課程。

3、計算方法的計算對象是微積分，線性代數(shù)，常微分方程中的數(shù)學(xué)問題。內(nèi)容包括：插值和擬合、數(shù)值微分和數(shù)值積分、求解線性方程組的直接法和迭代法、計算矩陣特征值和特征向量和常微分方程數(shù)值解等問題。2 Gauss求積公式2.1 基本原理求積公式(2.1)含有2n+2個待定參數(shù)， 當(dāng)為等距節(jié)點時得到的插值求積公其代數(shù)精度至少為n次，如果適當(dāng)選取，有可能使求積公式(2.1)具有2n+1次代數(shù)精度，這類求積公式稱為高斯求積公式。為具有一般性，研究帶權(quán)積分，這里為權(quán)函數(shù)，求積公式為(2.2)為不依賴于的求積系數(shù)，為求積節(jié)點，可適當(dāng)選取及，使(2.2)具有2n+1次代數(shù)精度。如果求積公式(2.2)具有2n+1次代數(shù)

4、精度，則稱其節(jié)點為高斯點，相應(yīng)求積公式(2.2)稱為高斯求積公式。根據(jù)定義要使(2.2)式具有2n+1次代數(shù)精度，只要對，()，令(2.2)式精確成立，即 .(2.3)當(dāng)給定權(quán)函數(shù)，求出右端積分，則可由(2.3)式解得及2.2 程序?qū)崿F(xiàn)建立gaussl.m文件，寫入如下內(nèi)容：function s=gaussl(a,b,n)h=(b-a)/n; s=0.0;for m=0:(1*n/2-1) s=s+h*(gaussf(a+h*(1-1/sqrt(3)+2*m)+gaussf(a+h*(1+1/sqrt(3)+2*m);end 2.3 實例分析例計算積分 解建立gaussf.m文件以調(diào)用gaus

5、sl.m文件中的函數(shù)，再寫入如下內(nèi)容：function y=gaussf(x)y=sqrt(x)*log(x);再在命令行中輸入：>> s=gaussl(0,1,20)得出如下結(jié)果：s = -0.44563 高斯-勒讓德求積公式3.1 基本原理在高斯求積公式(2.1)中，若取權(quán)函數(shù)，區(qū)間為，則得公式(3.1)由于勒讓德多項式是區(qū)間 上的正交多項式，因此，勒讓德多項式的零點就是求積公式(3.1)的高斯點。形如(3.1)式的高斯公式稱為高斯-勒讓德求積公式。3.2 程序?qū)崿F(xiàn)建立guasslegendre.m文件，寫入如下內(nèi)容：function ql,Ak,xk=guasslegendr

6、e(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|nargin>5 error('The Number of Input Arguments Is Wrong!');end% 計算求積節(jié)點syms xp=sym2poly(diff(x2-1)(n+1),n+1)/(2n*factorial(n);tk=roots(p); % 求積節(jié)點% 計算求積系數(shù)Ak=zeros(n+1,1);fo

7、r 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;% 檢驗積分函數(shù)fun的有效性fun=fcnchk(fun,'vectorize');% 計算變量代換之后積分函數(shù)的值fx=fun(xk)*(b-a)/2;% 計算積分值ql=sum(Ak.*fx);參數(shù)說明：fun：積分表達(dá)式，可以是函數(shù)句柄a,b：積分上下限n：積分階

8、數(shù)tol：積分精度，默認(rèn)1e-6ql：積分結(jié)果Ak:積分系數(shù)xk：求積節(jié)點，滿足ql=sum(Ak.*fun(xk)3.3 實例分析例用4點的高斯-勒讓德公式求解定積分的近似值。解：打開guasslegendre.m文件，并在命令行中輸入如下內(nèi)容>> syms x;>> fun=inline(cos(x)*x2);>> ql,Ak,xk=guasslegendre(fun,0,pi/2,4)得出結(jié)果：ql = 0.4674Ak = 0.5689 0.2369 0.4786 0.2369 0.4786xk = 0.7854 0.0737 0.3625 1.49

9、71 1.2083即的4點的高斯-勒讓德積分結(jié)果為ql=0.4674。4 復(fù)化Simpson求積公式4.1 基本原理 復(fù)化Simpson公式是一種比較實用的積分方法，可以給出誤差估計。首先將區(qū)間a,b N等分，子區(qū)間的長度為 (4.1)在每個子區(qū)間上采用Simpson公式，在用Simpson公式時，還要將子區(qū)間再二等分，因此有2N+1個分點。即 (4.2)經(jīng)推導(dǎo)得到， (4.3)稱為為復(fù)化Simpson值，稱（4.3）式為復(fù)化Simpson公式。4.2 程序?qū)崿F(xiàn) 編寫復(fù)化Simpson求積函數(shù)（函數(shù)名：s_quad.m）Function I=S_quad(x,y)% 復(fù)化求積公式% x為被積函

10、數(shù)自變量的等距節(jié)點；y為被積函數(shù)在節(jié)點處的函數(shù)值。n=length(x);m=length(y); % 積分自變量的節(jié)點數(shù)應(yīng)與它的函數(shù)值個數(shù)相同；if n=m error ('The length of X and Y must be equal'); return;endif rem(n-1,2)=0 % 如果n-1不能被2整除，則調(diào)用復(fù)化公式 error ('節(jié)點數(shù)不滿足要求'); return;endN=(n-1)/2;h=(x(n)-x(1)/N;a=zeros (1,n);for k=1:Na(2*k-1)=a(2*k-1)+1;a(2*k)=a(2*

11、k)+4; a(2*k+1)=a(2*k+1)+1;endI=h/6*sum(a.*y);然后調(diào)用s_quad函數(shù)，來實現(xiàn)復(fù)化Simpson公式法。建立一個文件SPS，內(nèi)容如下：clearx=input('請輸入積分上下限及點間的間隔（例如-1:0.1:1）：');y=input('請輸入被積公式：y=');I=S_quad(x,y);disp('得出積分值I=')disp(I);4.3 實例分析例1 用復(fù)化Simpson公式求積分，在積分區(qū)間中點與點之間的間隔取為 0.1。解：運行程序，按照提示輸入積分上下限、點間的間隔及被積公式，如下所示：請

12、輸入積分上下限及點間的間隔（例如-1:0.1:1）：-1:0.1:1請輸入被積公式：y=exp(-x.2)得出積分值I= 1.4936 真值為：1.4937例2 計算積分，將區(qū)間8等分。解：運行程序，按照提示輸入積分上下限、等分后的區(qū)間長度及被積公式，如下 所示：請輸入積分上下限及點間的間隔（例如-1:0.1:1）：0:0.125:1請輸入被積公式：y=x./(4+x.2)得出積分值I= 0.1116真值為：0.1115724.4 結(jié)果分析復(fù)化Simpson計算所得的結(jié)果誤差較小，精度較高，更適合科學(xué)計算與應(yīng)用，且公式具有收斂性，穩(wěn)定性良好。5 數(shù)值方法的實際應(yīng)用 在實際問題中，往往會遇到一些

13、困難。有些函數(shù)找不到用初等函數(shù)表示的原函數(shù)，例如，對于積分 (5.1)而言，不存在用初等函數(shù)表示的原函數(shù)。而有些函數(shù)雖然能找到原函數(shù)，但計算過于復(fù)雜，例如，橢圓型積分 (5.2) 而有些情況下，只能知道某些點處的函數(shù)值，并沒有函數(shù)的具體表達(dá)式。這些情況，使我們有必要研究積分的數(shù)值計算問題。下面我們就以梯形公式為例做以說明。 所謂梯形求積公式就是用梯形面積來近似曲邊梯形面積，利用梯形公式和連續(xù)增加a,b的區(qū)間數(shù)來逼近： (5.3)第j次循環(huán)在個等距節(jié)點處對采樣。5.1 實例分析 衛(wèi)星軌道是一個橢圓，橢圓周長計算公式是，這里a是橢圓半長軸，c是地球中心與軌道中心（橢圓中心）的距離，記h為近地點距離

14、，H為遠(yuǎn)地點距離，R=637km為地球半徑，則我國第一顆人造衛(wèi)星近地點距離h=439km，遠(yuǎn)地點距離H=2384，試求衛(wèi)星軌道的周長。解: 第一步：先利用Matlab畫出被積函數(shù)的圖形。 輸入程序如下：clearH=2384;h=439;R=6371;a=(2*R+H+h)/2c=(H-h)/2x=0:0.1:pi/2;y=sqrt(1-(c/a)2*(sin(x).2);plot(x,y,'-')title('梯形法則');xlabel('x');ylabel('y');輸出結(jié)果：a = 7.782500000000000e+0

15、03c = 9.725000000000000e+002輸出圖形：圖5.1 被積函數(shù)的圖形 第二步：應(yīng)用數(shù)值積分梯形公式。 首先建立一個名為trapezg.m的M文件，程序如下：function I=trapezg(f_name3,a,b,n)format long%輸出用15位數(shù)字表示n=n;h=(b-a)/n; x=a+(0:n)*h;f=feval(f_name3,x);I=h/2*(f(1)+f(n+1);if n>1 I=I+h*sum(f(2:n);endh1=(b-a)/100;xc=a+(0:100)*h1; fc=feval(f_name3,xc);plot(xc,fc

16、,'r');hold onxlabel('x');ylabel('y');plot(x,f)title('數(shù)值積分梯形效果圖');plot(x,zeros(size(x),'.')for i=1:n;plot(x(i),x(i),0,f(i),end 然后建立一個名為f_name3.m的M文件定義函數(shù)，Matlab命令如下：function y=f_name3(x)y=sqrt(1-(9.725000e+002/7.782500e+003)2*(sin(x).2)-0.99; 輸入命令程序： >> tr

17、apezg('f_name3',0,pi/2,30) 輸出結(jié)果： ans = 0.00955791054630輸出圖形：圖5.2 數(shù)值積分效果圖積分結(jié)果為：0.00955791054630+0.99=0,99955791054630第三步：計算最后結(jié)果： 第四步：考慮誤差。clearn=1;format longfprintf('n Extended Trapezoidal Rulen');fprintf('n n I Errorn');I2=0.00955791054630;for k=1:8n=n*2;I1=trapezg('f_na

18、me3',0,pi/2,n);format longif k=1;fprintf('%3.0f %10.8f %10.8fn', n, I1, I1-I2);endpauseend計算7步輸出結(jié)果：Extended Trapezoidal Rulen I Error4 0.00956 0.00000 8 0.00956 0.0000016 0.00956 0.0000032 0.00956 -0.0000064 0.00956 -0.00000128 0.00956 0.00000256 0.00956 -0.00000初始狀態(tài)圖：圖5.3 初始狀態(tài)圖計算一步結(jié)果圖：圖5.4 計算一步結(jié)果圖計算四步結(jié)果圖：圖5.