數(shù)值計(jì)算與MATLAB教材

數(shù)值計(jì)算與MATLAB0前言實(shí)際問題→數(shù)學(xué)模型→數(shù)值計(jì)算方法→程序設(shè)計(jì)→上機(jī)計(jì)算→求出結(jié)果

“數(shù)值計(jì)算”、“數(shù)值分析”為同義詞，這門學(xué)科有如下特點(diǎn)：面向計(jì)算機(jī)：可實(shí)現(xiàn)性有可靠的理論基礎(chǔ)：理論明確要有好的計(jì)算復(fù)雜性：可行性可實(shí)現(xiàn)數(shù)值驗(yàn)證：可驗(yàn)證性算法誤差明確、并可滿足要求：可信性

主要內(nèi)容：MATLAB基礎(chǔ)（方程組、最小二乘法、回歸分析）、MATLAB數(shù)據(jù)可視化、常用函數(shù)、復(fù)數(shù)、向量計(jì)算、插值法、曲線擬和、數(shù)值微積分、多項(xiàng)式、級(jí)數(shù)展開、常微分方程求解、符號(hào)計(jì)算等，其中有些問題包含在各章節(jié)中，主要以MATLAB為線索，進(jìn)行展開。1MATLAB基礎(chǔ)

1.1基本數(shù)據(jù)類型與運(yùn)算MATLAB數(shù)據(jù)類型及常數(shù)表達(dá)名稱：字母開頭,大小寫敏感，缺省變量：ans

實(shí)數(shù)：科學(xué)計(jì)數(shù)法，特殊數(shù)表達(dá)eps:浮點(diǎn)數(shù)的最小單位；Inf無窮大；NaN不是一個(gè)數(shù)值；pi;realmax;realmin復(fù)數(shù)：a+bi,或a+bj

，a,b為實(shí)數(shù),i、j虛數(shù)單位矩陣：[………],其中各項(xiàng)之間用空格、逗號(hào)、分號(hào)分開，分號(hào)換行1.1基本數(shù)據(jù)類型與運(yùn)算特殊矩陣：eyeonesrandrandnzeros…(n,m)/(n):nbym,nbyn線性向量：線性a:k:b

產(chǎn)生向量示例：1）用360個(gè)弧度表示一圓周；

2）產(chǎn)生200*200個(gè)元素為（0，1）正太分布隨機(jī)數(shù)矩陣；1.2運(yùn)算參與運(yùn)算的基本單位：矩陣單個(gè)實(shí)數(shù)、復(fù)數(shù)、向量都是矩陣，都可以參加運(yùn)算。矩陣運(yùn)算：+-*/\^數(shù)組運(yùn)算：

.*./.^其中矩陣的除法

/:B/A——B*A-1\:A\B——A-1B

如果XA=B則:X=B/A

如果AX=B則:X=A\B‘轉(zhuǎn)置.’非共軛轉(zhuǎn)置矩陣與常數(shù)運(yùn)算示例生成200*200個(gè)元素為復(fù)數(shù)元素矩陣，復(fù)數(shù)實(shí)部與虛部均為（10，5）正態(tài)分布隨機(jī)數(shù)產(chǎn)生2000個(gè)角度為一周均勻分布，直徑為正態(tài)分布的點(diǎn)云，方差為10，圓心為（0，0）點(diǎn)。exp求解1.3矩陣基本操作矩陣的轉(zhuǎn)置矩陣元素訪問及矩陣重構(gòu)A（向量，向量）幾個(gè)常用函數(shù)

fliplr

flipudrot90repmat(A,m,n)

reshape(A,m,n)示例：

A:100*100,B:40*40,將B嵌入到A中心位置

A:300*400,B:300*400,合成“塊交錯(cuò)”的C，塊的大小為5*5，奇數(shù)塊來至A,偶數(shù)塊來至于B1.4線性方程組、最小二乘法與回歸分析線性方程一般形式AX=BX=A\BA、B可以是復(fù)數(shù)？復(fù)數(shù)方程，可以做！A列數(shù)與X行數(shù)，B行數(shù)一致，未知數(shù)與方程數(shù)A行數(shù)少于列數(shù)A行數(shù)多于列數(shù)線性方程一般形式AX=B1）B列數(shù)大于1，相當(dāng)于一個(gè)A對(duì)應(yīng)多個(gè)B的方程2）A行數(shù)小于列數(shù)，相當(dāng)于方程數(shù)小于未知數(shù)個(gè)數(shù)3）A行數(shù)大于列數(shù)，相當(dāng)于方程數(shù)多于未知數(shù)個(gè)數(shù)1.4線性方程組、最小二乘法與回歸分析未知數(shù)多，方程數(shù)少：得到一個(gè)基本解未知數(shù)少，方程數(shù)多：采用最小二乘法，得到一個(gè)最綜合的解，也就是使得“=”左端計(jì)算的值，與B的差最小。最小二乘法：對(duì)于y=ax+b

形式的線性方程，如果已知點(diǎn)(x1,y1)和(x2,y2)滿足方程，則同時(shí)滿足y1=a*x1+b;y2=a*x2+b，就可以建立方程即可求出(a,b)值,如果多于2個(gè)點(diǎn)，則應(yīng)該找到一條直線與點(diǎn)最近，實(shí)際上是求出線到點(diǎn)的y距離的總和最小，數(shù)學(xué)上利用y距離的平方和確定總誤差。即E=SUM(y*-y)2,E中僅包含兩個(gè)未知數(shù)a,b.求E最小可以利用偏導(dǎo)等于0的方程。1.4線性方程組、最小二乘法與回歸分析X=1:5’X=[xoese(size(x))]Y=[1030324548]’A=x\y，即為線性方程的兩個(gè)系數(shù)誤差多大？怎么求？（相對(duì)誤差、絕對(duì)誤差）

如何做非線性回歸x=1:10;y=[20212424.62630.129.7384353]方程形式為y=ax2+bx+c實(shí)際上1.4線性方程組、最小二乘法與回歸分析MATLAB實(shí)現(xiàn):x,y

轉(zhuǎn)換為列向量X=[x.^2xones(size(x))]A=x\y,a(1)->aa(2)->ba(3)->c如何理解ones(size(x))?Size返回什么其他方程形式回歸！1.5環(huán)境管理Clc：清楚控制臺(tái)（操作窗口）Clear：清楚內(nèi)存變量，其他形式？清楚一個(gè)（一批）特定變量help內(nèi)容QuitSAVEfnameXYZ-ASCII-DOUBLeload(filename)whowhos輸出格式設(shè)定formatshort/long/shorte/longe/shortg/longg/formathex/bank/rat2常用函數(shù)與數(shù)據(jù)可視化2.1常用函數(shù)

Magic(n)魔方陣pascal(n)陣hadamard(n)n為2的密正交陣

calendar(yyyy,mm)月歷；clock時(shí)間向量

log2log10logsqrtexpabsimagrealangleconjfix向0取整floorceilround舍入

sign數(shù)單位值

Sincostanasin

acos

atanatan22.2平面曲線繪制plot(y)

如果y為實(shí)數(shù)向量則繪制各點(diǎn)折線，復(fù)數(shù)向量以實(shí)部為x,虛部為y繪制曲線，y為矩陣則以列為線繪制折線plot(X1,Y1,...)plot(X1,Y1,LineSpec,...)繪制（x1,y1）線并指定線的特性線型：-實(shí)線--劃線:點(diǎn)線-.點(diǎn)劃線數(shù)據(jù)標(biāo)志：ox+sd*.^v><ph

線的顏色：ym(agenta)c(yan)rgbw(blac)k

2.2平面曲線繪制繪制120度差的三條sin曲線；繪制多邊形；曲柄連桿機(jī)構(gòu)軌跡；繪制動(dòng)刀軌跡平面幾何變換：點(diǎn)陣的旋轉(zhuǎn)、縮放、平移2.2平面曲線繪制Bar(x,y,width)Polar(theta,rho)Pie(x,explode)SemilogxSemilogyLoglog示例：繪制螺旋線；繪制傳遞函數(shù)的幅頻特性曲線，相頻特性曲線（利用半對(duì)數(shù)曲線)2.3視圖管理與控制1）坐標(biāo)軸axis([xmin

xmax

ymin

ymax

zmin

zmax])axisautoaxismanualaxistightaxisfillaxisequalaxissquareaxisoffaxison2）分割/設(shè)定繪圖區(qū)

Subplot(m,n,p)分割成m行n列，當(dāng)前設(shè)為p列3)Holdon/off4）titlexlabel

ylabel

zlabel

指定標(biāo)題平面曲線示例3數(shù)據(jù)分析k=convhull(x,y):包圍所有點(diǎn)(x,y)的凸多邊形的提取，返回值為向量in=inpolygon(X,Y,xv,yv)多邊形內(nèi)點(diǎn)的檢測(cè),返回0,13.1常規(guī)統(tǒng)計(jì)分析常規(guī)分析f=factor(n)返回n整數(shù)的素?cái)?shù)因子primes(n)

返回小于等于n的所有素?cái)?shù)P=perms(v)

所有排列sum,mean,std,median,mode,:

格式：函數(shù)（a,dim）a為矩陣dim=2為列間計(jì)算max,min

格式：函數(shù)（a,b,dim）a,b相同大小[n,x]=hist(y,m):m為組數(shù)，n為個(gè)組計(jì)數(shù)，x為組中值polyarea(X,Y)多邊形面積

3.2有限差分Diff(x):差分Cumsum（x):累加.求導(dǎo)、積分示例：一點(diǎn)，在X軸上運(yùn)動(dòng)，每3秒觀察一次，數(shù)據(jù)為03255648108912849493，求速度3.3FFT變換傅里葉積分變換實(shí)質(zhì)：FFT(x):返回值包含兩部分，幅值與相位，并且是左右對(duì)稱。示例：s=10*sin(t)+7*sin(3*t+0.8)+5*sin(5*t+1.2)+3*sin(8*t+1.3)3.4數(shù)值濾波-IIF濾波infiniteimpulseresponse數(shù)值濾波處理y(n)=b(1)*x(n)+b(2)*x(n-1)+...+b(nb+1)*x(n-nb)-a(2)*y(n-1)-...-a(na+1)*y(n-na)y=filter(b,a,X)其中a的第1項(xiàng)一定為1如：對(duì)一個(gè)序列向前3項(xiàng)平均處理:x=[001111]y=filter([111]/3,1,x);結(jié)果為：000.330.66711*指數(shù)平滑：y(k)=ax(k)+(1-a)y(k-1)y=filter(0.7,[1-0.3],x);*聲音的ECHO處理，設(shè)回聲時(shí)差為0.3s?,回聲強(qiáng)度0.6?，采樣頻率320004多項(xiàng)式與插值4.1多項(xiàng)式表達(dá)采用實(shí)數(shù)行向量表示降冪多項(xiàng)式系數(shù)，如：p=[3425]表示為4.2多項(xiàng)式運(yùn)算與求解

conv(p1,p2)多項(xiàng)式相乘，卷積計(jì)算

roots(p)求解多項(xiàng)式的根

polyval(p,x)多項(xiàng)式求值

polyfit(x,y,n)擬和多項(xiàng)式Poly(a):如果a為矩陣則返回矩陣A的特征多項(xiàng)式p；如果a為向量，返回以向量中的元素為根的多項(xiàng)式。[q,r]=deconv(u,v)：u/v，q商多項(xiàng)式，r余數(shù)多項(xiàng)式4多項(xiàng)式與插值4.3插值

yi=interp1(x,Y,xi,method)一維插值，method可為：

'nearest‘;'linear';'spline';‘cubic’

*插值：利用插入點(diǎn)X附近的數(shù)據(jù)，按一定方法求出一個(gè)表達(dá)式，再利用表達(dá)式計(jì)算該點(diǎn)

spline:樣條插值；cubic:立方插值

zi=interp2(X,Y,Z,XI,YI,metho):二維插值示例：數(shù)據(jù)插值、曲面插值（10個(gè)值）插值示例Sinsin面z=sin(sqrt(x.^2+y.^2)/2)4多項(xiàng)式與插值4.4

數(shù)值網(wǎng)格化[XI,YI,ZI]=griddata(x,y,z,xi,yi)

。x,y,z為實(shí)測(cè)數(shù)值，xi,yi為網(wǎng)格坐標(biāo)，返回網(wǎng)格坐標(biāo)值，如：x=rand(300,1)*4-2;y=rand(300,1)*4-2;z=x.*exp(-x.^2-y.^2);[xi,yi]=meshgrid(-2:0.2:2);zi=griddata(x,y,z,xi,yi);5三維與特殊形式繪圖5.1三維曲線

plot3(X1,Y1,Z1,LineSpec,...)：三維空間曲線

Bar3(x,y,width)

cylinder(r,n)

5.2三維曲面1）產(chǎn)生網(wǎng)格數(shù)據(jù)[X,Y]=meshgrid(x,y)2）繪制網(wǎng)格曲面

mesh(X,Y,Z)mesh(Z)mesh(...,C)meshc(...)meshz(...)3）繪制曲面surf(Z)surf(X,Y,Z)surf(X,Y,Z,C)surfc(...)4）繪制等值線contour(Z,n)r=sqrt(x.^2+y.^2)+0.1;z=sin(r)./rPeaks(N)5.3曲面平滑與顏色

view(az,el)：觀察方向1）平滑shadingflat/faceted/interp

2）利用colormap(map)

設(shè)定曲面顏色，map為64by3的矩陣，為顏色索引，系統(tǒng)提供的map為：AutumnbonecoolcopperflagGrayhotjetlinesprismspringSummerwinter3）Colormap([r,g,b]);6邏輯函數(shù)邏輯關(guān)系：

<><=>===~=運(yùn)算：&|~All(a)TesttodetermineifallelementsarenonzeroAny(a)Testforanynonzeros

k=find(x);[i,j]=find(X)Findindicesandvaluesofnonzeroelements示例1）將a中大于0的元素加和2）將a中小于0的元素加57程序設(shè)計(jì)基礎(chǔ)7.1函數(shù)的定義

Function[mean,stdev]=stat(x)n=length(x);mean=sum(x)/n;

stdev=sqrt(sum((x-mean).^2/n));其中：Function[mean,stdev]=stat(x)為函數(shù)首部，“=”左側(cè)為返回值形式變量，右側(cè)為函數(shù)名，括號(hào)內(nèi)為入口的形式參數(shù)；其中可以用“%”引導(dǎo)注釋內(nèi)容。文件名應(yīng)和函數(shù)名相同。7.2程序結(jié)構(gòu)控制分支結(jié)構(gòu)

ifexpressionstatementselse

statementsend開關(guān)結(jié)構(gòu)switch表達(dá)式（標(biāo)量或字符串）

case值1

語句組ACase值2

語句組B…….Otherwise

語句組Nend掃描循環(huán)

forvariable=expressionstatementsendExpression：向量邏輯循環(huán)whileexpressionstatementsend示例：返回處入口向量升序、及降序、及相鄰項(xiàng)相等的個(gè)數(shù)function[abc]=myf1(x)l=length(x);%%求x中元素?cái)?shù)量

a=0;b=0;c=0;fork=1:(l-1)ifx(k)>x(k+1)a=a+1;elseifx(k)<x(k+1)b=b+1;elsec=c+1;endendend常規(guī)程序設(shè)計(jì)方法太復(fù)雜，如何運(yùn)用MATLAB的功能？設(shè)計(jì)出更簡(jiǎn)潔程序？function[abc]=myf1(x)

dx=diff(x);a=sum(dx>0);b=sum(dx<0);c=sum(dx==0);EndFuncitony=myf1(x);

dx=sign(diff(x))+2y=accumarray(dx',ones(size(dx)));End;十進(jìn)制乘法計(jì)算，可以使用卷積，如何十進(jìn)制進(jìn)位，字符操作？7.3幾個(gè)函數(shù)Accumarray(s,v):s列向量，v行向量，s,v元素?cái)?shù)量相同根據(jù)s指定，累加v，生成新的向量，其中s元素指出累加的位置，如：v=11:15,s=[21241]’,結(jié)果為：27；24；0；14circshift(A,shiftsize)：循環(huán)移動(dòng)矩陣的行和列，shiftsize

：[l,r]l：移動(dòng)行數(shù)，r：移動(dòng)的列數(shù)，正：右下；負(fù)：左上

a=1:5;a=a’;b=circshift(a,2),結(jié)果為：[45123]’Sort(a,dim)Colormap(v):v是個(gè)n*3向量；hot;gray;cool等，表示colorindexImage(a)：以colorindex形式顯示矩陣Globalv1v2Clf：clearfigurev=input(‘prompt')Ginput(n)Pause(d):控制程序停止ds，d可以是小數(shù)程序設(shè)計(jì)示例分析顯示移動(dòng)的彩色塊移動(dòng)的SIN函數(shù)旋轉(zhuǎn)的多邊形曲柄連桿機(jī)構(gòu)的運(yùn)動(dòng)一種“繪畫”技術(shù)functiona=test1;n=input('Interthelengthof1thline');

clf;holdonaxisequal;axisoffsubp(0,n*exp(i*pi/2));endfunctionb=subp(p,l)

plot([p;p+l]);ifabs(l)>5subp(p+l,0.71*l*exp(i*pi/4));subp(p+l,0.71*l*exp(-i*pi/4));End;end8微分方程數(shù)值解與數(shù)值積分8.1微分方程數(shù)值解

可以計(jì)算出一個(gè)序列的y，即為數(shù)值解，可以使用Cumsum(x)函數(shù)實(shí)現(xiàn)累加過程，但先要給定初值。例如：y’=cos(x),y=?EulerR-K方法微分方程數(shù)值解實(shí)際上是前后兩點(diǎn)斜率的平均值。上式叫二階Runge—Kutta法上式叫四階Runge—Kutta法MATLAB微分方程數(shù)值解法[t，x]=ode23(f，ts，x0，options)[t，x]=ode45(f，ts，x0，options)分別為二階，四階R-K方法，其中f微分方程（組）的描述函數(shù)；ts:為微分點(diǎn)的時(shí)間序列，求出的x,為與ts對(duì)應(yīng)的解序列，x0為初值，options為求解過程誤差控制參數(shù)，可以省略示例11）定義方程函數(shù)

functiondx=mf(t,x);

dx=cos(t);end2）求解

ts=0:0.01:10;[t,x]=ode45(‘mf’,ts,0);Plot(x);%%顯示曲線1）定義方程函數(shù)

functiondy=mf(t,y);

dy=-t*y.^2;end2）求解

ts=0:0.01:1;[t,y]=ode45(‘mf’,ts,2);Plot(y);%%顯示曲線示例2其中：a=20;b=40;c=15;t的范圍[0,0.5],x(0)=0;y(0)=0注意：為二元微分方程1）編輯函數(shù)functiondx=mf(t,x)a=20;b=40;c=15;s=sqrt((c-x(1))^2+(a*t-x(2))^2);

dx=[b*(c-x(1))/s;b*(a*t-x(2))/s];end函數(shù)返回值[dx/dt;dy/dt]2）ts=0:0.01:0.5;x0=[00];

[t,x]=ode45(‘mf’,ts,x0);%%采用缺省的精度求解X(1)——x,x(2)——y示例31)定義方程函數(shù)functiondu=mf(t,u);du=[-10*u(1)+9*u(2);10*u(1)-11*u(2)];end2)求解

ts=0:0.01:0.5;[t,u]=ode45(‘mf’,ts,[12]);示例4定義微分方程函數(shù)functiondy=mf(t,y);

dy=[0;0];dy(1)=y(2);dy(2)=(7*y(2)-5*y(1)+10)/2;end求解：[t,y]=ode45('mf',0:0.01:1,[21])實(shí)際上y(2)’，應(yīng)該與y(1)相同微分simulink解法上個(gè)方程8.2數(shù)值積分trapz(x,y)：梯形法求積分，求面積q=quad(@fun,a,b)：

Simpson法積分示例：定義被積函數(shù)functionyy=mf(x);

yy=1./(x.^3-2*x-5);end積分運(yùn)算Q=quad(@mf,0,2)9函數(shù)極值與優(yōu)化問題fplot(fun,limits):繪制函數(shù)曲線如：fplot(@sin,[0,2*pi]);x=fminbnd(fun,x1,x2):返回單變量函數(shù)fun，在[x1,x2]區(qū)間上的最小值，如：fminbed(@sin,-1,1),結(jié)果為-1x=fminsearch(fun,x0)：返回多變量函數(shù)fun，在x0附近的最小值如：functionf=myf(x);f=100*(x(2)-x(1).^2).^2+(1-x(1)).^2;endx=fminsearch(@myf,[01])9函數(shù)極值與優(yōu)化問題x=bintprog(f,A,b,Aeq,beq)：求多變量函數(shù)線性函數(shù)f，(1/0)規(guī)劃的極小值，約束條件：AX<=b;AeqX=beq;X取值范圍1/0,其中f為向量表示的線性目標(biāo)函數(shù)如：求f(x)=-9x1-5x2-6x3-4x4極小值，滿足

6x1+3x2+5x3+2x4<=9;x3+x4<=1;*x=fgoalattain(fun,x0,goal,weight,A,b,Aeq,beq,lb,ub)：目標(biāo)規(guī)劃x=fmincon(fun,x0,A,b,Aeq,beq,lb,ub)：條件極值9函數(shù)極值與優(yōu)化問題x=fsolve(fun,x0,options)：返回fun在x0的解，如：在[-5-5]附近的解，定義函數(shù)

functionf=myf1(x)f=[2*x(1)-x(2)-exp(-x(1));-x(1)+2*x(2)-exp(-x(2))];end求解：x=fsolve(@myf1,[-5;-5])options=optimset('param1',value1,'param2',value2,...)Display:Levelofdisplay.'off'displaysnooutput;'iter'displaysoutputateachiteration9函數(shù)極值與優(yōu)化問題Fsolve示例分析：四桿機(jī)構(gòu)求解已知量：l1,l2,l3,l4;q1,q4未知量：q2,q3function[bc]=jg();globall1l2l3l4q1q4;l1=20;l2=40;l3=30;l4=50;q4=0;x0=[11];b=[];c=[];op=optimset('Display','off');forq1=0:pi/180:2*pi;x=fsolve(@my,x0,op);x0=x;b=[b;l1*exp(i*q1)];c=[c;l2*exp(i*x(1))];endendfunctionf=my(x)globall1l2l3l4q1q4;%%x(1)-q2;x(2)-q3f=[l1*cos(q1)+l2*cos(x(1))+l3*cos(x(2))-l4*cos(q4);l1*sin(q1)+l2*s