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

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

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

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

4、精度，則稱其節(jié)點(diǎn)為高斯點(diǎn)，相應(yīng)求積公式(2.2)稱為高斯求積公式。根據(jù)定義要使(2.2)式具有2n+1次代數(shù)精度，只要對(duì)，()，令(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 實(shí)例分析例計(jì)算積分 解建立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)由于勒讓德多項(xiàng)式是區(qū)間 上的正交多項(xiàng)式，因此，勒讓德多項(xiàng)式的零點(diǎn)就是求積公式(3.1)的高斯點(diǎn)。形如(3.1)式的高斯公式稱為高斯-勒讓德求積公式。3.2 程序?qū)崿F(xiàn)建立guasslegendre.m文件，寫入如下內(nèi)容：function ql,Ak,xk=guasslegendre(fun,a,

6、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(

7、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ù)的值fx=fun(xk)*(b-a)/2;% 計(jì)算積分值ql=sum(Ak.*fx);參數(shù)說明：fun：積分表達(dá)式，可以是函數(shù)句柄a,b：積分上下限n：積分階數(shù)tol：積分精度，默認(rèn)1e-6ql：積分結(jié)果Ak:積分系數(shù)xk

8、：求積節(jié)點(diǎn)，滿足ql=sum(Ak.*fun(xk)3.3 實(shí)例分析例用4點(diǎn)的高斯-勒讓德公式求解定積分的近似值。解：打開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.4971 1.2083即的4點(diǎn)的高斯-勒讓德積分結(jié)果為ql=0.4674。4 復(fù)化Simpson求積公式4.1 基本

9、原理 復(fù)化Simpson公式是一種比較實(shí)用的積分方法，可以給出誤差估計(jì)。首先將區(qū)間a,b N等分，子區(qū)間的長(zhǎng)度為 (4.1)在每個(gè)子區(qū)間上采用Simpson公式，在用Simpson公式時(shí)，還要將子區(qū)間再二等分，因此有2N+1個(gè)分點(diǎn)。即 (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為被積函數(shù)自變量的等距節(jié)點(diǎn)；y為被積函數(shù)在節(jié)點(diǎn)處的函數(shù)值。n=length(x);m=length(y); % 積分自變

10、量的節(jié)點(diǎn)數(shù)應(yīng)與它的函數(shù)值個(gè)數(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é)點(diǎn)數(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*k)+4; a(2*k+1)=a(2*k+1)+1;endI=h/6*sum(a.*y);然后調(diào)用s_quad函數(shù)，來實(shí)現(xiàn)復(fù)化Simpson公式法。建立

11、一個(gè)文件SPS，內(nèi)容如下：clearx=input(請(qǐng)輸入積分上下限及點(diǎn)間的間隔（例如-1:0.1:1）：);y=input(請(qǐng)輸入被積公式：y=);I=S_quad(x,y);disp(得出積分值I=)disp(I);4.3 實(shí)例分析例1 用復(fù)化Simpson公式求積分，在積分區(qū)間中點(diǎn)與點(diǎn)之間的間隔取為 0.1。解：運(yùn)行程序，按照提示輸入積分上下限、點(diǎn)間的間隔及被積公式，如下所示：請(qǐng)輸入積分上下限及點(diǎn)間的間隔（例如-1:0.1:1）：-1:0.1:1請(qǐng)輸入被積公式：y=exp(-x.2)得出積分值I= 1.4936 真值為：1.4937例2 計(jì)算積分，將區(qū)間8等分。解：運(yùn)行程序，按照提示輸入

12、積分上下限、等分后的區(qū)間長(zhǎng)度及被積公式，如下 所示：請(qǐng)輸入積分上下限及點(diǎn)間的間隔（例如-1:0.1:1）：0:0.125:1請(qǐng)輸入被積公式：y=x./(4+x.2)得出積分值I= 0.1116真值為：0.1115724.4 結(jié)果分析復(fù)化Simpson計(jì)算所得的結(jié)果誤差較小，精度較高，更適合科學(xué)計(jì)算與應(yīng)用，且公式具有收斂性，穩(wěn)定性良好。5 數(shù)值方法的實(shí)際應(yīng)用 在實(shí)際問題中，往往會(huì)遇到一些困難。有些函數(shù)找不到用初等函數(shù)表示的原函數(shù)，例如，對(duì)于積分 (5.1)而言，不存在用初等函數(shù)表示的原函數(shù)。而有些函數(shù)雖然能找到原函數(shù)，但計(jì)算過于復(fù)雜，例如，橢圓型積分 (5.2) 而有些情況下，只能知道某些點(diǎn)處的

13、函數(shù)值，并沒有函數(shù)的具體表達(dá)式。這些情況，使我們有必要研究積分的數(shù)值計(jì)算問題。下面我們就以梯形公式為例做以說明。 所謂梯形求積公式就是用梯形面積來近似曲邊梯形面積，利用梯形公式和連續(xù)增加a,b的區(qū)間數(shù)來逼近： (5.3)第j次循環(huán)在個(gè)等距節(jié)點(diǎn)處對(duì)采樣。5.1 實(shí)例分析 衛(wèi)星軌道是一個(gè)橢圓，橢圓周長(zhǎng)計(jì)算公式是，這里a是橢圓半長(zhǎng)軸，c是地球中心與軌道中心（橢圓中心）的距離，記h為近地點(diǎn)距離，H為遠(yuǎn)地點(diǎn)距離，R=637km為地球半徑，則我國第一顆人造衛(wèi)星近地點(diǎn)距離h=439km，遠(yuǎn)地點(diǎn)距離H=2384，試求衛(wèi)星軌道的周長(zhǎng)。解: 第一步：先利用Matlab畫出被積函數(shù)的圖形。 輸入程序如下：clear

14、H=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+003c = 9.725000000000000e+002輸出圖形：圖5.1 被積函數(shù)的圖形 第二步：應(yīng)用數(shù)值積分梯形公式。 首先建立一個(gè)名為trapezg.m的M文件，程序如下：function I=trapezg(f_name3,a,b,n)format long%輸出用15位數(shù)字表示n=

15、n;h=(b-a)/n; x=a+(0:n)*h;f=feval(f_name3,x);I=h/2*(f(1)+f(n+1);if n1 I=I+h*sum(f(2:n);endh1=(b-a)/100;xc=a+(0:100)*h1; fc=feval(f_name3,xc);plot(xc,fc,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 然后建立一個(gè)名為f_name3.m的M文件定義函數(shù)，Matlab命令

16、如下：function y=f_name3(x)y=sqrt(1-(9.725000e+002/7.782500e+003)2*(sin(x).2)-0.99; 輸入命令程序： trapezg(f_name3,0,pi/2,30) 輸出結(jié)果： ans = 0.00955791054630輸出圖形：圖5.2 數(shù)值積分效果圖積分結(jié)果為：0.00955791054630+0.99=0,99955791054630第三步：計(jì)算最后結(jié)果： 第四步：考慮誤差。clearn=1;format longfprintf(n Extended Trapezoidal Rulen);fprintf(n n I Er

17、rorn);I2=0.00955791054630;for k=1:8n=n*2;I1=trapezg(f_name3,0,pi/2,n);format longif k=1;fprintf(%3.0f %10.8f %10.8fn, n, I1, I1-I2);endpauseend計(jì)算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

18、0.00956 -0.00000初始狀態(tài)圖：圖5.3 初始狀態(tài)圖計(jì)算一步結(jié)果圖：圖5.4 計(jì)算一步結(jié)果圖計(jì)算四步結(jié)果圖：圖5.5 計(jì)算四步結(jié)果圖最后計(jì)算八步結(jié)果圖：圖5.6 計(jì)算八步結(jié)果圖5.2 結(jié)果分析 數(shù)值微積分在科學(xué)計(jì)算上有很大應(yīng)用，復(fù)化梯形公式極大地簡(jiǎn)化了人們實(shí)際生活中的運(yùn)算復(fù)雜性。不僅算法簡(jiǎn)明，幾何意義明確，而且迭代結(jié)果準(zhǔn)確科學(xué)，具有更好的收斂性和廣泛的實(shí)用性，且精度高收斂速度快。本題就是利用復(fù)化積分原理，確定軌道，并加以計(jì)算，得到高精度的結(jié)果。結(jié)論通過本次課程設(shè)計(jì)，學(xué)習(xí)到了如何應(yīng)用數(shù)學(xué)軟件進(jìn)行大量的數(shù)值求解方法，對(duì)該課程的理解進(jìn)一步加強(qiáng)，對(duì)課堂知識(shí)的模糊點(diǎn)有了一個(gè)清晰的認(rèn)識(shí)。在本文敘述的幾種方法中，高斯-勒讓德求積公式可以用較少的節(jié)點(diǎn)數(shù)得到高精度的計(jì)算結(jié)果，是現(xiàn)實(shí)生活中經(jīng)常用到的數(shù)值積分方法，但是當(dāng)積分區(qū)間較大時(shí)，積分精度