matlab學(xué)習(xí)第四章數(shù)值計算

第四章數(shù)值計 計算指令之間的內(nèi)在聯(lián)系及區(qū)別。對于那些熟悉其他高級語言（FORTRAN，Pascal，C++）的讀者來說，通過本章，對于那些經(jīng)過大學(xué)基本數(shù)學(xué)的讀者來說，通過本章，精良完善的計算指令，自然易讀的程序?qū)⑹顾麄兏形颉啊睌?shù)學(xué)的基礎(chǔ)地位和局限性，看到從“理想化”簡對于那些熟悉基本指令的讀者來說，通過本章，圍繞基本數(shù)值問題展開的由于的基本運算單元是數(shù)組，所以本章內(nèi)容將從矩陣分析、線性代數(shù)的數(shù)值Fourier分析，和一般常微分方程初值、邊值問題。本章的最后稀疏矩陣的處理，因為 6.0LINPACKEISPACK庫為基LAPACKARPACK程序庫的新基礎(chǔ)上。為構(gòu)成空間的基向量時不唯一的，且新版的更。本書新版全部算例結(jié)果是在6.x版 6.0版起，提供了一組專門求微分方程“邊值問題”數(shù)值解的指令。適應(yīng)4.14.52個算例闡述求解細(xì)節(jié)。triplequad。本章新版對此作了相應(yīng)的改變。LU分解、行列式和恰定方程組的ans=ti=erireitd0erdred范數(shù)、條件數(shù)和方程解的精Hxb近似解和準(zhǔn)確解進(jìn)行比較。N=[681012fork=1:length(N)ndx=norm(x_approx-x_exact);nx=norm(x_approx);formatshortedisp('實際誤差disp('近似的最大可能誤差er_approx'),disp(er_approx),disp('')disp('最大可能誤差er_max'),disp(er_max),disp('') 實際誤差1.5410e-0101.7310e-近似的最大可能誤差er_approx3.3198e-0093.3879e-最大可能誤差7.9498e-0073.8709e-特征值和特征向量的--0.0430+0.8333+ 00.8333-D4010R21.0e-014 04.3.1-3eigeigs的比較。T_fullT_partd3.0140+ans3.0140+V_err=D_err=特征值問題的條件formatshortB1000010000110001c_equ=Warning:MatrixisclosetosingularorbadlyResultsmaybeinaccurate.RCOND=1.110223e-InD:\6P1\toolbox\\matfun\condeig.matline30c_eig=ADJ0100000100000100000100000DColumns1through 0-0.0119+ 00-0.0119- 00 0.0323+00 Column00000.0323-c_eig=1.0e+010*c_equ=復(fù)數(shù)特征值對角陣與實數(shù)塊特征值對角陣的轉(zhuǎn)4.3.3-14.3.1-1中的復(fù)數(shù)特征值對角陣DVR*DR/VR=AVR DR 矩陣的譜分解和矩陣4.3.4-1】數(shù)組乘方與矩陣乘方的比較。clear,A=[123;456;789]; Mp0.6962+ 0.4358+ 0.1755-0.63250.7309+0.8292-0.56881.0259-1.4830+-pA_A pA_MsinM 奇異值分解和矩陣結(jié)奇異值分解定矩陣結(jié)構(gòu)的奇異值分解描線性二乘問題矩陣除運算的廣義線性模型的最小二乘【例4.4.2.2-1】對于超定方程yAx，進(jìn)行三種解法比較。其中A取庫中的特A=gallery(5);A(:,1)=[];y=[1.77.56.30.83-Warning:MatrixisclosetosingularorbadlyResultsmaybeinaccurate.RCOND=5.405078e-x=xx=Warning:Rankdeficient,rank=3tol= xxx=0nx=nxxnxxxeeeeee法10----0-015偏導(dǎo)數(shù)和梯理論定數(shù)值計算指4.5.2.2-1】用一個簡單矩陣表現(xiàn)diffgradientF=123456789Dx333333=111111FX111111111FY333333333=222222222=666666666【例4.5.2.2-2】研究偶極子(Dipole)的電勢（Electricpotential）和電場強度（Electricfielddensity）。設(shè)在(ab處有電荷q，在(a,b處有電荷q

V(x,y)

(40

1

，E

r

(xa)2(yb)2,r

(x(xa)2(y

9

title('\fontname{隸書}\fontsize{22}偶極子的場'),holdonxlabel('x');ylabel('y'),hold偶極子的偶極子的6420---- - -0x 6多項式一元函數(shù)的零利 作圖指令獲取初步近似任意一元函數(shù)零點的精確【例4.6.2.2-1】通過求f(t)(sin2t)eatbt的零點，綜合敘述相關(guān)指令的用法 10--10-------0t54.6-令必須在指令窗中運行，并得到如圖4.6-2所示的局部放大圖及鼠標(biāo)操作線。zoomon[tt,yy]=ginput(5);zoomtt t4=y40exitflag=t3=y3exitflag=op=ActiveConstrTol:Display:TolX:2.2204e-TypicalX:ans=t4n=y4nexitflag=多元函數(shù)的零f1(x,y)sin(xy)4.6.3-1】求解二元函數(shù)方程組

f2(x,y)cos(xy)fv=[-0.2,0,holdon,contour(X,Y,F2,v),hold2210----- - - -0 4.6- Optimizationterminated-orderoptimalitylessthanOPTIONS.TolFun,andnonegative/zerocurvaturedetectedxy f1.0e-006 exit1functionff=fun(x)一元函數(shù)的極小值多元函數(shù)的極小值4.7.2-1f(xy100(yx221x)2Rosenbrock's的文件名為exm04072_1_1.m。）sx sfval=sexit1soutputiterations:funcCount:algorithm:'Nelder-MeadsimplexdirectWarning: beprovidedfortrust-regionusingline- ethodInD:\6P1\toolbox\optim\fminunc.matline211Optimizationterminatedsuccessfully:Currentsearchdirectionisadescentdirection,andmagnitudedirectionalderivativeinsearchdirectionlessthan2*options.TolFunux= sfval=uexit1uoutputiterations:funcCount:stepsize:1.2992orderopt:5.0020e-004algorithm:'medium-scale:Quasi-Newtonlinegrid1.0e-003hess820.4028--409.5496一元函數(shù)的數(shù)值積閉型數(shù)值積

1ex2dx，其精確值為0.746842040symsx;IS=int('exp(-x*x)','x',0,1)ISans= Isim=IL=Ig=Ssp=4.8-開型數(shù)值積functiong= x=[0.1488743390;0. w=t=.5*(b+a)+.5*(b-a)*[-flipud(x);x];W=[flipud(w);w];g=sum(W.*feval(fun,t))*(b-4.8.1.2-1f(x)cos(xsymsx;F=int('cos(x)','x',-Fans= aF=1【例4.8.1.2-2】求I ln1dx，準(zhǔn)確結(jié)果是1

88622692 symsWarning:ExplicitintegralcouldnotbeInD:\ 6P5\toolbox\symbolic\@sym\int.matline58In 6P5\toolbox\symbolic\@char\int.matlineIS （2）quad指令求積Warning:Dividebyzero.In \funfun\inlineeval.matlineInD:\ \funfun\@inline\feval.matline34InD:\ \funfun\quad.matline62IsimIg=多重數(shù)值積積分限為常數(shù)的二重積分指

21xydxdy

12xydxdysymsx

x

Warning:ExplicitintegralcouldnotbeInD:\6P5\toolbox\symbolic\@sym\int.matline58ssx01= ssx12 nsx01=nsx12內(nèi)積分限為函數(shù)的二重積function%function%ifisnumeric(x1)==1;xx1=x1*ones(size(y));elsexx1=feval(x1,y);endifisnumeric(x2)==1;xx2=x2*ones(size(y));elsexx2=feval(x2,y);endfori=1:n

4

(x2y2dxdyfunctionf=x_low(y)

1 Warning:Minimumstepsizereached;singularityInD:\ 6P5\toolbox\\funfun\quad.matline88InD:\ 6p5\work\G_yi.matline11InD:\ \funfun\@inline\feval.matline20InD:\ \funfun\quad.matline62InD:\ 6p5\work\double_int.matline8SS=Ssym= 卷“完整”離散序列的數(shù)值卷“截尾”離散序列的數(shù)值卷多項式乘法與離散卷積的算法同連續(xù)時間函數(shù)的數(shù)值卷積的實4.8.3.5-1A(n)

n和

B(n)

n。%axis([nc1-2,nc2+2,0,10]),grid,holdonstem(kc,c,'b'),stem(kx,cc(N+1:end),'g'),holdoff3210 3210 【例4.8.3.5-2】求函數(shù)u(t)etU(t)和h(t)tet/2U(t)的卷積。本例展示：（A）符號Laplace變換求卷積的理論表示；（B） symstao;t=sym('t','positive');HS1=yt1=4.8-dy124.8.3.5-3】用“零階”近似法求u(t)etU(t和h(t)tet2U(t點的確定。（D）離散卷積和連續(xù)卷積之間的關(guān)系。（E）conv的使用。（F）繪圖分% [3,11]段[11,14]段']),disp(dy13)003 003統(tǒng)計分布的數(shù)字特%%%XMoment1arithmetic:[10.80100.6246median:[10.6511geometric:[10.4691harmmonic:[0.09260.16821.8276]Moment2=[31.34290.2951[982.37600.0871[17.83980.2331[99.99000.9317樣本分布的頻數(shù)直方圖描4.9.2-1】hist指令的使用示例。%864520-050-05%98786564432210-010-01概率函數(shù)、分布函數(shù)、逆分布函數(shù)和隨機數(shù)的泊松分布(Poissontext(30,0.07,'\fontsize{12}{\mu}={\lambda}=20')==00正態(tài)分布（Normal%forsubplot(1,3,1),plot(xd,yd,'b');holdontext(mu-0.5*sigma,0.3,num2str(P(1))),holdoffsubplot(1,3,2),plot(xd,yd,'b');holdontext(mu-0.5*sigma,0.3,num2str(P(2))),holdoffsubplot(1,3,3),plot(xd,yd,'b');holdon 005000 52分布（Chi-squareplot(x,yd_c,'b'),holdontext(10,0.16,['\fontsize{16}x~{\chi}^2' text(1.5,0.08,'\fontname{隸書}\fontsize{22}置信水平0.9') holdoff 置信水平005多項式多項式擬合原擬合多項式階數(shù)的確多項式擬合 實for%xlabel('階次'),title('chi2與度')xlabel('階次'),title('Q0.5線')chichi2Q階0.5 8642 246階 階%errorbar(x,YE{3},D{3},'r');holdofftext(0.1,5.5,['chi2='num2str(chi2(3))'~'int2str(freedom(3))])text(0.1,5,['freedom='int2str(freedom(3))])text(0.6,1.7,['Q='num2str(Q(3))6 432101ans ans 非線性最小二乘估偽線性最小二借助fminsearch指令進(jìn)行非線性最小二乘估4.10.2.2-1y(x3e04x12e32x。x在[0,4]中取值；y120.3*(rand(n,1)-0.5)yaea3xaea4x12fminsearch從受污染數(shù)據(jù)中，估計出參數(shù)aa(1a(2a(3a(4a1a2a3a4functionfunction編寫M文件作為主文a0=[1111];Q=1-%char_y_est=[ych,a1,'*exp(-',a3,'*x)+',a2,'*exp(-',a4,'*x)'];plot(x,y,'b+');holdon,plot(x,y_est,'r');holdoff,axis([0,4,0,16])text(0.4,12,char_y_est),text(2.5,9,['chi2=',num2str(chi_est)])text(2.5,7,['freedom=',num2str(freedom)])text(2.5,5,['Q=',exm041022y=3*exp(-0.4*x)+12*exp(-yest=2.8967*exp(-0.38045*x)+11.983*exp(-86420 1234【例4.10.2.2-2ybea1xbea2x 計b b2，使用fminsearch估計a a2。x,y原始數(shù)據(jù)同上例functiona0=[1%whileifr<0.001;break;end%char_y_est=[ych,b1,'*exp(-',a1,'*x)+',b2,'*exp(-',a2,'*x)'];plot(x,y,'b+');holdon;plot(x,y_est,'r');holdtext(2.5,9,['chi2=',num2str(chi_est)])text(2.5,7,['freedom=',num2str(freedom)])text(2.5,5,['Q=',num2str(Q)])exm041022y=3*exp(-0.4*x)+12*exp(-yest=3.0187*exp(-0.39739*x)+11.8737*exp(-86420 1234非線性最小二乘估計124.10.2.3-14.10.2.2-112從受污染數(shù)據(jù)中，估計

yaea3x

ea4

的參functionE=twoexps3(a,x,y)a0=[1100.21];char_y_est=['yest=',a1,'*exp(-',a3,'*x)+',a2,'*exp(-',a4,'*x)'];disp(['原方程','y=3*exp(-0.4*x)+12*exp(-3.2*x)'])原方程y=3*exp(-0.4*x)+12*exp(-yest=2.8381*exp(-0.37236*x12.0338*exp(-一維插4.11.1-1y0.9510 t_ginput % % % % %<6>disp(['t_nearstt_linear disp([t_nearstt_linear t_nearstt_linear t_zero函數(shù)的插%-------5050--樣條插4.12.1-1】根據(jù)連續(xù)時間函數(shù)w(t)etspline重構(gòu)該連續(xù)函數(shù)，并plot(tt,ww,'b');holdonstem(t,w,'filled','r');holdofferror=110-054.12-y=[-0.51-0.5-10.51-0.5;0.510.5-1-0.510.5];plot(yy(1,:),yy(2,:),'b');holdonplot(y(1,:),y(2,:),'or');hold110------ -01樣條函數(shù)用于數(shù)值積分和微x【例4.12.2-1】對于函ysinx，很容易求S(x)0sinxdx1cosxxycosx%err_yy=err_interr_der% %<2>%DefiniteIntegral=bySpline:byTheory:0.9564Derivative=bySpline:-byTheory:-byDiference:-fnplt(pp,'b-');hold er_pp,'r--');hold210--01234567Fourier快速Fourier變換和逆變換指連續(xù)時間函數(shù)的Fourier展開系數(shù)的積分求取Fourier級數(shù)與DFT之間的數(shù)算法實現(xiàn)t 0.5t)

function%symstttifnargin<2;Nf=6;endifnargin<3;Nn=32;endfor%function%A_sym= B_sym function%%ifif(nargin<3|isempty(Nf));Nf=15;endforn=A(n+1)=B(n+1)=quadl(@sin_y,a,a+T,[],[],n,T)/T*2;functionwcos=cos_y(%functionwsin=sin_y(%function%%y=zeros(size(t));ii=find(t>=0.5&t<=1.5); %<1>A_quad=B_quad0asq*0bsq*view([46,38]),colormap(jet),shadingflat,light,lightinggouraud4.13.2.3-3FFT，按式(4.13.2.2-2)和(4.13.2.2-3)Fourier級數(shù)function%if(nargin<2|if %% n=0:1:(N- %% % C=[cn_SHcn_SH(1)]; ifNf>N/2;error([Nfint2str(N/2-1)]);endast=0bst利用DFT計算一般連續(xù)函數(shù)的Fourier變換CFT與DFT之間的數(shù)算法實現(xiàn)4.13.3.2-1FFT求取矩形脈沖w(t)

0t

functionif%%% %%n=0:1:(N-%if%%%if n=0:N-1; TimeTime1002468t---4.13-title('AliasingcausedbyAliasingcausedbybyAliasingcausedbybyFFT10- - - -0f---1124.13-初值常微分方程的解算指ODE的基本機solver解算指令的使用說ODE解算指令的使用解算指令簡潔格式使用示【例4.14.2.1-1】采用最簡潔格式的ODE文件和解算指令，研究圍繞地球旋轉(zhuǎn)的軌道function%%globalG%%Yd=[Vxy;-%globalG plot(X,Y,'b','Linewidth',2);holdon [XE,YE,ZE]=sphere(10); mesh(XE,YE,ZE),hold xx86420----05x4.14-一般說來，編寫程序時，應(yīng)盡量少用全局變量，以免引起。本例演示參數(shù)如何在指令間function%%switchcase error(['Unknownflag'''flag'''.']);解算指令較復(fù)雜格式的使用function%%%%%%%switchflagcase'' varargout{1}=f(t,Y,G,ME); case'init' [varargout{1:3}]=fi(tspan,Y0);%case'events'%[varargout{1:3}]=fev(t,Y,Y0);error(['Unknownflag'''flag'''.']);%functionYd=%X=Y(1:2);V=Y(3:4);r=sqrt(sum(X.^2));Yd=[V;-%function[ts,y0,options]=%ts=tspan;y0=options=odeset('Events','on','Reltol',1e-5,'Abstol',1e-%%function[value,isterminal,direction]=%dDSQdt=2*((Y(1:2)-Y0(1:2))'*value=[dDSQdt;dDSQdt];%direction=[1;-1];%isterminal=[1;0]; %%%[XE,YE,ZE]=sphere(10);RE=0.64e7;XE=RE*XE;YE=RE*YE;ZE=0*ZE;關(guān)于ODE文件的關(guān)于解算指令選項optionsoptionsoptions屬性處理和輸出函數(shù)使用演【例4.14.4.2-1】仍以軌道問題為例（原題見例4.14.2.1-1,4.14.2.1-2,4.14.2.2-1）。本例T具體目標(biāo)是：畫出解向量Y y4 v中由x,v構(gòu)成的T xfunctionvarargout=DYDt4(t,Y,flag,G,ME,tspan,Y0)switchflagcasevarargout{1}=f(t,Y,G,ME);case'init'[varargout{1:3}]=fi(tspan,Y0);error(['Unknownflag'''flag'''.']);%functionYd=f(t,Y,G,ME)X=Y(1:2);V=Y(3:4);r=sqrt(sum(X.^2));Yd=[V;-G*ME*X/r^3];%function[ts,y0,options]=fi(tspan,Y0)ts=tspan;y0=Y0;%options.OutputSel=[13];%clf,set(gca,'xlim',[-525]*1e7,'ylim',[-33]*1e3);% boxonholdon;ode45('DYDt4',[],[],[],G,ME,tspan,Y0);holdoff常微分方程的邊值問題bvp4c求解邊值問題的基本思求解邊值問題的基本配套指【例4.14.5.2-1zcz0z(0)0z(4)2c=1。本例的目的：（A）完整演示解題步驟；（B）functiondydx=twoode(x,y,c)dydx=[y(2),-funct