數(shù)值分析課程設(shè)計含代碼

1、下載可編輯成績評定表學生班級學號專業(yè)信息與計算課程設(shè)計題目數(shù)值分析算法案科學例評語組長簽字：成績?nèi)掌?0年月日.專業(yè) .整理 .下載可編輯課程設(shè)計任務(wù)書學院理學院專 業(yè)信息與計算科學學生班級學號課程設(shè)計題目數(shù)值分析算法案例實踐教學要求與任務(wù) :要求：格式以學校畢業(yè)論文格式要求為準，不準粘貼圖片，尤其公式。對每個試驗，要求有：實驗基本原理，實驗?zāi)康?，實驗容及?shù)據(jù)來源和實驗結(jié)論。以班級為單位統(tǒng)一裝訂封皮。6 月 25 日，十八周的周二交論文每人至少四個實驗，最少15 頁任務(wù)（實驗項目） ：線性方程組數(shù)值解法參考題目：（ 1)列主元 Gauss 消去法（ 2） LU分解法插值法和數(shù)據(jù)擬合參考題目：（

2、 1） Lagrange 插值（ 2） Newton 插值（ 3）最小二乘擬合數(shù)值積分參考題目：（ 1）復化 Simpon 積分（ 2）變步長的梯形積分公式（3）龍貝格求積公式常微分方程數(shù)值解Runge-Kutta方法數(shù)值方法實際應(yīng)用用數(shù)值方法解決實際問題（自選）工作計劃與進度安排:線性方程組數(shù)值解法（4 學時）插值法和數(shù)據(jù)擬合（4 學時）數(shù)值積分（4 學時）常微分方程數(shù)值解（4 學時）數(shù)值方法實際應(yīng)用（4 學時）答辯（4 學時）指導教師：專業(yè)負責人：學院教學副院長：201年月日201年月日201年月日.專業(yè) .整理 .下載可編輯摘 要實驗方法與理論方法是推動科學技術(shù)發(fā)展的兩大基本方法， 但有

3、局限性。許多研究對象， 由于空間或時間的限制， 既不可能用理論精確描述， 也不能用實驗手段實現(xiàn)。數(shù)值模擬或稱為科學計算突破了實驗和理論科學的局限， 在科技發(fā)展中起到越來越重要的作用。 可以認為，科學計算已于實驗、 理論一起成為科學方法上不可或缺的三個主要手段。計算數(shù)學的研究是科學計算的主要組成部分，而數(shù)值分析則是計算數(shù)學的核心。數(shù)值計算是研究使用計算機來解決各種數(shù)學問題的近似計算方法與理論， 其任務(wù)是提供在計算機上可解的、理論可靠的、計算復雜性低的各種常用算法。數(shù)值分析的主要容：1）、數(shù)值代數(shù)：求解線性和非線性方程組的解， 分直接方法和間接方法兩大類；2）、插值、曲線擬合和數(shù)值逼近；3）、數(shù)值

4、微分和數(shù)值積分；4）、常微分和偏微分方程數(shù)值解法。本文主要通過 Matlab 軟件，對數(shù)值分析中的一些問題進行求解，如列主元Gauss消去法， Lagrange 插值多項式，復化 Simpson 公式， Runge-Kutta 方法以及數(shù)值分析在實際問題中的應(yīng)用，并在求解的過程中更加熟識這門課程的主要容，以及加強對課程知識的掌握。在學習與設(shè)計計算方法時，從數(shù)學理論角度，學會分析方法的誤差、 收斂性和穩(wěn)定性， 保證計算方法的準確性； 從實際應(yīng)用的角度出發(fā)，掌握計算方法的結(jié)構(gòu)與流程， 能夠把計算方法轉(zhuǎn)換為可在計算機上直接處理的程序，保證算法的可用性。關(guān)鍵詞：列主元 Gauss 消去法；Lagran

5、ge 插值；復化 Simpson 公式；Runge-Kutta目 錄.專業(yè) .整理 .下載可編輯實驗一列主元 Gauss 消去法11.1實驗?zāi)康?1.2基本原理11.3實驗容21.4實驗結(jié)論3實驗二拉格朗日插值多項式42.1實驗?zāi)康?2.2基本原理42.3實驗容42.4實驗結(jié)論8實驗三復化 Simpson 求積公式93.1實驗?zāi)康?3.2基本原理103.3實驗容103.4實驗結(jié)論12實驗四龍格 - 庫塔 (Runge-Kutta) 方法124.1實驗?zāi)康?24.2基本原理124.3實驗容134.4實驗結(jié)論15實驗五數(shù)值方法實際應(yīng)用155.1實驗?zāi)康?55.2基本原理155.3實驗容165.4實

6、驗結(jié)論22參考文獻22.專業(yè) .整理 .下載可編輯實驗一列主元 Gauss 消去法1.1 實驗?zāi)康?) 理解列主元消去法的原理；2) 熟悉列主元消去法的計算步驟，能用代碼編寫；3) 解決實際問題。1.2 基本原理在順序 Gauss 消去法中，必須要求 a(k)0(k1,2, , ) ；否則無法進行計算。kkn即使 akk(k )0 ，但其絕對值 akk(k) 很小，由于舍入誤差的影響，也可能會引起很大的誤差，從而使上述方法失效。為了使消元過程中減小舍入誤差和不至于中斷，可以按照不同的自然順序進行消元。在第k 步消元時，增廣矩陣為a11(1)a12(1)a1(1k)1a22(2)a2(k2) 1

7、A( k )B( k )( k 1)ak 1k 1a1(1k)a1(1n)a2(2k)a2(2n)( k1)( k1)a k1kak1nakk(k )akn(k )ank(k )ann(k )b1(1)b2(2)( k 1)（1.1 ）bk 1bk(k)bn(k)(k1)( k 1)(k1)(k 1)中選取絕對值最大的作不一定選取 akk作為主元，而從同列 akk, ak1,k , ank為主元素，即ar(kk 1)max aik(k1)（1.2 ）ki n若 ark(k )0，此時矩陣不可逆，方程的解不確定，則停止計算；否則，當 r>k時，則其增廣矩陣中交換第k 行和第 r 行，即ak

8、j(k)arj(k ) jk, k 1, nbk(k )br( k)（1.3 ）使 ark(k ) 成為主元。然后再按 Gauss消去法進行消元運算。 于是就得到列主元 Gauss 消去法。.專業(yè) .整理 .下載可編輯1.3 實驗容程序來源首先建立一個 gaussMethod.m 的文件，用來實現(xiàn)列主元的消去方法。 文件容如下：function x=gaussMethod(A,b)%高斯列主元消去法，要求系數(shù)矩陣非奇異的n = size(A,1);if abs(det(A)<= 1e-8error('系數(shù)矩陣是奇異的 '); return;endfor k=1:nak =

9、 max(abs(A(k:n,k);index = find(A(:,k)=ak);if length(index) = 0index = find(A(:,k)=-ak);end%交換列主元temp = A(index,:);A(index,:) = A(k,:);A(k,:) = temp;temp = b(index);b(index) = b(k); b(k) = temp; %消元過程for i=k+1:nm=A(i,k)/A(k,k);%消除列元素A(i,k+1:n)=A(i,k+1:n)-m*A(k,k+1:n);b(i)=b(i)-m*b(k);endend% 回代過程x(n)

10、=b(n)/A(n,n);for k=n-1:-1:1;x(k)=(b(k)-A(k,k+1:n)*x(k+1:n)')/A(k,k);.專業(yè) .整理 .下載可編輯end; end然后調(diào)用gaussMethod 函數(shù)，來實現(xiàn)列主元的高斯消去法。建立一個文件gauss，容如下：cleardisp('*')x=gaussMethod(input('請輸入系數(shù)矩陣： '),input('請輸入常數(shù)列：')disp('*')實例分析例：在 Matlab 上，利用列主元法求線性方程組的解:x12x2x34x4132x10x24x33

11、x4284x12x22x3x4203x1x23x32x46解：運行程序，按照提示輸入方程的系數(shù)矩陣及常數(shù)列，如下所示：*請輸入系數(shù)矩陣： 1 2 1 4;2 0 4 3;4 2 2 1;-3 1 3 2請輸入常數(shù)列： 13;28;20;6x =3-142*即該方程的解為：x3, 1,4,21.4 實驗結(jié)論把向量計算得到的解帶入方程組，驗證正確性， 和其他的方法比較， 列主元具有一定的簡單性， 比較容易實現(xiàn)。 避免使用其他方法的誤差或不能進行性。而列主元也有一定的限制，要求行列式的值不為0。.專業(yè) .整理 .下載可編輯實驗二拉格朗日插值多項式2.1 實驗?zāi)康?）熟悉簡單的拉格朗日插值多項式的基本

12、概念；2）熟悉 Lagrange 公式及源代碼，會利用它來計算基本函數(shù)；3）能構(gòu)造出正確的插值多項式；2.2 基本原理設(shè)函數(shù) yf ( x) 在區(qū)間 a,b上有定義，且已知在點 ax0 x1 xxn b上的值 y0 , y1, yn, 若存在一個次數(shù)不超過n 的多項式Ln (x)a0a1xan xn（2.1 ）使其滿足Ln ( xk )yk ,k0,1, ,n（2.2 ）則稱 Ln (x) 為 f ( x) 的 n次插值多項式，稱點 x k (k 0,1,n) 為插值節(jié)點，稱條件（ 2.2 ）為插值條件。包含插值節(jié)點的區(qū)間成為插值區(qū)間。通過平面上不同的兩點可以確定一條直線經(jīng)過這兩點， 就是拉格

13、朗日線性插值問題，對于不在同一直線的三點得到的插值多項式為拋物線。 拉格朗日是比較基礎(chǔ)的方法，本身比較容易實現(xiàn)，容易理解。給定 n+1 個不同節(jié)點，構(gòu)造f x0 , x1, x2 , xi 的 n 次拉格朗日插值多項式：L ( x)nn 1(xx j )（2.3 ）yi li ( x) ， l i ( x)( xi,i 1,2, , n 1i 0j1xj )ji2.3 實驗容程序來源首先建立一個 Lagrange.m 的文件，用來實現(xiàn)Lagrange 插值。文件容如下：%輸入： x 是插值節(jié)點橫坐標向量；y 是插值節(jié)點對應(yīng)縱坐標向量%輸出： C是拉格朗日插值多項式的系數(shù)矩陣；L 是插值函數(shù)系數(shù)

14、矩陣functionC,L=Lagrange(x,y).專業(yè) .整理 .下載可編輯w=length(x);n=w-1;L=zeros(w,w);for k=1:n+1V=1;for j=1:n+1if k=jV=conv(V,poly(x(j)/(x(k)-x(j);endendL(k,:)=V;endC=y*L然后調(diào)用 Lagrange 函數(shù)，來實現(xiàn)Lagrange 插值法。建立一個文件Lg，容如下：cleardisp('*')x=input('請輸入已知點的橫坐標組：');y=input('請輸入已知點的縱坐標：');C,L=Lagrange

15、(x,y);yi=polyval(C,input('請輸入需要計算得橫坐標組：')xx=1.5:0.05:6.5;yy=polyval(C,xx);plot(xx,yy,x,y,'o')disp('*')實例分析例 1 有 4 對數(shù)據(jù)（ 1.6,3.3 ），（2.7,4.22 ），（3.9,5.61 ），（ 5.6,2.94 ），寫出這 4個數(shù)據(jù)點的 Lagrange 插值公式，并計算出橫坐標組 xi=2.101,4.234 時對應(yīng)的縱坐標值。解： 4 個數(shù)據(jù)點的 Lagrange 插值公式為：.專業(yè) .整理 .下載可編輯（x）3.3*( x2.

16、7) * (x3.9) * (x5.6)( x1.6) * ( x3.9) * ( x 5.6)L 3(1.62.7) * (1.63.9) * (1.64.22*(2.71.6) * (2.73.9) * (2.7 5.6)5.6)3.9 *(x1.6) * ( x2.7) * (x 5.6)( x1.6) * (x2.7) * ( x 3.9)1.6) * (3.92.7) * (3.9 5.6)2.94 *1.6)* (5.62.7) * (1.6 3.9)(3.9(5.6運行程序，按照提示輸入已知點的橫坐標組、 縱坐標組及需要計算得橫坐標組，如下所示：*請輸入已知點的橫坐標組：1.6,

17、2.7,3.9,5.6請輸入已知點的縱坐標：3.3,1.22,5.61,2.94C =-1.053911.0551 -34.493334.5053請輸入需要計算得橫坐標組：2.101,4.234yi =1.05966.6457*即y10539x3110551234 4933x345053- .x- .輸出圖形：1050-5-10-1522.533.544.555.566.51.5圖 2.1輸出擬合曲線例 2 將區(qū)間 -5 ， 5 等分 5 份、 10 份，求函數(shù) y1拉格朗日差值多項式，x21做出函數(shù)原圖像，觀察龍格現(xiàn)象。解：首先將區(qū)間五等分，取各端點坐標擬合曲線，輸入：.專業(yè) .整理 .下載

18、可編輯clearx=-5:2:5;y=1./(1+x.2);C,L=Lagrange(x,y);輸出結(jié)果：C =0.0000 0.0019 -0.0000-0.06920.00000.5673即y0.0019x40.0692x20 5673.然后輸出圖形，比較擬合效果。在matlab 中輸入：xx=-5:0.1:5;yy=polyval(C,xx);hold onplot(xx,yy,'-',x,y,'o')xp=-5:0.01:5;z=1./(1+xp.2);plot(xp,z,'r')legend(' 拉格朗日插值曲線 ',&

19、#39;插值點 ','原曲線 ')輸出五等分插值龍格現(xiàn)象圖形：1.2拉格朗日插值曲線插 值 點1原 曲 線0.80.60.40.20-0.2-5-4-3-2-1012345圖 2.2五等分插值龍格現(xiàn)象圖形十等分的過程與五等分基本相似，如下輸入：clear,clfx=-5:1:5;y=1./(1+x.2);.專業(yè) .整理 .下載可編輯C,L=Lagrange(x,y);輸出：C =-0.0000 -0.00000.0013 -0.0000-0.0244 0.00000.1974-0.0000 -0.6742 -0.00001.0000即y0.0013x80.0244x60

20、.1974x40.6742x 2 1然后輸出圖形，比較擬合效果。在matlab 中輸入：xx=-5:0.1:5;yy=polyval(C,xx);hold onplot(xx,yy,'-',x,y,'o')xp=-5:0.01:5;z=1./(1+xp.2);plot(xp,z,'r')legend('拉格朗日插值曲線 ','插值點 ','原曲線 ')輸出十等分插值圖形：2拉格朗日插值曲線插 值 點原 曲 線1.510.50-0.5-5-4-3-2-1012345圖 2.3十等分插值龍格現(xiàn)象圖形2.

21、4 實驗結(jié)論通過圖像可以得出：1）并不是插值節(jié)點越多，插值多項式逼近函數(shù)效果就越好；2）誤差較大的地方，是在插值區(qū)間兩端點附近出現(xiàn)；.專業(yè) .整理 .下載可編輯3）求插值多項式，不需要求解線性方程組，當數(shù)據(jù)比較多時，此公式才能顯示出優(yōu)越性；4）函數(shù)值可以用符號形式表示， 數(shù)據(jù)點未確定的縱坐標可以用多項式表示；實驗三復化 Simpson 求積公式3.1 實驗?zāi)康?）了解單獨的 Simpson 公式原理及幾何意義；2）學會復化 Simpson 求積公式的編程與應(yīng)用；3）掌握 Mat lab提供的計算積分的各種函數(shù)的使用方法；.專業(yè) .整理 .下載可編輯3.2 基本原理復化 Simpson 公式是一

22、種比較實用的積分方法，可以給出誤差估計。 首先將區(qū)間 a,b N等分，子區(qū)間的長度為ba（3.1）hnN在每個子區(qū)間上采用 Simpson 公式，在用 Simpson 公式時，還要將子區(qū)間再二等分，因此有 2N+1個分點。即xkx0k hN , k0,1,2N , x0a.（3.2 ）2經(jīng)推導得到，defhN ()( ) 2N 1()4N()bf ( x)dxffx2 kfx2 k 1（3.3）SNa6f abk 1k 1稱為 SN 為復化 Simpson 值，稱（ 3.3 ）式為復化 Simpson 公式。3.3 實驗容程序來源編寫復化 Simpson 求積函數(shù)（函數(shù)名： s_quad.m）

23、Function I=S_quad(x,y)% 復化求積公式% x 為被積函數(shù)自變量的等距節(jié)點； y 為被積函數(shù)在節(jié)點處的函數(shù)值。n=length(x);m=length(y); % 積分自變量的節(jié)點數(shù)應(yīng)與它的函數(shù)值個數(shù)相同； if n=merror ('The length of X and Y must be equal');return;endif rem(n-1,2)=0%如果 n-1 不能被 2 整除，則調(diào)用復化公式error ('節(jié)點數(shù)不滿足要求 ');return;end.專業(yè) .整理 .下載可編輯N=(n-1)/2;h=(x(n)-x(1)/N;

24、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ù)，來實現(xiàn)復化Simpson 公式法。建立一個文件SPS，容如下：cleardisp('*')x=input('請輸入積分上下限及點間的間隔（例如-1:0.1:1）：');y=input('請輸入被積公式： y=');I=S_quad(x,y);disp('得出積分值 I=')disp(I);disp('

25、*')實例分析：1 2例 1 用復化 Simpson公式求積分 e x dx ，在積分區(qū)間中點與點之間的間隔取為10.1。解：運行程序，按照提示輸入積分上下限、點間的間隔及被積公式，如下所示：*請輸入積分上下限及點間的間隔（例如-1:0.1:1）：-1:0.1:1請輸入被積公式： y=exp(-x.2)得出積分值 I=1.4936*.專業(yè) .整理 .下載可編輯真值為： 1.49371 x例 2 計算積分 0 4 2 dx ，將區(qū)間 8 等分。x解：運行程序，按照提示輸入積分上下限、等分后的區(qū)間長度及被積公式，如下所示：*請輸入積分上下限及點間的間隔（例如-1:0.1:1）：0:0.12

26、5:1請輸入被積公式： y=x./(4+x.2)得出積分值 I=0.1116*真值為： 0.1115723.4 實驗結(jié)論復化 Simpson 計算所得的結(jié)果誤差較小， 精度較高，更適合科學計算與應(yīng)用，且公式具有收斂性，穩(wěn)定性良好。實驗四龍格 - 庫塔 (Runge-Kutta)方法4.1 實驗?zāi)康?）熟悉 Runge-Kutta 常微分方程初值問題的基本原理2）了解 Runge-Kutta 常微分方程初值問題的計算3）能編程實現(xiàn) Runge-Kutta 常微分方程初值問題4.2 基本原理龍格 - 庫塔 (Runge-Kutta) 方法是一種在工程上應(yīng)用廣泛的高精度單步算法。由于此算法精度高，

27、采取措施對誤差進行抑制， 所以其實現(xiàn)原理也較復雜。該算法是構(gòu)建在數(shù)學支持的基礎(chǔ)之上的。對于一階精度的歐拉公式有：.專業(yè) .整理 .下載可編輯yn 1ynh * K1(4.1)K1f xn , yn當用點 xn 處的斜率近似值 K1 與右端點 xn 1 處的斜率 K 2的算術(shù)平均值作為平均斜率 K * 的近似值，那么就會得到二階精度的改進歐拉公式：yn 1ynh * K1 K 2 / 2K1f xn , yn(4.2)K 2f xn 1, yn h* K1依次類推，如果在區(qū)間xn , xn 1 多預(yù)估幾個點上的斜率值，并用他們的加權(quán)平均數(shù)作為平均斜率K * 的近似值，顯然能構(gòu)造出具有很高精度的高

28、階計算公式。經(jīng)數(shù)學推導、 求解，可以得四階龍格庫塔公式，也就是在工程中應(yīng)用廣泛的經(jīng)典龍格庫塔算法。四階龍格庫塔方法基本公式：yk 1 ykh2K 22 K 3K 4 )( K16K 1f ( tk , yk )K 2f (tkh , ykh K1 )(4.3)22K 3f ( tkh , ykh K2 )22K 4f (t kh, ykhK 3 )4.3 實驗容程序來源寫出四階 Runge-Kutta 方法源代碼的f (t, y) 得.m 文件：function R=rk4(f,a,b,ya,M)% a,b 分別為左右端點， ya 為初始值， M為步長幅數(shù)。h=(b-a)/M;T=zeros(

29、1,M+1);Y=zeros(1,M+1);T=a:h:b;Y(1)=ya;for j=1:Mk1=feval(f,T(j),Y(j);.專業(yè) .整理 .下載可編輯k2=feval(f,T(j)+h/2,Y(j)+k1*h/2);k3=feval(f,T(j)+h/2,Y(j)+k2*h/2);k4=feval(f,T(j)+h,Y(j)+h*k3);Y(j+1)=Y(j)+(k1+2*k2+2*k3+k4)/6;endR=T' Y'4.3.2 實例分析：例 用四階龍格庫塔方法求解初值問題y 't 2ty 取步長h=0.1 。 并將y(0)0a,b,ya,M 分別取定

30、0,1,0,10。解：首先寫出執(zhí)行函數(shù)的 tt.m文件，其中 tt表示 f ，即為求導函數(shù)：function f=tt(t,y)f=t2+t-y;再輸入調(diào)用四階Runge-Kutta 方法，執(zhí)行語句如下：>> R=rk4('tt',0,1,0,10)輸出：R =000.100000000000000.051627083333330.200000000000000.168478293489580.300000000000000.307517220780890.400000000000000.466661978888610.500000000000000.6458118

31、56562070.600000000000000.844961981894520.700000000000001.064112119207480.800000000000001.303262257100000.900000000000001.562412395020551.000000000000001.84156253294245.專業(yè) .整理 .下載可編輯4.4 實驗結(jié)論龍格 - 庫塔法具有精度高，收斂，穩(wěn)定（在一定條件下） ，計算過程中可以改變步長，不需要計算高階導數(shù)等優(yōu)點，但仍需計算在一些點上的值，如四階龍格- 庫塔方法每計算一步需要計算四次的值，這給實際計算帶來一定的復雜性。實驗五數(shù)

32、值方法實際應(yīng)用5.1 實驗?zāi)康?）了解梯形公式原理及幾何意義；2）學會復化梯形求積公式的編程與應(yīng)用；3）掌握 Matlab 提供的計算積分的各種函數(shù)的使用方法；5.2 基本原理在實際問題中， 往往會遇到一些困難。 有些函數(shù)找不到用初等函數(shù)表示的原函數(shù)，例如，對于積分1 sinxdx(5.1)0 x而言，不存在用初等函數(shù)表示的原函數(shù)。 而有些函數(shù)雖然能找到原函數(shù)， 但計算過于復雜，例如，橢圓型積分x2bx cdx(5.2)ax2x1而有些情況下， 只能知道某些點處的函數(shù)值， 并沒有函數(shù)的具體表達式。 這些情況，使我們有必要研究積分的數(shù)值計算問題。下面我們就以梯形公式為例做.專業(yè) .整理 .下載可

33、編輯以說明。所謂梯形求積公式就是用梯形面積來近似曲邊梯形面積， 利用梯形公式和連續(xù)增加 a,b 的區(qū)間數(shù)來逼近：bf (x)dxh2j f (xk 1) f (xk )a2 k1(5.3)第 j 次循環(huán)在 2 j1個等距節(jié)點處對f x 采樣。5.3 實驗容衛(wèi)星軌道是一個橢圓，橢圓周長計算公式是2（c 2sin2d ，）S 4a1-0 a這里 a 是橢圓半長軸， c 是地球中心與軌道中心（橢圓中心）的距離，記h 為近地點距離， H為遠地點距離， R=637km為地球半徑，則a(2RHh) / 2,c( Hh) / 2,我國第一顆人造衛(wèi)星近地點距離h=439km，遠地點距離 H=2384，試求衛(wèi)星

34、軌道的周長。解 : 第一步：先利用 Matlab 畫出被積函數(shù)的圖形。輸入程序如下：clear H=2384; h=439; R=6371; a=(2*R+H+h)/2 c=(H-h)/2 x=0:0.1:pi/2;y=sqrt(1-(c/a)2*(sin(x).2); plot(x,y,'-')title(' 梯形法則 '); xlabel('x');.專業(yè) .整理 .下載可編輯ylabel('y');輸出結(jié)果：a =7.782500000000000e+003c =9.725000000000000e+002輸出圖形：梯形法則

35、1.00110.9990.9980.997y0.9960.9950.9940.9930.99200.511.5x圖 5.1被積函數(shù)的圖形第二步：應(yīng)用數(shù)值積分梯形公式。首先建立一個名為trapezg.m 的 M文件，程序如下：function I=trapezg(f_name3,a,b,n)format long%輸出用 15 位數(shù)字表示.專業(yè) .整理 .下載可編輯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:10

36、0)*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然后建立一個名為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;輸入命令程序：>> trapezg('f_name3',0,pi/2,30)輸出結(jié)果：ans =0.00955791054630輸出圖形：.專業(yè) .整理 .下載可編輯數(shù)值積分梯形效果圖0.0120.010.008y 0.0060.0040.002000.20.40.60.811.21.41.6x圖 5.2數(shù)值積分效果