數值積分用matlab實現

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

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

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

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

5、公式就是辛普森求積公式： （2.4）2.2 牛頓柯特斯公式的Matlab實現functionC, g = NCotes(a, b, n, m)% a，b分別為積分的上下限；% n是子區(qū)間的個數； % m是被調用第幾個被積函數； % 當n=1時計算梯形公式；當n=2時計算辛普森公式，以此類推；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 復合求積公式及Matlab實現3.1 復合梯形公式的基本原理將區(qū)間a, b劃分成n等分

7、，分點，,在每個子區(qū)間（）上采用梯形公式得： （3.1）記 （3.2）稱式(3.2)為復合梯形公式3.2 復合梯形公式的Matlab實現function s = trapr1(f, a, b, n)% f表示被積函數；% a，b表示積分上下限；% n是子區(qū)間的個數；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 復合辛普森求積公式的基本原理將區(qū)間a,b分等分，在每一個子區(qū)間上采用辛普森公式，若記，則得 （3.3）記（3.4）稱式(3.4)為復合辛普森求積公式3.4 復合辛普森求積公式的Matlab實現function s = simpr1(f, a, b, n) % f表示被積函數；% a，b表示積分上下限；% n是子區(qū)間的個數；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實現4.1 龍貝格算法的基本原理由梯形的遞推法可以看出，將積分區(qū)間等分時，用復化梯形公式計算的結果作為積分的近似值，其誤差近似值為 （4.1）可以設想，如果用這個誤差作為的一種補償，即將 （4.2）作為積分的近似值，可望提高其精確度直接根據復化求積公式，不難驗證 （4.3）這說明，將區(qū)間對分前后兩

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

11、區(qū)間分半，算出及的值，并根據(4.3)和(4.9)式計算和，再有公式(4.4)求出；第四步：再將區(qū)間分半，計算，, ,并根據公式(4.5)計算第五步：再將區(qū)間分半，類似上述過程計算，, ，重復以上步驟即可得到， ,一直到龍貝格序列中前后兩項的絕對值差不超過給定的誤差險為止4.2 龍貝格算法的Matlab實現function R, quad, err, h=romber(f, a, b, n, delta) % f表示被積函數；% a，b表示積分上下限；% n是子區(qū)間的個數% 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實現5.1 高斯勒讓德求積公式的基本原理在高斯求積公式 （5.1）中，若取權函數=1，區(qū)間-1,1，則得公式 （5.2）我們知道勒讓德多項式是區(qū)間-1,1上的正交多項式，因此，勒讓德多項式的零點就是求積公式(5.2)的高斯點，形如(5.1)式的高斯公式特別的稱為高斯-勒讓德求積公式如下表5.1所示為高斯-勒讓德求積公式的節(jié)點數和系數5.2 高斯勒讓

14、德求積公式的Matlab實現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é)點數和系數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 各個求積公式的應用舉例與比較分析6.1 簡單數值積分的解(精確值0.946083070367183)6.1.1 牛頓柯特斯當n=1時的梯形算法和n=2時的辛普森算法的結果解：先用M文件定義一個f1.m的函數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 復合梯形公式和復合辛普森求積公式的計算結果解：建立一個M文件定義一個f.m函數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、龍貝格求積公式的計算結果（取誤差不超過）解：建立一個M文件定義一個f.m函數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 高斯勒讓德求積公式的計

18、算結果(給定節(jié)點3)解：建立一個M文件定義一個f.m函數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 復雜數值積分的解(精確值-1.8785)6.2.1 牛頓柯特斯當n = 1時的梯形算法和n = 2時的辛普森算法的結果解：先用M文件定義一個f1.m的函數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 復合梯形公式和復合辛普森求積公式的計算結果解：先用M文件定義一個f1.m的函數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 龍貝格求積公式的計算結果（取誤差不超過）解：先用M文件定義一個f1.m的函數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、 各個求積公式的比較分析（以的各個積分結果為例） 表6.1各個數值積分的比較積分方法牛頓柯特斯（梯形）牛頓-柯特斯（辛普森）復合梯形復合辛普森龍貝格高斯-勒讓德精確值0.94610.94610.94610.94610.94610.9461實際值0.92700.94600.94580.94600.94600.8956誤差0.1890.00010.00030.00010.00010.0505代數精度1233105牛頓-柯特斯方法是一種利用插值多項式來構造數值積分的常用方法，這其中梯形積分方法的誤差最大，近似效果最差，辛普森方法的精度比梯形積分高了一個數量級，它的代數精度比梯形積分的代數精度高，能更好地近似積分值；牛頓-柯特斯積分方法的誤差比辛普森積分精度高兩個數量級復合梯形積分方法比單獨的梯形積分精度高，龍貝格方法收斂速度快、計算精度較高，但是計算量較大高斯求積方法積分精度高、數值穩(wěn)定、收斂速度較快，但是節(jié)點與系數的計算較麻煩，而且要求已知積分函數結 論本文主要討論了數值積分的計算方法并通過MATLAB軟件編程實現，通過前面的研究我們知道求數值積分近似值的計算