數(shù)值積分用matlab實(shí)現(xiàn)

1、東北大學(xué)秦皇島分校數(shù)值計(jì)算課程設(shè)計(jì)報(bào)告數(shù)值積分及Matlab實(shí)現(xiàn)學(xué) 院數(shù)學(xué)與統(tǒng)計(jì)學(xué)院專 業(yè)信息與計(jì)算科學(xué)學(xué) 號(hào)5133117姓 名楚文玉指導(dǎo)教師張建波 姜玉山成 績(jī)教師評(píng)語(yǔ)：指導(dǎo)教師簽字： 2015年07月14日1 緒論在科研計(jì)算中，經(jīng)常會(huì)碰到一些很難用公式定理直接求出精確解的積分問(wèn)題，對(duì)于這類問(wèn)題，我們一般轉(zhuǎn)化為數(shù)值積分問(wèn)題，用計(jì)算機(jī)來(lái)實(shí)現(xiàn)求解問(wèn)題1.1 課題的背景對(duì)于定積分在求某函數(shù)的定積分時(shí)，在一定條件下，雖然有牛頓-萊布里茨公式可以計(jì)算定積分的值，但在很多情況下的原函數(shù)不易求出或非常復(fù)雜被積函數(shù)的原函數(shù)很難用初等函數(shù)表達(dá)出來(lái)，例如等；有的函數(shù)的原函數(shù)存在，但其表達(dá)式太復(fù)雜，計(jì)算量太大，

2、有的甚至無(wú)法有解析表達(dá)式因此能夠借助牛頓-萊布尼茲公式計(jì)算定積分的情形是不多的另外，許多實(shí)際問(wèn)題中的被積函數(shù)往往是列表函數(shù)或其他形式的非連續(xù)函數(shù)，對(duì)這類函數(shù)的定積分，也不能用不定積分方法求解，只能設(shè)法求其近似值因此，探討近似計(jì)算的數(shù)值積分方法是有明顯的實(shí)際意義的，即有必要研究定積分的數(shù)值計(jì)算方法，以解決定積分的近似計(jì)算而數(shù)值積分就是解決此類問(wèn)題的一種有效的方法，它的特點(diǎn)是利用被積函數(shù)在一些節(jié)點(diǎn)上的信息求出定積分的近似值微積分的發(fā)明是人類科學(xué)史上一項(xiàng)偉大的成就，在科學(xué)技術(shù)中，積分是經(jīng)常遇到的一個(gè)重要計(jì)算環(huán)節(jié)數(shù)值積分是數(shù)學(xué)上重要的課題之一，是數(shù)值分析中重要的內(nèi)容之一隨著計(jì)算機(jī)的出現(xiàn)，近幾十年來(lái)，對(duì)

3、于數(shù)值積分問(wèn)題的研究已經(jīng)成為一個(gè)很活躍的研究領(lǐng)域現(xiàn)在，數(shù)值積分在計(jì)算機(jī)圖形學(xué)，積分方程，工程計(jì)算，金融數(shù)學(xué)等應(yīng)用科學(xué)領(lǐng)域都有著相當(dāng)重要的應(yīng)用，所以研究數(shù)值積分問(wèn)題有著很重要的意義國(guó)內(nèi)外眾多學(xué)者在數(shù)值積分應(yīng)用領(lǐng)域也提出了許多新方法在很多實(shí)際應(yīng)用中，只能知道積分函數(shù)在某些特定點(diǎn)的取值，比如天氣測(cè)量中的氣溫、濕度、氣壓等，醫(yī)學(xué)測(cè)量中的血壓、濃度等等通過(guò)這個(gè)課題的研究，我們將會(huì)更好地掌握運(yùn)用數(shù)值積分算法求出特殊積分函數(shù)的定積分的一些基本方法、理論基礎(chǔ)；并且通過(guò)Matlab軟件編程的實(shí)現(xiàn)，應(yīng)用于實(shí)際生活中1.2 課題的主要內(nèi)容框架1.2.1 數(shù)值積分各求積公式簡(jiǎn)介簡(jiǎn)介牛頓-柯特斯求積公式及其辛普森求積公

4、式，龍貝格求積公式，高斯求積公式的基本理論基礎(chǔ)和方法1.2.2 求積公式的代碼實(shí)現(xiàn)通過(guò)理解各種數(shù)值積分求積公式的原理方法，通過(guò)Matlab軟件編程，實(shí)現(xiàn)以上求積公式1.2.3 應(yīng)用舉例通過(guò)簡(jiǎn)單舉例，自建一個(gè)相對(duì)簡(jiǎn)單和復(fù)雜的函數(shù)，用上面編寫的Matlab源程序來(lái)解決實(shí)際問(wèn)題，體會(huì)數(shù)值積分和Matlab的優(yōu)勢(shì)2 牛頓柯特斯公式及Matlab實(shí)現(xiàn)2.1 牛頓柯特斯公式的基本原理方法設(shè)將積分區(qū)間a, b劃分為n等分，步長(zhǎng)為，選取等距節(jié)點(diǎn)構(gòu)造出的差值型求積公式， （2.1）稱為牛頓-柯特斯公式，式中稱為柯特斯系數(shù)根據(jù)， （2.2）引進(jìn)變量代換，則有 （2.3）當(dāng)n = 2時(shí)，此時(shí)柯特斯系數(shù)為，相應(yīng)的求積

5、公式就是辛普森求積公式： （2.4）2.2 牛頓柯特斯公式的Matlab實(shí)現(xiàn)functionC, g = NCotes(a, b, n, m)% a，b分別為積分的上下限；% n是子區(qū)間的個(gè)數(shù)； % m是被調(diào)用第幾個(gè)被積函數(shù)； % 當(dāng)n=1時(shí)計(jì)算梯形公式；當(dāng)n=2時(shí)計(jì)算辛普森公式，以此類推；I = n; h = (b - a) / i; z = 0; for j = 0 : i x(j + 1) = a + j * h; s = 1; if j = 0 s = s; elsefor k=1 : j s =s * k; endendr = 1; if i - j = 0 r = r; else

6、for k = 1 : (I - j) r = r * k; end end if mod(I - j), 2) = 1 q = -(I * s * r); else q = i * s * r; endy = 1; for k = 0 : i if k = j y = y * (sym('t') - k); endend l = int(y, 0 , i); C(j + 1)= l / q; z = z + C(j + 1)*f1(m, x(j + 1); end g=(b - a)*z3 復(fù)合求積公式及Matlab實(shí)現(xiàn)3.1 復(fù)合梯形公式的基本原理將區(qū)間a, b劃分成n等分

7、，分點(diǎn)，,在每個(gè)子區(qū)間（）上采用梯形公式得： （3.1）記 （3.2）稱式(3.2)為復(fù)合梯形公式3.2 復(fù)合梯形公式的Matlab實(shí)現(xiàn)function s = trapr1(f, a, b, n)% f表示被積函數(shù)；% a，b表示積分上下限；% n是子區(qū)間的個(gè)數(shù)；h = (b - a) / n; s = 0; for k = 1 : (n - 1) x = a + h * k; s = s + feval('f', x); end format long s = h*(feval('f', a) + feval('f', b) / 2 + h

8、* s;3.3 復(fù)合辛普森求積公式的基本原理將區(qū)間a,b分等分，在每一個(gè)子區(qū)間上采用辛普森公式，若記，則得 （3.3）記（3.4）稱式(3.4)為復(fù)合辛普森求積公式3.4 復(fù)合辛普森求積公式的Matlab實(shí)現(xiàn)function s = simpr1(f, a, b, n) % f表示被積函數(shù)；% a，b表示積分上下限；% n是子區(qū)間的個(gè)數(shù)；h = (b - a) / (2 * n); s1 = 0; s2 = 0; for k = 1 : n x = a + h * (2*k - 1); s1 = s1 + feval('f', x); end for k = 1 : (n -

9、1) x = a + h * 2 * k; s2 = s2 + feval('f', x); end s = h*(feval('f', a) + feval('f', b) + 4 * s1 + 2 * s2) / 3;4 龍貝格求積公式及Matlab實(shí)現(xiàn)4.1 龍貝格算法的基本原理由梯形的遞推法可以看出，將積分區(qū)間等分時(shí)，用復(fù)化梯形公式計(jì)算的結(jié)果作為積分的近似值，其誤差近似值為 （4.1）可以設(shè)想，如果用這個(gè)誤差作為的一種補(bǔ)償，即將 （4.2）作為積分的近似值，可望提高其精確度直接根據(jù)復(fù)化求積公式，不難驗(yàn)證 （4.3）這說(shuō)明，將區(qū)間對(duì)分前后兩

10、次復(fù)化梯形公式的值，按式 （4.4）作線性組合恰好等于復(fù)合辛普森公式的值，它比更接近于近似值 同樣，根據(jù) （4.5）用于作線性組合會(huì)得到比更精確的值，且通過(guò)直接驗(yàn)證可得 （4.6）再由 （4.7）用與作線性組合，又可得到比更精確的值，通常記為，即 （4.8）此式（4.8）就稱為龍貝格求積公式 上述用若干積分近似值推算出更為精確的積分近似值得方法，稱為外推法我們將序列分別稱為梯形序列，辛普森序列，柯特斯序列和龍貝格序列由龍貝格序列求積的算法稱為龍貝格算法具體步驟為：第一步：算出和的值，根據(jù)公式 （4.9）求出;第二步：將區(qū)間a,b分半，算出的值，并根據(jù)(4.3)和(4.9)式計(jì)算和;第三步：再將

11、區(qū)間分半，算出及的值，并根據(jù)(4.3)和(4.9)式計(jì)算和，再有公式(4.4)求出；第四步：再將區(qū)間分半，計(jì)算，, ,并根據(jù)公式(4.5)計(jì)算第五步：再將區(qū)間分半，類似上述過(guò)程計(jì)算，, ，重復(fù)以上步驟即可得到， ,一直到龍貝格序列中前后兩項(xiàng)的絕對(duì)值差不超過(guò)給定的誤差險(xiǎn)為止4.2 龍貝格算法的Matlab實(shí)現(xiàn)function R, quad, err, h=romber(f, a, b, n, delta) % f表示被積函數(shù)；% a，b表示積分上下限；% n是子區(qū)間的個(gè)數(shù)% delta是誤差限M = 1; h = b - a; err = 1 J = 0; R = zeros(4, 4); R

12、(1, 1) = h * (feval('f', a) + feval('f', b) / 2 while (err > delta) & (J < n) ) | (J < 4) J = J + 1; h = h / 2; s = 0; for p = 1 : M x = a + h * (2*p - 1); s = s + feval('f', x); end R(J + 1, 1) = R(J, 1) / 2 + h*s; M = 2*M;for K = 1 : J R(J + 1,K + 1) = R(J + 1,

13、 K) + (R(J + 1,K) - R(J, K) / (4K - 1); end err = abs(R(J, J) - R(J + 1,K + 1); end quad = R(J + 1, J + 1)；5 高斯勒讓德求積公式及Matlab實(shí)現(xiàn)5.1 高斯勒讓德求積公式的基本原理在高斯求積公式 （5.1）中，若取權(quán)函數(shù)=1，區(qū)間-1,1，則得公式 （5.2）我們知道勒讓德多項(xiàng)式是區(qū)間-1,1上的正交多項(xiàng)式，因此，勒讓德多項(xiàng)式的零點(diǎn)就是求積公式(5.2)的高斯點(diǎn)，形如(5.1)式的高斯公式特別的稱為高斯-勒讓德求積公式如下表5.1所示為高斯-勒讓德求積公式的節(jié)點(diǎn)數(shù)和系數(shù)5.2 高斯勒讓

14、德求積公式的Matlab實(shí)現(xiàn)function quad = gauss8(f ,a, b, x, A)N = length(x);T = zeros(1, N);T = (a + b) / 2 + (b - a) / 2) * x;quad = (b - a) / 2) * sum(A. * feval('f', T); 表5.1 高斯-勒讓德求積公式的節(jié)點(diǎn)數(shù)和系數(shù)nn00.00000002.000000030.86113630.33998100.34785480.652145210.57735031.000000040.90619780.53846930.00000000.2

15、3692690.47862870.568888920.77459670.00000000.55555560.888888950.93246950.66120940.23861920.17132450.36076160.46791396 各個(gè)求積公式的應(yīng)用舉例與比較分析6.1 簡(jiǎn)單數(shù)值積分的解(精確值0.946083070367183)6.1.1 牛頓柯特斯當(dāng)n=1時(shí)的梯形算法和n=2時(shí)的辛普森算法的結(jié)果解：先用M文件定義一個(gè)f1.m的函數(shù)function f = f1(i, x) g(1) = sqrt(x);if x = 0 g(2) = 1; else g(2) = sin(x) / x;

16、 end f = g(i); 輸入>>Ncotes(0, 1, 1, 2) 回車得到；ans = 0.9270輸入>>Ncotes（0, 1, 2, 2） 回車得到：ans = 0.94616.1.2 復(fù)合梯形公式和復(fù)合辛普森求積公式的計(jì)算結(jié)果解：建立一個(gè)M文件定義一個(gè)f.m函數(shù)function y = f(x) if x = 0y = 1; else y = sin(x) / xend輸入>>trapr1(f, 0, 1, 10) 回車得到：ans = 0.9458輸入>>simpr1(f, 0, 1, 10)ans =0.94616.1.3

17、龍貝格求積公式的計(jì)算結(jié)果（取誤差不超過(guò)）解：建立一個(gè)M文件定義一個(gè)f.m函數(shù)function y = f(x) if x = 0y = 1; else y = sin(x) / xend輸入>>romber（f, 0, 1, 10，0.5*10(- 10)） 回車得到：quad = 0.9461ans = 0.9207 0 0 0 0 0.9398 0.9461 0 0 0 0.9445 0.9461 0.9461 0 0 0.9457 0.9461 0.9461 0.9461 00.9460 0.9461 0.9461 0.9461 0.94616.1.4 高斯勒讓德求積公式的計(jì)

18、算結(jié)果(給定節(jié)點(diǎn)3)解：建立一個(gè)M文件定義一個(gè)f.m函數(shù)function y = f(x) ；y = sin(x) / x;輸入>>gauss8（f,0,1,-0.7745966692,0.7745966692,0,0.5555555556,0.5555555556,0.8888888888）回車得到：ans = 0.89566.2 復(fù)雜數(shù)值積分的解(精確值-1.8785)6.2.1 牛頓柯特斯當(dāng)n = 1時(shí)的梯形算法和n = 2時(shí)的辛普森算法的結(jié)果解：先用M文件定義一個(gè)f1.m的函數(shù)function f = f1(i, x) g(1) = sqrt(x);if x = 0 g(2

19、) = 1; else g(2) = cos(x) 1 / (1 + x.2) 1 / (4.*sqrt(4 + x.2);endf = g(i); 輸入>>Ncotes(0, 2*pi, 1, 2)回車得到ans = -1.8692輸入>>Ncotes(0, 2*pi, 2, 2) 回車得到ans = -1.87046.2.2 復(fù)合梯形公式和復(fù)合辛普森求積公式的計(jì)算結(jié)果解：先用M文件定義一個(gè)f1.m的函數(shù)function y = f(x) if x = 0y = 1; else y = cos(x) 1 / (1 + x.2) 1 / (4.*sqrt(4 + x.2

20、);end輸入>>trapr1(f, 0, 2*pi, 10) 回車得到：ans = -1.8748輸入>>simpr1(f, 0, 2*pi, 10)ans = -1.86946.2.3 龍貝格求積公式的計(jì)算結(jié)果（取誤差不超過(guò)）解：先用M文件定義一個(gè)f1.m的函數(shù)function y = f(x) if x = 0y = 1; else y = cos(x) 1 / (1 + x.2) 1 / (4.*sqrt(4 + x.2);end輸入>>romber（f, 0, 2*pi ,10，0.5*10(-10)） 回車得到：quad = -1.87646.3

21、 各個(gè)求積公式的比較分析（以的各個(gè)積分結(jié)果為例） 表6.1各個(gè)數(shù)值積分的比較積分方法牛頓柯特斯（梯形）牛頓-柯特斯（辛普森）復(fù)合梯形復(fù)合辛普森龍貝格高斯-勒讓德精確值0.94610.94610.94610.94610.94610.9461實(shí)際值0.92700.94600.94580.94600.94600.8956誤差0.1890.00010.00030.00010.00010.0505代數(shù)精度1233105牛頓-柯特斯方法是一種利用插值多項(xiàng)式來(lái)構(gòu)造數(shù)值積分的常用方法，這其中梯形積分方法的誤差最大，近似效果最差，辛普森方法的精度比梯形積分高了一個(gè)數(shù)量級(jí)，它的代數(shù)精度比梯形積分的代數(shù)精度高，能更好地近似積分值；牛頓-柯特斯積分方法的誤差比辛普森積分精度高兩個(gè)數(shù)量級(jí)復(fù)合梯形積分方法比單獨(dú)的梯形積分精度高，龍貝格方法收斂速度快、計(jì)算精度較高，但是計(jì)算量較大高斯求積方法積分精度高、數(shù)值穩(wěn)定、收斂速度較快，但是節(jié)點(diǎn)與系數(shù)的計(jì)算較麻煩，而且要求已知積分函數(shù)結(jié) 論本文主要討論了數(shù)值積分的計(jì)算方法并通過(guò)MATLAB軟件編程實(shí)現(xiàn)，通過(guò)前面的研究我們知道求數(shù)值積分近似值的計(jì)算