山大MATLAB編程指導(dǎo)課件第4章 MATLAB 的數(shù)值計算

1、2022/8/19第1頁第4章 MATLAB 的數(shù)值計算 4.1 數(shù)值微積分4.2 矩陣和代數(shù)方程4.3 概率分布和統(tǒng)計分析4.4 多項式運算和卷積2022/8/19第2頁4.1數(shù)值微積分4.1.1 近似數(shù)值極限及導(dǎo)數(shù)4.1.2 數(shù)值求和與近似數(shù)值積分4.1.3 計算精度可控的數(shù)值積分4.1.4 函數(shù)極值的數(shù)值求解4.1.5 常微分方程的數(shù)值解2022/8/19第3頁在MATLAB數(shù)值計算中，既沒有專門的求極限指令，也沒有專門的求導(dǎo)指令。但MATLAB提供了與“求導(dǎo)”概念有關(guān)的“求差分”指令。dx=diff(X) %計算向量X的前向差分FX=gradient(F) %求一元(函數(shù))梯度FX,

2、FY =gradient(F) %求二元（函數(shù)）梯度對diff而言，當X是向量時，dx= X(2:n)-X(1:n-1) ;當X是矩陣時，dx= X(2:n, :)-X(1:n-1, :) 。 dx的長度比x的長度少1個元素。4.1.1 近似數(shù)值極限及導(dǎo)數(shù)2022/8/19第4頁dx=diff(X) %計算向量X的前向差分FX=gradient(F) %求一元(函數(shù))梯度FX, FY =gradient(F) %求二元（函數(shù)）梯度對gradient而言，當F是向量時，F(xiàn)X(1) = F(2)-F(1), FX(2:end-1) = (F(3:end)-F(1:end-2)/2 , FX(end

3、) = F(end)-F(end-1) ; FX長度與F的長度相同當F是矩陣時， FX, FY是與F同樣大小的矩陣。 FX的每行對應(yīng)F相應(yīng)行元素間的梯度 ; FY的每列對應(yīng)F相應(yīng)列元素間的梯度 ; 4.1.1 近似數(shù)值極限及導(dǎo)數(shù)2022/8/19第5頁數(shù)值極限和導(dǎo)數(shù)的應(yīng)用應(yīng)十分謹慎x=eps;L1=(1-cos(2*x)/(x*sin(x), L2=sin(x)/x,L1 = 0L2 = 1syms tf1=(1-cos(2*t)/(t*sin(t);f2=sin(t)/t;Ls1=limit(f1,t,0)Ls2=limit(f2,t,0) Ls1 = 2Ls2 = 1 x=pi/1000;

4、 %可得到正確結(jié)果2022/8/19第6頁數(shù)值極限和導(dǎo)數(shù)的應(yīng)用應(yīng)十分謹慎%（1）自變量的增量取得過?。╡ps數(shù)量級）d=pi/100;t=0:d:2*pi; x=sin(t);dt=5*eps; x_eps=sin(t+dt);dxdt_eps=(x_eps-x)/dt;plot(t,x,LineWidth,5)hold onplot(t,dxdt_eps)hold offlegend(x(t),dx/dt)xlabel(t) 【例4.1-2】已知, 求該函數(shù)在區(qū)間 中的近似導(dǎo)函數(shù)。數(shù)值導(dǎo)數(shù)受計算中有限精度困擾，當增量dt過小時，f(t+dt)與f(t)的數(shù)值十分接近，高位有效數(shù)字完全相同，

5、df =f(t+dt)-f(t) 造成高位有效數(shù)字消失，精度急劇變差。2022/8/19第7頁數(shù)值極限和導(dǎo)數(shù)的應(yīng)用應(yīng)十分謹慎%（2）自變量的增量取得適當x_d=sin(t+d);dxdt_d=(x_d-x)/d;plot(t,x,LineWidth,5)hold onplot(t,dxdt_d)hold offlegend(x(t),dx/dt)xlabel(t) 【例4.1-2】已知, 求該函數(shù)在區(qū)間 中的近似導(dǎo)函數(shù)。d=pi/100;2022/8/19第8頁d=pi/100; t=0:d:2*pi;x=sin(t);dxdt_diff=diff(x)/d;dxdt_grad=gradien

6、t(x)/d;【例4.1-3】已知 采用diff和gradient計算該函數(shù)在區(qū)間 中的近似導(dǎo)函數(shù)。subplot(1,2,1); plot(t,x,b);hold onplot(t,dxdt_grad,m,LineWidth,8)plot(t(1:end-1),dxdt_diff,.k,MarkerSize,8)axis(0,2*pi,-1.1,1.1); title(0, 2pi)legend(x(t),dxdt_grad,dxdt_diff,Location,North)xlabel(t), box off;hold offsubplot(1,2,2)kk=(length(t)-10):

7、length(t); hold on; plot(t(kk),dxdt_grad(kk),om,MarkerSize,8)plot(t(kk-1),dxdt_diff(kk-1),.k,MarkerSize,8)title(end-10, end)legend(dxdt_grad,dxdt_diff,Location,SouthEast)xlabel(t),box off; hold off 宏觀上, diff和gradient結(jié)果大致相同細節(jié)上, diff和gradient數(shù)值有差異， diff沒有給出最后一點導(dǎo)數(shù)2022/8/19第9頁2022/8/19第10頁Sx=sum(X) % su

8、m按列向求和得(1n)數(shù)組Scs=cumsum(X)%沿X列向求累計和, 仍是數(shù)組, 第(i, k)個元素是X數(shù)組第k列前i個元素的和。最后一行等于Sx4.1.2 數(shù)值求和與近似數(shù)值積分St=trapz(t,X) 或 St=dt*trapz(X) %梯形法求積分Sct=cumtrapz(t,X) 或 Sct=dt*cumtrapz(X) %梯形法沿列向求X關(guān)于x的累計積分,最后一個值等于StS=dt*sum(X), S=sum(t,X) %近似矩形法求積分Scs=cumsum(t,X) =dt*cumsum(X) 2022/8/19第11頁clear; d=pi/8; t=0:d:pi/2;

9、y=0.2+sin(t);s=sum(y); s_sa=d*s;% s_sa=sum(t, y), 近似矩形法積分s_ta=trapz(t,y); %梯形法積分【例4.1-4】求積分 , 其中disp(sum求得積分,blanks(3),trapz求得積分)disp(s_sa, s_ta)t2=t,t(end)+d; y2=y,nan;stairs(t2,y2,:k);hold onplot(t,y,r,LineWidth,3)h=stem(t,y,LineWidth,2);set(h(1),MarkerSize,10)axis(0,pi/2+d,0,1.5)hold off; shg sum

10、求得積分 trapz求得積分 1.5762 1.30132022/8/19第12頁4.1.3 計算精度可控的數(shù)值積分一重積分(quadrature精度可控)：S1=quad(fun,a,b,tol) %自適應(yīng)Simpson法S2=quadl(fun,a,b,tol) %自適應(yīng)羅巴托 Lobatlo法二重積分(精度可控) ：S3=dblquad(fun,xmin,xmax,ymin,ymax,tol)三重積分(精度可控) ： S4=triplequad(fun,xmin,xmax,ymin,ymax, zmix,zmax,tol) fun:可為字符串,內(nèi)聯(lián)對象,匿名函數(shù),函數(shù)句柄,注意乘除法運算

11、一定加點.用數(shù)組運算 tol: 默認積分的絕對精度為10-62022/8/19第13頁（1）syms xIsym=vpa(int(exp(-x2),x,0,1) % x.2數(shù)組乘方亦可Isym = 0.74682413281242702539946743613185 例4.1-5：求 . % （2） 梯形法積分format longd=0.001;x=0:d:1;Itrapz=d*trapz(exp(-x.*x) % x.*x必須為數(shù)組乘Itrapz = 0.746824071499185 （3）fx=exp(-x.2); %一定用數(shù)組乘Ic=quad(fx,0,1,1e-8) Ic = 0.

12、746824132854452 %實際精度控制到1e-102022/8/19第14頁（1）符號計算法syms x ys=vpa(int(int(xy,x,0,1),y,1,2)Warning: Explicit integral could not be found. s = 0.40546510810816438197801311546435 例4.1-6：求 . （2）數(shù)值積分法format longs_n=dblquad(x,y)x.y,0,1,1,2) %匿名函數(shù)表示被積函數(shù)s_n = 0.405466267243508 s_n=dblquad(x.y,0,1,1,2) %字符串表示被

13、積函數(shù)s_n=dblquad(inline(x.y),0,1,1,2) %內(nèi)聯(lián)函數(shù)表示被積函數(shù)% 一定為數(shù)組乘2022/8/19第15頁4.1.4 函數(shù)極值的數(shù)值求解 x,fval,exitflag,output=fminbnd(fun,x1,x2,options) %一元函數(shù)在區(qū)間bound(x1,x2)中極小值 x,fval: 極值點坐標和對應(yīng)目標函數(shù)極值x,fval,exitflag,output=fminsearch(fun,x0,options) %單純形法求搜索起點x0附近多元函數(shù)極值點 % x每列代表一個候選極值點，各列按目標函數(shù)極小值遞增順序 %x(:,1)對應(yīng)的目標函數(shù)極小值

14、點由fval給出 fun: 字符串,內(nèi)聯(lián)函數(shù),匿名函數(shù),函數(shù)句柄,注意乘除法運算一定加點.用數(shù)組運算 options: 配置優(yōu)化參數(shù),可略 exitflag: 給出大于0的數(shù),則成功搜索到極值點 output: 給出具體的優(yōu)化算法和迭代次數(shù)2022/8/19第16頁【例4.1-7 】已知 ，在-10 x10區(qū)間，求函數(shù)的最小值。（1）用“導(dǎo)數(shù)為零法”求極值點syms xy=sin(x)2*exp(-0.1*x)-0.5*sin(x)*(x+0.1);yd=diff(y,x); %求導(dǎo)函數(shù)xs0=solve(yd,x) %求導(dǎo)函數(shù)為零的根，即極值點yd_xs0=vpa(subs(yd,x,xs0

15、) %驗證極值點處導(dǎo)函數(shù)為零y_xs0=vpa(subs(y,x,xs0) %求極值點處極值xs0 = matrix(0.050838341410271656880659496266968)yd_xs0 = 2.2958874039497802890014385492622*10(-41)y_xs0 = -0.001263317776974196724544154344118 無法判斷是否最小值2022/8/19第17頁【例4.1-7 】已知 ，在-10 x10區(qū)間，求函數(shù)的最小值。（2）采用優(yōu)化算法求極小值x1=-10;x2=10;yx=(x)(sin(x)2*exp(-0.1*x)-0.5

16、*sin(x)*(x+0.1); xn0,fval,exitflag,output=fminbnd(yx,x1,x2)xn0 = 2.514797840754235fval = %比“導(dǎo)數(shù)為零法”求得的極值更小, 更可能是最小值 -0.499312445280039exitflag = 1output = iterations: 13 funcCount: 14 algorithm: golden section search, parabolic interpolation message: 1x112 char 2022/8/19第18頁【例4.1-7 】已知 ，在-10 x10區(qū)間，求函

17、數(shù)的最小值。（4）據(jù)圖形觀察，重設(shè)fminbnd的搜索區(qū)間x11=6;x2=10;yx=(x)(sin(x)2*exp(-0.1*x)-0.5*sin(x)*(x+0.1); xn00,fval,exitflag,output=fminbnd(yx,x11,x2)xn00 = 8.023562824723015fval = -3.568014059128578 %最小值exitflag = 1output = iterations: 9 funcCount: 10 algorithm: golden section search, parabolic interpolation message

18、: 1x112 char y=sin(x)2*exp(-0.1*x)-0.5*sin(x)*(x+0.1);（3）繪圖觀察最小值xx=-10:pi/200:10;yxx=subs(y,x,xx);plot(xx,yxx)xlabel(x),grid on 2022/8/19第19頁2022/8/19第20頁例4.1-8: f(x,y)=100(y-x2)2+(1-x)2在區(qū)間-5,5的極小值ff=(x)100*(x(2)-x(1).2).2+(1-x(1).2; %匿名函數(shù) x0=-5,-2,2,5;-5,-2,2,5; %4個搜索起點 sx,sfval,sexit,soutput=fmins

19、earch(ff,x0)其理論極小值點為x=1, y=1%sx給出一組使優(yōu)化函數(shù)值非減的局部極小點 sx = 0.99998 -0.68971 0.41507 8.0886 0.99997 -1.9168 4.9643 7.8004sfval = 2.4112e-0102.4112e-010 5.7525e+002 2.2967e+003 3.3211e+005 format short e %取5位科學(xué)計數(shù)法 disp(ff(sx(:,1),ff(sx(:,2),ff(sx(:,3),ff(sx(:,4) 用x的二元向量表示x,yx0 =-5 -2 2 5 -5 -2 2 52022/8/1

20、9第21頁4.1.5 常微分方程Ordinary Differential Equation的數(shù)值解t,Y=ode45(odefun,tspan,y0) % 4階龍格庫塔 數(shù)值法odefun: 待求解一階微分方程組的函數(shù)文件句柄tspan: 自變量微分二元區(qū)間t0, tfy0: 一階微分方程組的(n*1)初值列向量matlab為解常微分方程初值問題提供了一組配套齊全,結(jié)構(gòu)嚴整的指令, 包括: ode45, ode23, ode113, ode23t, ode15s, ode23s, ode23tb.在此只介紹最常用的ode45的基本使用方法.ode45使用方法:t二元區(qū)間的點系列Y原函數(shù)在微分

21、區(qū)間點系列上的函數(shù)值2022/8/19第22頁例4.1-9求解：%解算微分方程tspan=0,30; y0=1;0;tt,yy=ode45(DyDt,tspan,y0); figure(1)plot(tt,yy(:,1)xlabel(t),title(x(t)解:令，上式寫成一階微分方程組形式%畫相平面圖(函數(shù)和其導(dǎo)數(shù)勾畫的曲線稱為相軌跡)figure(2)plot(yy(:,1),yy(:,2)xlabel(位移),ylabel(速度) function ydot=DyDt(t,y)mu=2; ydot=y(2);mu*(1-y(1)2)*y(2)-y(1);據(jù)以上方程組，編寫M函數(shù)文件Dy

22、Dt.m2022/8/19第23頁2022/8/19第24頁4.2矩陣和代數(shù)方程4.2.1 矩陣運算和特征參數(shù)4.2.2 矩陣的變換和特征值分解4.2.3 線性方程的解4.2.4 一般代數(shù)方程的解2022/8/19第25頁4.2.1 矩陣運算和特征參數(shù)矩陣與標量之間的四則運算與數(shù)組運算相同矩陣和矩陣之間的四則運算矩陣和矩陣之間的加減運算與數(shù)組運算相同設(shè) A 是一個 mn 矩陣，B 是一個 pq 矩陣，當 np 時，兩個矩陣可以相乘，乘積為 mq 矩陣。矩陣乘法不可逆。在 MATLAB 中，矩陣乘法由“*”實現(xiàn)。矩陣除法在實際中主要用于求解線性方程組矩陣轉(zhuǎn)置：符號“ ”實現(xiàn)矩陣的轉(zhuǎn)置操作。對于實

23、數(shù)矩陣， “ ”表示矩陣轉(zhuǎn)置，對于復(fù)數(shù)矩陣，“ ”實現(xiàn)共軛轉(zhuǎn)置。對于復(fù)數(shù)矩陣，如果想要實現(xiàn)非共軛轉(zhuǎn)置，可以使用符號“ . ”。2022/8/19第26頁format rat %有理格式顯示A=magic(2) + j*pascal(2) A = 1 + 1i 3 + 1i 4 + 1i 2 + 2iA1=A A1=1 - 1i 4 - 1i %共軛轉(zhuǎn)置 3 - 1i 2 - 2iA2=A. A2=1 + 1i 4 + 1i %非共軛轉(zhuǎn)置，數(shù)組運算操作 3 + 1i 2 + 2i例4.2-2 矩陣和數(shù)組轉(zhuǎn)置操作的差別2022/8/19第27頁計算矩陣標量特征參數(shù)秩,跡,行列式的指令rank(A)

24、 %求秩（Rank）det(A) %求行列式（Determinant）trace(A) %求跡（Trace），即矩陣主對角元素的和2022/8/19第28頁 A=reshape(1:9,3,3); r=rank(A) %求秩 d3=det(A) %非滿秩矩陣的行列式一定為零 d2=det(A(1:2,1:2) %求子式的行列式 t=trace(A)【例4.2-3】矩陣標量特征參數(shù)計算示例。A = 1 4 7 2 5 8 3 6 9r = 2d3 =0d2 = -3t = 15 2022/8/19第29頁【例4.2-4】矩陣標量特征參數(shù)的性質(zhì)。format short; rand(twister

25、,0)A=rand(3,3); B=rand(3,3);C=rand(3,4); D=rand(4,3);tAB=trace(A*B) %任何符合矩陣乘法規(guī)則的兩個矩陣乘積tBA=trace(B*A) %的“跡”不變。同階乘積“跡”不變tCD=trace(C*D) %兩個“內(nèi)維”相等矩陣的乘積“跡”不變tDC=trace(D*C)tAB = 2.6030tBA = 2.6030tCD = 4.1191tDC = 4.1191 dCD=det(C*D)dDC=det(D*C) dCD = 0.0424 dDC = -2.6800e-018d_A_B=det(A)*det(B)dAB=det(A*

26、B) dBA=det(B*A) 同階矩陣乘積行列式等于各矩陣行列式之乘積d_A_B = 0.0094dAB = 0.0094dBA = 0.0094 非同階矩陣乘積行列式不等于各矩陣行列式之乘積2022/8/19第30頁dx=diff(X) %計算向量X的前向差分FX=gradient(F) %求一元(函數(shù))梯度FX, FY =gradient(F) %求二元（函數(shù)）梯度4.1.1 近似數(shù)值極限及導(dǎo)數(shù)S1=sum(X,1) % sum按列向求和得(1n)數(shù)組, =sum(X) (X多行)S2=sum(X,2) % sum按行向求和得(n1)數(shù)組Scs=cumsum(X)%沿X列向求累計和, 仍

27、是數(shù)組, 最后一行等于Sx4.1.2 數(shù)值求和與近似數(shù)值積分St=trapz(t,X) 或 St=dt*trapz(X) %梯形法求積分Sct=cumtrapz(t,X) 或 Sct=dt*cumtrapz(X) %梯形法沿列向求X關(guān)于x的累計積分,最后一個值等于StX= 1 2; 3 4 S1= 4 6 S2= 3 7sum(S1) = 10sum(S1,1) = 4 6Review2022/8/19第31頁4.1.3 計算精度可控的數(shù)值積分一重積分(quadrature精度可控)：S1=quad(fun,a,b,tol) %自適應(yīng)Simpson法S2=quadl(fun,a,b,tol)

28、%自適應(yīng)羅巴托 Lobatlo法二重積分(精度可控) ：S3=dblquad(fun,xmin,xmax,ymin,ymax,tol)三重積分(精度可控) ： S4=triplequad(fun,xmin,xmax,ymin,ymax, zmix,zmax,tol) fun:可為字符串,內(nèi)聯(lián)對象,匿名函數(shù),函數(shù)句柄,注意乘除法運算一定加點.用數(shù)組運算 tol: 默認積分的絕對精度為10-6Review2022/8/19第32頁4.1.4 函數(shù)極值的數(shù)值求解 x,fval,exitflag,output=fminbnd(fun,x1,x2,options) %一元函數(shù)在區(qū)間bound(x1,x2

29、)中極小值 x,fval: 極值點坐標和對應(yīng)目標函數(shù)極值x,fval,exitflag,output=fminsearch(fun,x0,options) %單純形法求搜索起點x0附近多元函數(shù)極值點 % x每列代表一個候選極值點，各列按目標函數(shù)極小值遞增順序 %x(:,1)對應(yīng)的目標函數(shù)極小值點由fval給出 fun: 字符串,內(nèi)聯(lián)函數(shù),匿名函數(shù),函數(shù)句柄, 注意乘除法運算一定加點.用數(shù)組運算 options: 配置優(yōu)化參數(shù),可略 exitflag: 給出大于0的數(shù),則成功搜索到極值點 output: 給出具體的優(yōu)化算法和迭代次數(shù)Review2022/8/19第33頁rank(A) %求秩（R

30、ank）det(A) %求行列式（Determinant）trace(A) %求跡(Trace)，即矩陣主對角元素的和Reviewt,Y=ode45(odefun,tspan,y0) % 4階龍格庫塔 數(shù)值法odefun: 待求解一階微分方程組的函數(shù)文件句柄tspan: 自變量微分二元區(qū)間t0, tfy0: 一階微分方程組的(n*1)初值列向量ode45使用方法:t二元區(qū)間的點系列Y原函數(shù)在微分區(qū)間點系列上的函數(shù)值4.1.5 常微分方程Ordinary Differential Equation的數(shù)值解4.2.1 矩陣運算和特征參數(shù)2022/8/19第34頁4.2.2 矩陣的變換和特征值分解R

31、, ci=rref(A) %借助初等變換將A縮減行變成行階梯矩陣R。B=orth(A)可得到矩陣A的正交基，B的列與A的列可張成相同的空間，而且B的列是正交的，因此B*B=eye(rank(A)，B的列數(shù)正好是A的秩。 X=null(A) %A矩陣零空間的全部正交基，滿足AX=0, X*X=I。B=orth(A)V, D=eig(A) %A矩陣的特征值、特征向量分解，使AV=VD。ci 是行數(shù)組, 其元素表示A中線性獨立列的序號。length(ci)=rank(A)2022/8/19第35頁【例4.2-5】行階梯陣簡化指令rref計算結(jié)果的含義。A=magic(4)R,ci=rref(A) %

32、行階梯分解A = 16 2 3 13 5 11 10 8 9 7 6 12 4 14 15 1R = 1 0 0 1 0 1 0 3 0 0 1 -3 0 0 0 0ci = 1 2 3 r_A=length(ci) % =rank(A)r_A = 3 aa=A(:,1:3)*R(1:3,4) %A前三列線形組合err=norm(A(:,4)-aa) aa = 13 8 12 1err = 0 2022/8/19第36頁A=reshape(1:15,5,3) X=null(A)S=A*XS =1.0e-014 * 0 -0.1776 -0.2665 -0.3553 -0.5329A =1 6

33、11 2 7 12 3 8 13 4 9 14 5 10 15n= 3l= 1Rank_A=2X =0.4082 -0.8165 0.4082ans = 1 例4.2-6 矩陣零空間及含義。設(shè) 是矩陣 的零空間,即n=size(A,2)l=size(X,2)Rank_A=rank(A)n-l=rank(A)2022/8/19第37頁【例4.2-7】簡單實陣的特征分解。eig，cdf2rdf Complex Diagonal Form to Real-block Diagonal Form%（1）A=1,-3;2,2/3V,D=eig(A) A =1.0000 -3.0000 2.0000 0.

34、6667V = 0.7746 0.7746 0.0430 - 0.6310i 0.0430 + 0.6310iD =0.8333 + 2.4438i 0 0 0.8333 - 2.4438i %（2）VR,DR=cdf2rdf(V,D) VR =0.7746 0 0.0430 -0.6310DR =0.8333 2.4438 -2.4438 0.8333 %（3）A1=V*D/V %因計算誤差有很小虛部A1_1=real(A1) %去除虛部后等于AA2=VR*DR/VRerr1=norm(A-A1,fro)err2=norm(A-A2,fro) A1 =1.0000 + 0.0000i -3.

35、0000 2.0000 - 0.0000i 0.6667 A1_1 =1.0000 -3.0000 2.0000 0.6667A2 =1.0000 -3.0000 2.0000 0.6667err1 =6.7532e-016err2 =4.4409e-016 Frobenius矩陣范數(shù) 實數(shù)塊對角陣實陣矩陣Frobenius范數(shù),類似向量2范數(shù)，不同于矩陣2范數(shù)2022/8/19第38頁4.2.2 線性方程的解對于方程組Amnx=b (m個方程，n個未知數(shù))， 當向量b在矩陣A列向量所張空間中, rank(A,b) =rank(A)= r a) 若n =r，則解唯一。 b) 若n r，則解不唯

36、一。 當向量b不在矩陣A列向量所張空間中, 則無準確解但存在最小二乘解。1. 線性方程解的一般討論matlab定義的左除運算可以很方便地解上述方程組： x=Ab (x=A-1b只適用于A非奇異時)2.除法運算解方程解 當m=n時, “恰定”方程 當mn時, “超定”方程 當mr, 解不唯一% （3）求特解和通解，并對由之構(gòu)成的全解進行驗算xs=Ab; xg=null(A); % xg是齊次方程Ax=0的解c=rand(1); ba=A*(xs+c*xg) % ba是A乘 “一個隨機的全解”Warning: Rank deficient, rank = 2, tol = 1.8757e-014.

37、ba =13.0000 14.0000 15.0000 16.0000 例4.2-8: 求方程 的解norm(ba-b) ans = 1.1374e-014 2022/8/19第40頁例4.2-9: “逆陣法”和”左除法”解恰定方程的性能對比。測試陣產(chǎn)生指定異常值和特殊帶寬的隨機陣3. 矩陣的逆A_1=inv(A) %求非奇異陣A的逆randn(state,0); A=gallery(randsvd,300,2e13); %產(chǎn)生條件數(shù) 為2e13的300階隨機矩陣;x=ones(300,1); %定義真解b=A*x;cond(A) %驗算矩陣條件數(shù), 結(jié)果a1.9978e+013, 值越大,

38、陣越病態(tài),%求逆法tic; xi=inv(A)*b; ti=toceri=norm(x-xi) rei=norm(A*xi-b)/norm(b) %左除法tic; xd=Ab;td=tocerd=norm(x-xd)red=norm(A*xd-b)/norm(b)ti =0.0185eri =0.0883rei =0.0051td =0.0066erd =0.0298red =8.7810e-015矩陣計算對于誤差越敏感, 數(shù)值穩(wěn)定性差 2022/8/19第41頁求解任意函數(shù)f(x)=0(可能無解,單解,多解)的步驟:作圖獲取初步近似解 觀察f(x)與橫軸的交點坐標，用zoom放大，用ginp

39、ut得較精確些的交點坐標值。4.3 一般代數(shù)方程的解 用”泛函”指令求精確解 x,favl = fzero(fun,x0) %一元函數(shù)求零點 x,fval = fsolve(fun,x0) %解非線性方程組 fun: 內(nèi)聯(lián)對象,匿名函數(shù),函數(shù)句柄,字符串;被解函數(shù)自變量一般用x, 注意乘除法運算一定加點.用數(shù)組運算 x0: 零點初始猜測值. x0為標量時取與之最靠近的零點; x0取a,b時,在此區(qū)間內(nèi)尋找一個零點. x: 所求零點的自變量值 fval: 函數(shù)值2022/8/19第42頁例4.2-10: 求f(t)=(sin2t)e-0.1t0.5|t|的零點y_C=inline(sin(t).

40、2.*exp(-0.1*t)-0.5*abs(t),t); t=-10:0.01:10;Y=y_C(t); clf, plot(t,Y,r); hold onplot(t,zeros(size(t),k);xlabel(t);ylabel(y(t) t4 = 0.5993y4 = 0tt =-2.0039 -0.5184 -0.0042 0.6052 1.6717 t3 = -0.5198y3 = 5.5511e-017zoom ontt,yy=ginput(5),zoom offt4,y4=fzero(y_C,tt(4) t3,y3=fzero(y_C,tt(3) %t=0處沒穿越橫軸202

41、2/8/19第43頁4.3 概率分布和統(tǒng)計分析4.3.1 概率函數(shù)、分布函數(shù)、逆分布函數(shù)和隨機數(shù)的發(fā)生4.3.2 隨機數(shù)發(fā)生器和 統(tǒng)計分析指令2022/8/19第44頁4.3 概率分布和統(tǒng)計分析1. 二項分布（Binomial distribution ）4.3.1 概率函數(shù)、分布函數(shù)、逆分布函數(shù)和隨機數(shù)的發(fā)生pk=binopdf(k,N,p) 事件A發(fā)生k次的概率Fk=binocdf(k,N,p) 事件A發(fā)生次數(shù)不大于k的概率R=binornd(N,p) 產(chǎn)生符合二項分布B(N,p)的(mn) 隨機數(shù)組（元素值為事件A可能發(fā)生的次數(shù)）。Binomial Cumulative Distribu

42、tion FunctionBinomial Probability Density Function2022/8/19第45頁N=100;p=0.5;k=0:N;pdf=binopdf(k,N,p); cdf=binocdf(k,N,p);h=plotyy(k,pdf,k,cdf);set(get(h(1),Children),Color,b,Marker,.,MarkerSize,13)set(get(h(1),Ylabel),String,pdf)set(h(2),Ycolor,1,0,0)set(get(h(2),Children),Color,r,Marker,+,MarkerSize,4)set(get(h(2),Ylabel),String,cdf)xlabel(k) grid on 例4.3-1. 畫出N=100，p=0.5情況下的二項分布概率特性曲線。4.3.1 概率函數(shù)、分布函數(shù)、逆分布函數(shù)和隨機數(shù)的發(fā)生2022/8/19第46頁px=normpdf(x,M