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

第四章數(shù)值計(jì) 計(jì)算指令之間的內(nèi)在聯(lián)系及區(qū)別。對(duì)于那些熟悉其他高級(jí)語(yǔ)言（FORTRAN，Pascal，C++）的讀者來(lái)說(shuō)，通過(guò)本章，對(duì)于那些經(jīng)過(guò)大學(xué)基本數(shù)學(xué)的讀者來(lái)說(shuō)，通過(guò)本章，精良完善的計(jì)算指令，自然易讀的程序?qū)⑹顾麄兏形颉啊睌?shù)學(xué)的基礎(chǔ)地位和局限性，看到從“理想化”簡(jiǎn)對(duì)于那些熟悉基本指令的讀者來(lái)說(shuō)，通過(guò)本章，圍繞基本數(shù)值問(wèn)題展開(kāi)的由于的基本運(yùn)算單元是數(shù)組，所以本章內(nèi)容將從矩陣分析、線性代數(shù)的數(shù)值Fourier分析，和一般常微分方程初值、邊值問(wèn)題。本章的最后稀疏矩陣的處理，因?yàn)?6.0LINPACKEISPACK庫(kù)為基LAPACKARPACK程序庫(kù)的新基礎(chǔ)上。為構(gòu)成空間的基向量時(shí)不唯一的，且新版的更。本書(shū)新版全部算例結(jié)果是在6.x版 6.0版起，提供了一組專門求微分方程“邊值問(wèn)題”數(shù)值解的指令。適應(yīng)4.14.52個(gè)算例闡述求解細(xì)節(jié)。triplequad。本章新版對(duì)此作了相應(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('實(shí)際誤差disp('近似的最大可能誤差er_approx'),disp(er_approx),disp('')disp('最大可能誤差er_max'),disp(er_max),disp('') 實(shí)際誤差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=特征值問(wèn)題的條件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ù)特征值對(duì)角陣與實(shí)數(shù)塊特征值對(duì)角陣的轉(zhuǎn)4.3.3-14.3.1-1中的復(fù)數(shù)特征值對(duì)角陣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)的奇異值分解描線性二乘問(wèn)題矩陣除運(yùn)算的廣義線性模型的最小二乘【例4.4.2.2-1】對(duì)于超定方程yAx，進(jìn)行三種解法比較。其中A取庫(kù)中的特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ù)值計(jì)算指4.5.2.2-1】用一個(gè)簡(jiǎn)單矩陣表現(xiàn)diffgradientF=123456789Dx333333=111111FX111111111FY333333333=222222222=666666666【例4.5.2.2-2】研究偶極子(Dipole)的電勢(shì)（Electricpotential）和電場(chǎng)強(qiáng)度（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{隸書(shū)}\fontsize{22}偶極子的場(chǎng)'),holdonxlabel('x');ylabel('y'),hold偶極子的偶極子的6420---- - -0x 6多項(xiàng)式一元函數(shù)的零利 作圖指令獲取初步近似任意一元函數(shù)零點(diǎn)的精確【例4.6.2.2-1】通過(guò)求f(t)(sin2t)eatbt的零點(diǎn)，綜合敘述相關(guān)指令的用法 10--10-------0t54.6-令必須在指令窗中運(yùn)行，并得到如圖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-開(kāi)型數(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=多重?cái)?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ù)值卷多項(xiàng)式乘法與離散卷積的算法同連續(xù)時(shí)間函數(shù)的數(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）符號(hào)Laplace變換求卷積的理論表示；（B） symstao;t=sym('t','positive');HS1=yt1=4.8-dy124.8.3.5-3】用“零階”近似法求u(t)etU(t和h(t)tet2U(t點(diǎn)的確定。（D）離散卷積和連續(xù)卷積之間的關(guān)系。（E）conv的使用。（F）繪圖分% [3,11]段[11,14]段']),disp(dy13)003 003統(tǒng)計(jì)分布的數(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ù)和隨機(jī)數(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{隸書(shū)}\fontsize{22}置信水平0.9') holdoff 置信水平005多項(xiàng)式多項(xiàng)式擬合原擬合多項(xiàng)式階數(shù)的確多項(xiàng)式擬合 實(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ù)中，估計(jì)出參數(shù)aa(1a(2a(3a(4a1a2a3a4functionfunction編寫(xiě)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 計(jì)b b2，使用fminsearch估計(jì)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非線性最小二乘估計(jì)124.10.2.3-14.10.2.2-112從受污染數(shù)據(jù)中，估計(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í)間函數(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】對(duì)于函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í)間函數(shù)的Fourier展開(kāi)系數(shù)的積分求取Fourier級(jí)數(shù)與DFT之間的數(shù)算法實(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級(jí)數(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計(jì)算一般連續(xù)函數(shù)的Fourier變換CFT與DFT之間的數(shù)算法實(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的基本機(jī)solver解算指令的使用說(shuō)ODE解算指令的使用解算指令簡(jiǎn)潔格式使用示【例4.14.2.1-1】采用最簡(jiǎn)潔格式的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-一般說(shuō)來(lái)，編寫(xiě)程序時(shí)，應(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)于解算指令選項(xiàng)optionsoptionsoptions屬性處理和輸出函數(shù)使用演【例4.14.4.2-1】仍以軌道問(wèn)題為例（原題見(jiàn)例4.14.2.1-1,4.14.2.1-2,4.14.2.2-1）。本例T具體目標(biāo)是：畫(huà)出解向量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常微分方程的邊值問(wèn)題bvp4c求解邊值問(wèn)題的基本思求解邊值問(wèn)題的基本配套指【例4.14.5.2-1zcz0z(0)0z(4)2c=1。本例的目的：（A）完整演示解題步驟；（B）functiondydx=twoode(x,y,c)dydx=[y(2),-funct