實驗目的：
1.熟悉Matlab矩陣操作
2.用Matlab實現求積公式：梯形，辛普森，復合梯形，復合辛普森等算法
3.學會數值積分有關應用
實驗要求：
1.掌握梯形和辛普森算法
2.掌握復合梯形和復合辛普森算法
3.掌握判斷迭代停止的手段
實驗內容：
1.矩陣的四則運算、冪運算；矩陣元素的查找、排序、求和、求積差分。
2.用Matlab實現求積公式：梯形，辛普森，復合梯形，復合辛普森等算法
3.學會數值積分有關應用
（補充2，不在實驗過程中展示）
實驗步驟：
1.矩陣的四則運算，進行加減法的前提是參與運算的兩個矩陣或多個矩陣必須具有相同的行數和列數；或者其中有一個或多個矩陣為標量。其加減法定義：
，當其中含有標量x時，
，且矩陣的加法運算滿足“交換律”、“結合律”、“存在零元”、“存在負元”的運算律。矩陣加減法示例：
數與矩陣到乘法、矩陣與矩陣的乘法、矩陣的除法。
C的元素就是用書x乘矩陣A對應的元素而得到，數與矩陣的乘法滿足下列運算律：
， 
， 
， 
數乘示例：
A為m*h矩陣，B為h*n矩陣，則兩矩陣的乘積
為一個矩陣，且
且滿足“結合律”、“左分配律”、“右分配律”、“單位矩陣的存在性”。矩陣乘法示例：
矩陣的除法是乘法的逆運算，左除“\”和右除“/”。若A和B是標量，A/B和B\A等價。對一般的二維矩陣A和B，當進行A\B運算時，要求A的行數與B的行數相等；當進行A/B運算時，要求A到列數與B的列數相等。矩陣除法示例：
矩陣的冪運算，當矩陣A為方陣是可進行矩陣的冪運算，其定義為：
方陣冪運算示例：
矩陣元素的查找，函數find()的作用是進行矩陣元素的查找，它通常與干洗函數和邏輯運算相結合。語法格式：ind=find(X):該函數查找矩陣X中的非零元素，返回這些元素的單下標；[row,col]=find(X,...):該函數查找矩陣X中的非零元素，函數返回這些元素的雙下標i和j。
矩陣元素的排序，函數sort()的作用是按升序排序，排序后的矩陣和眼函數矩陣位數一致，語法：B=sort(A):該函數對矩陣A進行升序排序，A可為矩陣或向量；B=sort(A,dim)：當dim=1時矩陣A按列升序，當dim=2時矩陣按行升序；B=sort(...,mode):該函數按mode指定的方式排序（ascend升序，descend降序）。
矩陣元素的求和、求積和差分。函數sum()和cumsum()的作用是對矩陣元素求和，函數prod()和cumprod()的作用是對矩陣元素求積，函數diff()的作用計算矩陣的差分示例：
2.求積公式，梯形公式示例：
梯形求積代碼：
1 function t=rctrap(fun,a,b,n) 2 %復化梯形公式 3 %n等分 4 %a,b區間的左右端點 5 h=(b-a)/n; 6 t=0; 7 for i=1:n-1 8 t=t+h*feval(fun,i*h+a); 9 end 10 t=t+0.5*h*(feval(fun,a+eps)+feval(fun,b))
運行部分截圖：
辛普森求積代碼：
1 function s=sptrap(fun,a,b,n) 2 % n，對應的等分點為2n 3 h=(b-a)/(2*n); 4 s=0; 5 for i=1:n 6 s=s+feval(fun,(2*i-1)*h+a)*4; 7 end 8 for i=1:n-1 9 s=s+feval(fun,(2*i)*h+a)*2; 10 end 11 s=s+feval(fun,a+eps)+feval(fun,b); 12 s=s*h/3
運行部分截圖：
根據梯形公式和估計誤差公式編寫程序計算積分，分別取
，用復合梯形公式計算定積分
並與精確值比較.然后觀察 對計算結果的有效數字和絕對誤差的影響.
1 h=pi/8000;a=0;b=pi/2;x=a:h:b;n=length(x), y=exp(sin(x)); 2 z1=(y(1)+y(n))*h/2; z2=sum(y(2:n-1))*h; z8000=z1+z2, 3 syms t 4 f=exp(sin(t)); intf=int(f,t,a,b), Fs=double(intf), 5 Juewucha8000=abs(z8000-Fs)
運行部分示例：
用復合梯形公式計算定積分
，並與精確值比較.然后觀察 對計算結果的有效數字和絕對誤差的影響.
代碼：
1 function T=rctrap1(fun1,a,b,m) 2 n=1;h=b-a; T=zeros(1,m+1); x=a; T(1)=h*(feval(fun1,a)+feval(fun1,b))/2; 3 for i=1:m 4 h=h/2; n=2*n; s=0; 5 for k=1:n/2 6 x=a+h*(2*k-1); s=s+feval(fun1,x); 7 end 8 T(i+1)=T(i)/2+h*s; 9 end 10 T=T(1:m);
運行：
用復合辛普森公式計算定積分
取n=20001個等距節點，並與精確值比較.然后再取n=13，觀察 n對誤差的影響.解 由n=2m+1=20001,得m=10000.根據辛普森公式編寫並輸入下面的程序。代碼：
1 a=0;b=1;m=10000 2 m=6; h=(b-a)/(2*m); x=a:h:b; 3 y=exp((-x.^2)./2)./(sqrt(2*pi)); 4 z1=y(1)+y(2*m+1); z2=2*sum(y(2:2:2*m)); 5 z3=4*sum(y(3:2:2*m)); 6 z=(z1+z2+z3)*h/3, syms t,f=exp((-t^2)/2)/(sqrt(2*pi)); 7 intf=int(f,t,a,b), Fs=double(intf), Juewucha=abs(z-Fs)
運行：
復合辛普森數值積分的MATLAB程序comsimpson（fun,a,b,n)。用comsimpson.m和quad.m分別計算定積分
取精度為10e-4，並與精確值比較.。代碼：
1 [Q1,FCNT14] = quad(@fun2,0,1,1.e-4,3), 2 Q2 =comsimpson (@fun2,0,1,10000) 3 syms x 4 fi=int(exp( (-x^2)/2)/(sqrt(2*pi)),x,0, 1); 5 fs=double(fi) 6 wQ1= double (abs(fi-Q1)), wQ2= double (abs(fi-Q2))
運行：
3.數值分析的應用。利用復合梯形公式求
。代碼：
1 function I=tquad(x,y) 2 n=length(x);m=length(y); 3 if n~=m 4 error('向量x,y的長度必須一致'); 5 end 6 h=(x(n)-x(1))/(n-1); 7 a=[1 2*ones(1,n-2) 1]; 8 I=h/2*sum(a.*y);
運行：
小結：
在參考PPT和書本的MATLAB編程之后，感覺解決問題使用計算機的方法，就好像使用一種特殊的語言與一堆無機物交流，並通過一些物理的方法使它回應我們。解數學題的編程的關鍵，就是熟悉數學規律，並使用該編程語言簡潔地表達出來。