1、第第3章章 MATLAB在高等數(shù)學中的應用在高等數(shù)學中的應用 3.1 矩陣分析矩陣分析3.2 多項式運算多項式運算 3.3 數(shù)據(jù)的分析與統(tǒng)計數(shù)據(jù)的分析與統(tǒng)計 3.4 函數(shù)分析與數(shù)值積分函數(shù)分析與數(shù)值積分 3.1 矩陣分析矩陣分析 1 1矢量范數(shù)和矩陣范數(shù)矢量范數(shù)和矩陣范數(shù) 矩陣范數(shù)是對矩陣的一種測度。矢量的p范數(shù)和矩陣A的p范數(shù)分別定為： 當p2時為常用的歐拉范數(shù)，一般p還可取l和。這在MATLAB中可利用norm函數(shù)實現(xiàn)，p缺省時為p=2。格式：格式：nnorm(A)功能：計算矩陣A的最大奇異值，相當于n=max(svd(A)。格式：格式：nnorm(A，p)功能：norm函數(shù)可計算幾種不同

2、類型的矩陣范數(shù)，根據(jù)p的不同可得到不同的范數(shù) pipipxxpipipxxppxpxAxAmax2 2矩陣求逆及行列式值矩陣求逆及行列式值 矩陣求逆函數(shù)矩陣求逆函數(shù)inv及行列式值函數(shù)及行列式值函數(shù)det逆矩陣的定義：逆矩陣的定義：對于任意階 nn 方陣A，如果能找到一個同階的方陣V，使得滿足：A*V=I。其中I為n階的單位矩陣eye(n)。則V就是A的逆矩陣。數(shù)學符號表示為：V=A-1。逆矩陣V存在的條件是A的行列式不等于0。格式：格式：V=inv(A)功能：返回方陣A的逆矩陣V。格式：格式：X=det(A)功能：計算方陣A的行列式值。 pipipxx偽逆矩陣函數(shù)偽逆矩陣函數(shù)pinv偽逆矩陣

3、的偽逆矩陣的MATLAB定義：定義：從數(shù)學意義上講，當矩陣A為非方陣時，其矩陣的逆是不存在的。在MATLAB中，為了求線性方程組的需要，把inv(A*A)*A的運算定義為偽逆函數(shù)pinv，這樣對非方陣，利用偽逆函數(shù)pinv可以求得矩陣的偽逆，偽逆在一定程度上代表著矩陣的逆。格式：格式：C=pinv(A)功能：計算非方陣A的偽逆矩陣。3 3線性代數(shù)方程求解線性代數(shù)方程求解 寫成矩陣形式可表示為：AXB 或 XAB。其中系數(shù)矩陣A的階數(shù)為mn。在MATLAB中，引入矩陣除法求解。 (1)求解方程求解方程AX=B 格式：格式：X=AB 條件：矩陣A與矩陣B的行數(shù)必須相等。 (2)求解方程求解方程XA

4、=B 格式：格式： X=B/A 條件：矩陣A與矩陣B的列數(shù)必須相等。pipipxxmnmnmmnnnnbxaxaxabxaxaxabxaxaxa12211212222121111212111一般線性方程組的4矩陣的分解 (1)三角三角(LU)分解函數(shù)分解函數(shù)lu 所謂三角解就是將一個方陣表示成兩個基本三角陣的乘積(A=LU），其中一個為下三角矩陣L，另一個為上三角形矩陣U，因而矩陣的三角分解又叫LU分解或叫LR分解。矩陣 分解的兩個矩陣分別可表示為： pipipxxnnijaA100100012121nnlllLnnnnuuuuuuU00022211211格式一：格式一：L,U=lu(A) 功

5、能：返回一個上三角矩陣U和一個置換下三角矩陣L(即下三角矩陣與置換矩陣的乘積)，滿足A=L*U。格式二：格式二：L,U,P=lu(A)功能：返回上三角矩陣U，真正下三角矩陣L，及一個置換矩陣P(用來表示排列規(guī)則的矩陣)，滿足L*U=P*A；如果P為單位矩陣，滿足A=L*U。(2)正交正交(QR)分解函數(shù)分解函數(shù) 將矩陣A分解為一個正交矩陣與另一個矩陣的乘積稱為矩陣A的正交分解。格式一：格式一：Q, R=qr(A)功能：產(chǎn)生與A同維的上三角矩陣R和一個實正交矩陣或復歸一化矩陣Q，滿足：A=Q*R，Q*Q=I。格式二：格式二：Q,R,E=qr(A)功能：產(chǎn)生一個置換矩陣E，一個上三角矩陣R(其對角

6、線元素降序排列)和一個歸一化矩陣Q，滿足A*E=Q*R；5 5奇異值分解奇異值分解 矩陣A的奇異值和相應的一對奇異矢量u、v滿足：同樣利用奇異值構(gòu)成對角陣，相應的奇異矢量作為列構(gòu)成兩個正交矩陣U、V，則有：其中AT表示轉(zhuǎn)置矩陣。由于U和V正交，因此可得奇異值分解：格式一：格式一：U,S,V=svd(x)功能：返回3個矩陣，使得X=U*S*V。其中S為與X相同維數(shù)的矩陣，且其對角元素為非負遞減。格式二：格式二：S=svd(A)功能：返回奇異值組成的向量。pipipxxuAvvuATUAV TTVUATVUA6 6矩陣的特征值分析矩陣的特征值分析 矩陣A的特征值 和特征矢量 ，滿足：以特征值構(gòu)成對

7、角陣 ，相應的特征矢量作為列構(gòu)成矩陣V，則有：如果V為非奇異，則上式就變成了特征值分解：格式一：格式一：d=eig(A)功能：返回方陣A的全部特征值所構(gòu)成的向量。格式二：格式二：V,D=eig(A)功能：返回矩陣V和D。其中對角陣D的對角元素為A的特征值，V的列向量是相應的特征向量，使得A*V=V*D。pipipxxvvAvVAV1VVA7 7矩陣的冪次運算矩陣的冪次運算: : Ap 在MATLAB中，矩陣的冪次運算是指以下兩種情況： 1、矩陣為底數(shù)，指數(shù)是標量的運算操作； 2、底數(shù)是標量，矩陣為指數(shù)的運算操作。 兩種情況都要求矩陣是方陣，否則，將顯示出錯信息。(1) 矩陣的正整數(shù)冪矩陣的正整

8、數(shù)冪 如果A是一個方陣，p是一個正整數(shù)，那么冪次表示A自己乘p次。(2) 矩陣的負數(shù)冪矩陣的負數(shù)冪 如果A是一個非奇異方陣，p是一個正整數(shù)，那么A(p)表示inv(A)自己乘p次。 (3) 矩陣的分數(shù)冪矩陣的分數(shù)冪如果A是一個方陣，p取分數(shù)，它的結(jié)果取決于矩陣的特征值的分布。(4) 矩陣的元素冪、按矩陣元素的冪矩陣的元素冪、按矩陣元素的冪 利用運算符“A.p”實現(xiàn)矩陣的元素冪或按矩陣元素的冪運算。 pipipxx8 8矩陣結(jié)構(gòu)形式的提取與變換矩陣結(jié)構(gòu)形式的提取與變換 (1) 矩陣左右翻轉(zhuǎn)函數(shù)矩陣左右翻轉(zhuǎn)函數(shù)fliplr( )格式：格式：X=fliplr(A)(2) 矩陣上下翻轉(zhuǎn)函數(shù)矩陣上下翻轉(zhuǎn)

9、函數(shù)flipud格式：格式：X=flipud(A)(3) 矩陣階數(shù)重組函數(shù)矩陣階數(shù)重組函數(shù)reshape格式一：格式一：X=reshape(A,n,m)功能：將矩陣A中的所有元素按列的秩序重組成nm階矩陣X，當A中沒有mn個元素時會顯示出錯信息。格式二格式二 ：X=reshape(A,m,n,p,.) 或或 X=reshape(A,m,n,p,.)功能：從A中形成多維陣列(mnp.)。pipipxx(4) 矩陣整體反時針旋轉(zhuǎn)函數(shù)矩陣整體反時針旋轉(zhuǎn)函數(shù)rot90( )格式一：格式一：X=rot90(A)功能：將矩陣按反時針旋轉(zhuǎn)90o。格式二：格式二：X=rot90(A, k)功能：將矩陣按反時針

10、旋轉(zhuǎn)k*90o，其中k應為整數(shù)。(5) 對角矩陣和矩陣的對角化函數(shù)對角矩陣和矩陣的對角化函數(shù)diag( )格式一：格式一：X=diag(A,k)功能：當A為n元向量時，可得n+abs(k)階的方陣X，其A的元素處于第k條對角線上；k=0表示主對角線，k0表示在主對角線之上，k0表示主對角線之上，k0表示主對角線之上，k=1.0e-6)&(n=1e6)x=x1;a=f(x); b=df(x);x1=h(a,b,x)n=n+1;endx1, n 要迭代的函數(shù)為要迭代的函數(shù)為其導數(shù)為其導數(shù)為編制的編制的M文件如下：文件如下：33ex即該題的解為即該題的解為 2.8214，與例，與例2-2-6的結(jié)果相

11、的結(jié)果相同。同。 取初始值取初始值x0=3, 運行結(jié)果為運行結(jié)果為x1 = 2.8244x1 = 2.8214x1 = 2.8214x1 = 2.8214n = 3 n=input(n=); %example n=5L=linspace(-2,2,320); X,Y=meshgrid(L); Z=X+i*Y; Z(abs(Z)eps)=1; for k=1:90; % 90次迭代次迭代 Z=Z-(Z.n-1)./(n.*Z.(n-1); %牛頓迭代式牛頓迭代式 end rand(state,0) Co=rand(n,3); R=ones(size(Z); G=R; B=G; z=solve(x

12、,num2str(n),-1); z=double(z); for k=1:n; R(abs(Z-z(k).20.03)=Co(k,1); G(abs(Z-z(k).20.03)=Co(k,2); B(abs(Z-z(k).20時，時，tail=1。 輸出參數(shù)輸出參數(shù)h=0，表示接受原假設；若，表示接受原假設；若h=1，則拒絕原假設；，則拒絕原假設；p表示在原假設下樣本均值出現(xiàn)的概率，表示在原假設下樣本均值出現(xiàn)的概率，p越大，原假設越可越大，原假設越可信；信；ci是是0的置信區(qū)間。的置信區(qū)間。【例【例2-3-6】仍以例】仍以例2-3-4中的樣本數(shù)據(jù)為例，設總體方差已中的樣本數(shù)據(jù)為例，設總體方差

13、已知，假設檢驗中的知，假設檢驗中的0分別取分別取50.00和和50.37，分別取分別取0.60和和0.62，利用假設檢驗函數(shù)觀測樣本接受那種假設。利用假設檢驗函數(shù)觀測樣本接受那種假設?！窘狻吭诿畲爸墟I入：【解】在命令窗中鍵入：x=50.6,50.8,49.9,50.3,50.4,51.0, 49.7,51.2,51.4,50.5,49.3,49.6,50.6,50.2,50.9,49.6;h,p,ci=ztest(x,50.00,0.60,0.05,0) 返回：返回：h = 1p = 0.0124ci = 50.0810 50.6690 h=1,說明拒絕接受原假設，且可看出原假設下樣本均值出

14、現(xiàn)說明拒絕接受原假設，且可看出原假設下樣本均值出現(xiàn)的概率的概率p很小。很小。若取若取0=50.37，=0.62，則得到，則得到h =0，p =0.9743，ci =50.0712 ， 50.6788， 說明接受原假設，且原假設下樣本說明接受原假設，且原假設下樣本均值出現(xiàn)的概率均值出現(xiàn)的概率p較大。較大。 兩個樣本均值差異的檢驗兩個樣本均值差異的檢驗以上檢驗是對單一總體而言，實踐中我們還可能遇以上檢驗是對單一總體而言，實踐中我們還可能遇到對兩個總體進行假設檢驗的問題。對于這類問題，到對兩個總體進行假設檢驗的問題。對于這類問題，我們可利用我們可利用MATLAB提供的函數(shù)提供的函數(shù)ttest2(兩樣

15、本的兩樣本的t檢檢驗驗)來檢驗兩樣本均值之間是否存在顯著差異。該函來檢驗兩樣本均值之間是否存在顯著差異。該函數(shù)的調(diào)用格式為：數(shù)的調(diào)用格式為： h,p,ci=ttest2(x,y,alpha,tail)其中，其中，x和和y是兩個樣本數(shù)據(jù)。是兩個樣本數(shù)據(jù)。【例【例2-3-7】設甲、乙車床生產(chǎn)同一種零件?，F(xiàn)從這兩個車床生】設甲、乙車床生產(chǎn)同一種零件?，F(xiàn)從這兩個車床生產(chǎn)的產(chǎn)品中分別抽取樣品，測得其外徑（單位：產(chǎn)的產(chǎn)品中分別抽取樣品，測得其外徑（單位：mm）為：）為：甲：甲：15.0 14.5 15.2 15.5 14.8 15.1 15.2 14.8乙：乙：15.2 15.0 14.8 15.2 15

16、.0 15.0 14.8 15.1 14.8假定外徑都服從正態(tài)分布，試問兩車床加工的精度是否有顯著假定外徑都服從正態(tài)分布，試問兩車床加工的精度是否有顯著差異？差異？【解】首先假定【解】首先假定1=2，tail取取0。若接受該假設，則應返回若接受該假設，則應返回h=0。在命令窗中鍵入：在命令窗中鍵入：x=15.0 14.5 15.2 15.5 14.8 15.1 15.2 14.8;y=15.2 15.0 14.8 15.2 15.0 15.0 14.8 15.1 14.8; h,p,ci=ttest2(x,y,0.01,0) 運行結(jié)果：運行結(jié)果：h = 0p = 0.8435ci = -0.3

17、227 0.3700 可知，接受原假設，即認為兩車床加工的精度無明顯差異?？芍邮茉僭O，即認為兩車床加工的精度無明顯差異。 在物理實驗中，有時不能直接寫出自在物理實驗中，有時不能直接寫出自變量變量x和因變量和因變量y的函數(shù)關系式，而只能得的函數(shù)關系式，而只能得到函數(shù)在若干個點的函數(shù)值。當要求知道到函數(shù)在若干個點的函數(shù)值。當要求知道觀測點之外的值時，就需要根據(jù)觀測點的觀測點之外的值時，就需要根據(jù)觀測點的值來構(gòu)造一個函數(shù)關系值來構(gòu)造一個函數(shù)關系y=f(x)，并使該函數(shù)，并使該函數(shù)在觀測點的值等于或盡量接近已知的數(shù)值。在觀測點的值等于或盡量接近已知的數(shù)值。當然，構(gòu)造這樣的函數(shù)關系的方法有很多。當

18、然，構(gòu)造這樣的函數(shù)關系的方法有很多。這里，我們僅介紹經(jīng)常要用到的曲線擬合這里，我們僅介紹經(jīng)常要用到的曲線擬合和插值這兩種處理觀測數(shù)據(jù)的方法。和插值這兩種處理觀測數(shù)據(jù)的方法。 2.3.3 曲線擬合與插值曲線擬合與插值 此格式求已知數(shù)據(jù)此格式求已知數(shù)據(jù)x, y的的n階擬合多項式階擬合多項式f(x)的的系數(shù)系數(shù)p。其中，。其中，x必須是單調(diào)的。必須是單調(diào)的。(1)曲線擬合曲線擬合 已知離散點上的數(shù)據(jù)集已知離散點上的數(shù)據(jù)集(x1,y1),(x2,y2),(xn,yn)，求，求解析函數(shù)解析函數(shù)y=f(x)，并使，并使f(x)在原離散點在原離散點xi上盡可能接上盡可能接近給定近給定yi的值，這一過程叫作曲

19、線擬合。最常用的曲的值，這一過程叫作曲線擬合。最常用的曲線擬合是最小二乘法曲線擬合，擬合結(jié)果可使誤差線擬合是最小二乘法曲線擬合，擬合結(jié)果可使誤差的平方和最小。也就是說，該方法是尋求使誤差平的平方和最小。也就是說，該方法是尋求使誤差平方和最小的函數(shù)方和最小的函數(shù)f(x)。調(diào)用格式：調(diào)用格式：p=polyfit(x,y,n) 在求得系數(shù)在求得系數(shù)p以后，可以利用命令函數(shù)以后，可以利用命令函數(shù)polyval求求擬合多項式擬合多項式f(x)在某點（例如在某點（例如x1）的值）的值y1。 調(diào)用格式為：調(diào)用格式為： yi=polyval(p,xi) 。 求擬合多項式求擬合多項式 f(x)，畫出擬合曲，畫出

20、擬合曲線并與實驗數(shù)據(jù)點線并與實驗數(shù)據(jù)點比較。比較。例例【 2-3-8】設在實驗中得到下列兩組數(shù)據(jù)：】設在實驗中得到下列兩組數(shù)據(jù)：圖圖2-3-3 擬合曲線效果x：0.5,1.0,1.5,2.0,2.5,3.0y：1.75,2.45,3.81,4.80,7.00,8.60【解】在命令窗中鍵入：【解】在命令窗中鍵入：在命令窗中鍵入：在命令窗中鍵入：x=0.5 1.0 1.5 2.0 2.5 3.0; y=1.75 2.45 3.81 4.80 7.00 8.60;p=polyfit(x,y,2)x1=0.5:0.05:3.0;y1=polyval(p,x1);plot(x,y,*b,x1,y1,-r

21、) 返回：返回：p = 0.5614 0.8287 1.1560 擬合多項式為：擬合多項式為：f(x)0.5614x2+0.8287x+ 1.1560。繪制的圖形如圖繪制的圖形如圖2-3-3所示。所示。 圖圖2-3-3 擬合曲線效果擬合曲線效果 所謂插值，就是在已知離散點上的數(shù)據(jù)集所謂插值，就是在已知離散點上的數(shù)據(jù)集(x1,y1),(x2,y2),(xn,yn)情況下，找出一解析函數(shù)連接情況下，找出一解析函數(shù)連接自變量相鄰的兩個點自變量相鄰的兩個點(xi,xi+1)，并求得兩點間的數(shù)值。，并求得兩點間的數(shù)值。 調(diào)用格式：調(diào)用格式： yi=interp1(x,y,xi,method) 一維插值一

22、維插值 zi=interp2(x,y,z,xi,yi,method) 二維插值二維插值 vi=interp3(x,y,z,v,xi,yi,zi,method) 三維插值三維插值其中，其中，method表示插值方法：表示插值方法：nearest（最鄰近插（最鄰近插值）；值）；linear（線性插值，默認）；（線性插值，默認）；spline（三次樣條（三次樣條函數(shù)插值）；函數(shù)插值）；cubic（三次函數(shù)插值（三次函數(shù)插值 ）。）。 (2) 插值插值【例【例2-3-9】 仍以例仍以例2-3-8為例。在每兩個數(shù)據(jù)點為例。在每兩個數(shù)據(jù)點之間插入一個值。之間插入一個值?！窘狻吭诿畲爸墟I入：【解】在命令窗

23、中鍵入：x=0.5 1.0 1.5 2.0 2.5 3.0; y=1.75 2.45 3.81 4.80 7.00 8.60;x1=linspace(0.7,2.8,5) %插入插入5個點個點y1=interp1(x,y,x1, cubic) %三次函數(shù)插值三次函數(shù)插值plot(x,y,bo,x1,y1,r*) x1 =0.7000 1.2250 1.7500 2.2750 2.8000 y1 =1.9609 3.0265 4.2775 5.9643 8.0275 插入的這插入的這5個點的坐標為：個點的坐標為：圖圖2-3-4 例例2-3-9圖圖【例【例2-3-10】二維插值舉例?！慷S插值舉例

24、。 【解】在命令窗中鍵入：【解】在命令窗中鍵入：clear allx=1:6; y=1:4; %給出自變量數(shù)據(jù)給出自變量數(shù)據(jù)t=12,10,11,11,13,15;16,22,28,35,27,20;18,21,26,32,28,25;20,25, 30,33,32,30;subplot(1,2,1); mesh(x,y,t) %畫出插值前的圖畫出插值前的圖 x1=1:0.1:6; y1=1:0.1:4; %將將x、y細化為細化為51個點個點 x2,y2=meshgrid(x1,y1); %產(chǎn)生網(wǎng)絡數(shù)據(jù)點產(chǎn)生網(wǎng)絡數(shù)據(jù)點 t1=interp2(x,y,t,x2,y2,cubic); subplo

25、t(1,2,2); mesh(x1,y1,t1) %畫出插值后的圖畫出插值后的圖圖圖2-3-5 例例2-3-10圖圖0246123410152025303502461234010203040運行結(jié)果如圖運行結(jié)果如圖2-3-5所示。所示。 (3) 曲線擬合工具曲線擬合工具 現(xiàn)在現(xiàn)在, 我們來介紹方便實用的曲線擬合工具我們來介紹方便實用的曲線擬合工具和多項式擬合工具。和多項式擬合工具。 曲線擬合工具曲線擬合工具(Curve Fitting Tool)調(diào)用格式：調(diào)用格式： cftool (x,y) x , y是數(shù)據(jù)向量是數(shù)據(jù)向量【例【例2-3-11】利用曲線擬合工具對數(shù)據(jù)點進行擬】利用曲線擬合工具對

26、數(shù)據(jù)點進行擬合和處理。合和處理?！窘狻吭诿畲爸墟I入：【解】在命令窗中鍵入： 在命令窗中鍵入（也可編寫在命令窗中鍵入（也可編寫M文件）文件） x=0.5,1.0, 1.5, 2.0, 2.5, 3.0; y=1.75, 2.45, 3.81, 4.80, 7.00, 8.60; cftool(x,y) 運行后，彈出如圖運行后，彈出如圖2-3-6 所示的所示的cftool界面。單擊界面。單擊“Data”按鈕，彈出按鈕，彈出“Data”對話框，可用來對話框，可用來指定要分析的數(shù)據(jù)以及指定要分析的數(shù)據(jù)以及選擇平滑處理；單擊選擇平滑處理；單擊“Fitting”按鈕，打開按鈕，打開“Fitting”對話

27、框。點對話框。點“New fit”，填寫相應選，填寫相應選項后點擊項后點擊“Apply”按鈕，按鈕，進行曲線擬合，給出進行曲線擬合，給出 圖圖2-3-7 所示的圖形。所示的圖形。圖圖2-3-6 曲線擬合工具界面還有曲線擬合工具界面還有“Plotting”和和“Analysis”等按鈕，可以實現(xiàn)數(shù)據(jù)處理的其它功能，如數(shù)據(jù)等按鈕，可以實現(xiàn)數(shù)據(jù)處理的其它功能，如數(shù)據(jù)的統(tǒng)計分析以及給出殘差圖等，讀者可自行練習。的統(tǒng)計分析以及給出殘差圖等，讀者可自行練習。 多項式擬合工具多項式擬合工具調(diào)用格式：調(diào)用格式：polytool(x,y) 對列向量對列向量x和和y的數(shù)據(jù)進行擬合，并顯示結(jié)的數(shù)據(jù)進行擬合，并顯示結(jié)

28、果的交互圖。該圖為圖形用戶界面，可以觀察果的交互圖。該圖為圖形用戶界面，可以觀察擬合多項式最高階次改變帶來的影響。擬合多項式最高階次改變帶來的影響。polytool(x,y,n) 擬合一個最高階次為擬合一個最高階次為n的多項式的多項式polytool(x,y,n,alpha) 繪制預測值的繪制預測值的100(1-alpha) ,即置信區(qū)間。即置信區(qū)間。 圖圖2-3-7 擬合曲線擬合曲線圖圖2-3-6 cftool界面界面【例【例2-3-11】利用多項式擬合工具對數(shù)據(jù)點進】利用多項式擬合工具對數(shù)據(jù)點進行擬合和處理。行擬合和處理?！窘狻俊窘狻吭诿畲爸墟I入或編寫下列在命令窗中鍵入或編寫下列M文件：

29、文件：y=85,162,230,289,339,381,413,437,452,458,456,440,400,356;x=1,2,3,4,5,6,7,8,9,10,11,12,13,14; polytool(x,y) 圖圖2-3-8 線性擬合線性擬合圖圖2-3-9 三次多項式擬合三次多項式擬合 運行后彈出如圖運行后彈出如圖2-3-8所示的工具圖形界面。所示的工具圖形界面。在圖在圖2-3-8所示的所示的X軸編輯框中輸入數(shù)值或在圖中軸編輯框中輸入數(shù)值或在圖中拖拉垂向參考線，可以評價函數(shù)。拖拉垂向參考線，可以評價函數(shù)。 當鼠標越過垂線時其形狀由箭頭變?yōu)槭之斒髽嗽竭^垂線時其形狀由箭頭變?yōu)槭中?，?/p>

30、示該直線是可以拖拉的。當拖拉參考線時，形，表示該直線是可以拖拉的。當拖拉參考線時，y的預測值將會更新。的預測值將會更新。 變量變量n控制擬合多項式的階次，在圖形頂部控制擬合多項式的階次，在圖形頂部的彈出菜單中進行選擇可以改變多項式的階次。的彈出菜單中進行選擇可以改變多項式的階次。 完成以上操作后，單擊完成以上操作后，單擊“close”按鈕。按鈕。 圖圖2-3-8 所示的是默認時的輸出線性擬合圖。在所示的是默認時的輸出線性擬合圖。在“Export”下拉列表框中進行選擇，可以輸出相應下拉列表框中進行選擇，可以輸出相應的統(tǒng)計量。例如，在的統(tǒng)計量。例如，在“Degree”窗口中輸入窗口中輸入3，進行，

31、進行三次多項式擬合，得圖三次多項式擬合，得圖2-3-9。圖圖2-3-9 三次多項式擬合三次多項式擬合 圖圖2-3-8 線性擬合線性擬合 MATLAB提供的偏微分方程數(shù)值解工具箱提供的偏微分方程數(shù)值解工具箱(PDE Toolbox) 可以幫助我們較為方便地解決這可以幫助我們較為方便地解決這些復雜問題，它可以用指令進行操作，也可以在些復雜問題，它可以用指令進行操作，也可以在圖形用戶界面圖形用戶界面(GUI)中進行操作。本節(jié)主要介紹中進行操作。本節(jié)主要介紹圖形用戶界面中的操作。圖形用戶界面中的操作。2-4 偏微分方程的偏微分方程的GUI解法解法 物理學中有許多問題都可歸結(jié)為偏微分方程物理學中有許多問

32、題都可歸結(jié)為偏微分方程的求解問題。例如，彈性體的振動、電磁波的傳的求解問題。例如，彈性體的振動、電磁波的傳播、熱傳導、粒子擴散、流體運動和場的分布等。播、熱傳導、粒子擴散、流體運動和場的分布等。而要得到這些復雜問題的精確解則是非常困難的，而要得到這些復雜問題的精確解則是非常困難的，我們只能借助于數(shù)值解法。我們只能借助于數(shù)值解法。rhu c uqugnn Dirichlet條件條件 Neumann條件條件 其中，其中， 是是上的單位外法向矢量，上的單位外法向矢量，h、r、q和和g是定義在是定義在程的系數(shù)則可以依賴于時間程的系數(shù)則可以依賴于時間t。2.4.1 偏微分方程的基本類型及邊界條件偏微分方

33、程的基本類型及邊界條件(1) 方程的基本類型方程的基本類型在在PDE Toolbox工具箱求解的基本方程有橢圓型方程、工具箱求解的基本方程有橢圓型方程、拋物型方程、雙曲型方程、特征值方程、橢圓型方程組拋物型方程、雙曲型方程、特征值方程、橢圓型方程組以及非線性橢圓型方程，其形式見表以及非線性橢圓型方程，其形式見表2.4.1。(2) 邊界條件邊界條件一般情況下方程的邊界條件類型有下面幾種：一般情況下方程的邊界條件類型有下面幾種：上的函數(shù)，也可依賴于上的函數(shù)，也可依賴于u ，對于拋物型方程和雙曲型方，對于拋物型方程和雙曲型方表表2.4.1 偏微分方程的基本類型偏微分方程的基本類型1212111ruh

34、uh2222121ruhuh11112211 1122121122221 12222cucuq uq ugcucuq uq ugnnnn1212111ruhuh11112211 112211121122221 1222212cucuq uq ughcucuq uq ughnnnn Dirichlet條件為條件為 , Neumann條件為條件為 的計算要使得的計算要使得Dirichlet條件滿足。在有限元方法條件滿足。在有限元方法對方程組情況：對方程組情況： 混合邊界條件為混合邊界條件為 其中，其中，中，中，Dirichlet條件也稱為本質(zhì)邊界條件，條件也稱為本質(zhì)邊界條件，Neumann條件也稱

35、條件也稱為自然邊界條件。為自然邊界條件。 在進行有限元計算以前，要先對形體進行離散在進行有限元計算以前，要先對形體進行離散化處理，用簡單形體來逼近和代替實際的形體，即化處理，用簡單形體來逼近和代替實際的形體，即建立幾何模型。然后，根據(jù)要求解的問題賦予初始建立幾何模型。然后，根據(jù)要求解的問題賦予初始條件、邊界條件和方程參數(shù)等。組合各個方面的貢條件、邊界條件和方程參數(shù)等。組合各個方面的貢獻后，組成一個大型的矩陣方程，并采用一定的方獻后，組成一個大型的矩陣方程，并采用一定的方法求得方程的解。法求得方程的解。(3) 解偏微分方程的基本思路解偏微分方程的基本思路 MATLAB是采用有限元法來求解偏微分方

36、程是采用有限元法來求解偏微分方程的。該方法利用剖分插值把區(qū)域內(nèi)需求解的微分的。該方法利用剖分插值把區(qū)域內(nèi)需求解的微分方程離散成求解線性代數(shù)方程組，以近似解代替方程離散成求解線性代數(shù)方程組，以近似解代替精確解。精確解。 鑒于求解偏微分方程的復雜性，我們主要介紹鑒于求解偏微分方程的復雜性，我們主要介紹圖形用戶界面（圖形用戶界面（GUI）操作方法。）操作方法。2.4.2 利用利用GUI求解偏微分方程的一般步驟求解偏微分方程的一般步驟 在命令窗中鍵入在命令窗中鍵入pdetool后回車，返回下面的用戶界面：后回車，返回下面的用戶界面： 在在Options菜單中用鼠標指向菜單中用鼠標指向Applicati

37、on選項，選項，會彈出一個子菜單供選擇；也可直接在工具條中會彈出一個子菜單供選擇；也可直接在工具條中單擊單擊Generic Scalar下拉列表框，在其中進行選擇。下拉列表框，在其中進行選擇。表表2.4.2列出了可供選擇的模式。列出了可供選擇的模式。 一般而言，利用一般而言，利用PDE圖形用戶界面求解圖形用戶界面求解PDE問題的過程可分為幾個步驟：問題的過程可分為幾個步驟：（1）選擇應用模式）選擇應用模式 PDE工具箱根據(jù)實際問題的不同提供了很多工具箱根據(jù)實際問題的不同提供了很多應用模式，用戶可以選擇適當?shù)哪Ｐ瓦M行建模和應用模式，用戶可以選擇適當?shù)哪Ｐ瓦M行建模和分析。分析。表表2.4.2 PD

38、E工具箱中可供選擇的模式工具箱中可供選擇的模式序號序號模模 式式序號序號模模 式式1Generic Scalar 一般標量一般標量模式（默認）模式（默認）6Magnetostatics 靜磁學應用模式靜磁學應用模式2Generic System 一般系統(tǒng)一般系統(tǒng)模式模式7AC Power Electromagnetics交交流電電磁學流電電磁學3Structural Mech.,Plane Stress 結(jié)構(gòu)力學平面應力結(jié)構(gòu)力學平面應力8Conductive Media DC 導電介質(zhì)導電介質(zhì)直流電直流電4Structural Mech.,Plane Strain 結(jié)構(gòu)力學平面應變結(jié)構(gòu)力學平面

39、應變9Heat Tranfer 熱傳導應用模式熱傳導應用模式5Electrostatics 靜電學應用靜電學應用模式模式10Diffusion 擴散問題應用模式擴散問題應用模式建立幾何模型建立幾何模型 利用利用Draw菜單中的選菜單中的選項或工具條中的前五個項或工具條中的前五個工具按鈕可畫出工具按鈕可畫出PDE問問題的幾何模型。這五個題的幾何模型。這五個工具按鈕分別為：繪制工具按鈕分別為：繪制矩形或方形矩形或方形 ；繪制；繪制同心矩形或方形同心矩形或方形 ；繪橢圓或圓繪橢圓或圓 ；繪同；繪同心橢圓或圓心橢圓或圓 ；繪；繪多義線多義線 。圖圖2-4-1 建立幾何模型建立幾何模型 在在Option

40、s菜單中選擇菜單中選擇Grid選項可在圖中顯示網(wǎng)格。選項可在圖中顯示網(wǎng)格。單擊單擊Grid Spacing選項，可打選項，可打開開Grid Spacing 對話框調(diào)整對話框調(diào)整網(wǎng)絡間隔的大??；選擇網(wǎng)絡間隔的大??；選擇Snap選項，則鼠標在圖中點擊時選項，則鼠標在圖中點擊時將自動選擇離該點最近的網(wǎng)將自動選擇離該點最近的網(wǎng)絡交點。另外，在絡交點。另外，在“set”框中框中要消去原來的公式表達，根要消去原來的公式表達，根據(jù)問題的要求鍵入新的公式據(jù)問題的要求鍵入新的公式表達。作為例子，圖表達。作為例子，圖2-4-1畫畫出了兩個同心圓的幾何模型。出了兩個同心圓的幾何模型。 圖圖2-4-1 建立幾何模型建

41、立幾何模型(3)定義邊界條件定義邊界條件 在在Boundary菜單中選擇菜單中選擇Boundary mode選項或直接選項或直接在工具條中點擊按鈕，則顯示幾何模型的外邊界和內(nèi)在工具條中點擊按鈕，則顯示幾何模型的外邊界和內(nèi)邊界。在圖中單擊邊界，則邊界變成黑色，表示被選邊界。在圖中單擊邊界，則邊界變成黑色，表示被選中（按住中（按住Shift鍵，可以連續(xù)選擇）。然后，雙擊邊界鍵，可以連續(xù)選擇）。然后，雙擊邊界打開如圖打開如圖2-4-2所示的所示的 “Boundary Conditions”對話框。對話框。圖圖2-4-2 “Boundary Conditions”對話框?qū)υ捒?從圖中可知，最上面一行顯

42、示邊界條件方程，左從圖中可知，最上面一行顯示邊界條件方程，左邊給出條件類型的邊給出條件類型的Neumann和和Dirichlet兩種邊界條件兩種邊界條件的選擇按鈕；右邊的系數(shù)值欄中可輸入邊界條件公式的選擇按鈕；右邊的系數(shù)值欄中可輸入邊界條件公式中的系數(shù)值。中的系數(shù)值。 在在PDE菜單中選擇菜單中選擇PDE mode選項；圖形選項；圖形將顯示為將顯示為PDE模式；在該菜單中單擊模式；在該菜單中單擊PDE Specification選項，可以打開選項，可以打開“PDE Specification”對話框來選擇對話框來選擇PDE問題的類型；問題的類型；右邊的右邊的“Coefficient”欄用來設置

43、欄用來設置PDE方程的系方程的系數(shù)值。數(shù)值。（4）定義）定義PDE類型和類型和PDE系數(shù)系數(shù) 在工具條中單擊按在工具條中單擊按鈕鈕 ，或在，或在Mesh菜單中選菜單中選擇擇Initialize mesh選項，可選項，可以實現(xiàn)研究域的三角形網(wǎng)以實現(xiàn)研究域的三角形網(wǎng)格初始化，如圖格初始化，如圖2-4-3所示。所示。在工具條中單擊在工具條中單擊 ，可對可對初始網(wǎng)格進行細化。利用初始網(wǎng)格進行細化。利用細化后的網(wǎng)絡進行計算，細化后的網(wǎng)絡進行計算，可以獲得具有更高精度的可以獲得具有更高精度的解。解。（5）三角形網(wǎng)絡剖分）三角形網(wǎng)絡剖分圖圖2-4-3 初始化網(wǎng)絡初始化網(wǎng)絡 (6) 在在Mesh菜單中選擇菜單

44、中選擇Parameters選項，選項，打開打開“Mesh Parameters”對話框，可在其中進對話框，可在其中進行設置，明確與網(wǎng)絡剖分有關的參數(shù)。行設置，明確與網(wǎng)絡剖分有關的參數(shù)。 按鈕或在按鈕或在Solve菜單中選擇菜單中選擇Solve PDE選項，選項，可對前面定義的可對前面定義的PDE問題進行求解。問題進行求解。 在在Solve菜單中單擊菜單中單擊Parameters選項，打開選項，打開“Solve Parameters”對話框并進行設置，可以對話框并進行設置，可以確定求解方法和參數(shù)。解的圖形表達確定求解方法和參數(shù)。解的圖形表達PDE工具工具箱提供了解的多種圖形表達方式。箱提供了解的

45、多種圖形表達方式。 單擊按鈕單擊按鈕 或在或在Solve菜單中單擊菜單中單擊Parameters選項，可打開選項，可打開“Plot Selection”對對話框。該對話框包括下列選話框。該對話框包括下列選項：項： Plot type 控制列：該列控件控制圖形類型控制列：該列控件控制圖形類型 的選擇，包括色彩圖、等值線圖、矢量圖、的選擇，包括色彩圖、等值線圖、矢量圖、 變形網(wǎng)絡圖、三維圖和演示圖。變形網(wǎng)絡圖、三維圖和演示圖。 Property控制列：定義對接的那一部分進控制列：定義對接的那一部分進行行 圖形顯示。圖形顯示。 Plot in x-y grid 核選列：在核選列：在x-y網(wǎng)絡中繪圖。

46、網(wǎng)絡中繪圖。 show mesh 核選列：在當前圖顯示網(wǎng)絡。核選列：在當前圖顯示網(wǎng)絡。 Plot solution automatically核選框：系統(tǒng)核選框：系統(tǒng)自自 動繪制解的圖形。動繪制解的圖形。 下面，我們通過一個具體例子來說明如何下面，我們通過一個具體例子來說明如何利用利用GUI求解偏微分方程。求解偏微分方程?！纠纠?-4-1】 帶電云與大地之間存在一個勻強帶電云與大地之間存在一個勻強電場電場E，平行于大地的電纜相當于一根無限長導，平行于大地的電纜相當于一根無限長導體，問在平行于電場方向垂直于電纜的截面上的體，問在平行于電場方向垂直于電纜的截面上的電場如何？電場如何？【解】這是一

47、個靜電學問題。假設整個【解】這是一個靜電學問題。假設整個求解區(qū)域是一正方形，正方形邊界上的電勢為求解區(qū)域是一正方形，正方形邊界上的電勢為y, 再設導體上的電勢為再設導體上的電勢為0。根據(jù)電動力學知識可知，。根據(jù)電動力學知識可知，電勢滿足拉普拉斯方程：電勢滿足拉普拉斯方程： 200| 1,| 1=0 =0.2xyyr正方形界上 使用使用GUI求解這個問題：求解這個問題：(1) 逐層依次點擊逐層依次點擊ptionsApplicationElectrostatic， 即選擇即選擇Electrostatic模式。模式。(2) 區(qū)域設置：單擊菜單欄區(qū)域設置：單擊菜單欄Draw中的中的 Rectangle

48、/Square，在窗口拖拉出一個矩形；，在窗口拖拉出一個矩形； 雙擊雙擊 矩形區(qū)域，在矩形區(qū)域，在Object Dialog對話框中輸入對話框中輸入 Left為為-1，Bottom為為-1，Width為為2，Height為為2， 單擊單擊OK，顯示矩形區(qū)域，顯示矩形區(qū)域R1。點擊。點擊Draw中的中的 Ellipse/Circle，用同樣的辦法作圓，用同樣的辦法作圓C1。輸入。輸入 X-center為為0，Y-center為為0，Radius為為0.2, 單擊單擊 OK，顯示圓形區(qū)域，顯示圓形區(qū)域C1。然后，在。然后，在Set formulals 欄中銷掉原來默認的公式欄中銷掉原來默認的公式,

49、另外鍵入另外鍵入R1-C1。(3) 設置邊界條件：設置邊界條件：選擇選擇Boundary菜單項，單擊邊界模式菜單項，單擊邊界模式Boundary Mode，使邊界變紅，然后分別雙擊每段邊界，打，使邊界變紅，然后分別雙擊每段邊界，打開開Boundary Condition對話框設置邊界條件，矩形對話框設置邊界條件，矩形邊界條件設置為邊界條件設置為 Dirichlet條件，輸入條件，輸入h=1，r=y；圓的邊界在矩形內(nèi)部，自動由銜接條件決定，不圓的邊界在矩形內(nèi)部，自動由銜接條件決定，不必再設邊界條件。必再設邊界條件。(4) 設置方程類型：單擊菜單欄設置方程類型：單擊菜單欄PDE中的中的PDE Sp

50、ecification，設置方程類型為，設置方程類型為Elliptic，輸入，輸入 epsilon為為8.85419e-12，rho為為0，單擊，單擊OK。(5) 網(wǎng)格剖分：網(wǎng)格剖分： 單擊單擊Mesh中的中的Initialize Mesh建立和顯示初建立和顯示初始化三角形網(wǎng)格，再單擊始化三角形網(wǎng)格，再單擊Refine Mesh加密當前加密當前三角形網(wǎng)格。然后單擊三角形網(wǎng)格。然后單擊Jiggle Mesh優(yōu)化網(wǎng)格。優(yōu)化網(wǎng)格。(6) 繪圖：繪圖： 單擊單擊Plot中的中的Parameters，選擇，選擇Contour（等值線），（等值線），Arrows（-grad(u)箭頭），箭頭），Height（3-D plot）（三維圖形），并在）（三維圖形），并在Conto