1、Matlab数值计算第 4 章 数值计算 .1.1 近似数值极限及导数【例4.1-1】x=eps;L1=(1-cos(2*x)/(x*sin(x), %L2=sin(x)/x, % L1 = 0L2 = 1 syms 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 【例4.1-2】（1）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,LineW
2、idth,5)hold on plot(t,dxdt_eps)hold offlegend(x(t),dx/dt)xlabel(t) 图 4.1-1 增量过小引起有效数字严重丢失后的毛刺曲线（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 增量适当所得导函数比较光滑【例4.1-3】clfd=pi/100; %t=0:d:2*pi;x=sin(t);dxdt_diff=diff(x)/d; %dxdt_gr
3、ad=gradient(x)/d; %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 offhold offsubplot(1,2,2)kk=(length(t)-10):length(t);%thold onplot(t(kk),dxdt_gra
4、d(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 offhold off 图 4.1-3 diff和gradient求数值近似导数的异同比较.1.2 数值求和与近似数值积分【例 4.1-4】 cleard=pi/8; %t=0:d:pi/2; %y=0.2+sin(t); %s=sum(y); %s_sa=d*s; % s_ta=d*trapz(y); %
5、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 offshg sum求得积分 trapz求得积分 1.5762 1.3013图 4.1-4 sum 和trapz求积模式示意.1.3 计算精度可控的数值积分【例 4.1-5】（1）syms xIs
6、ym=vpa(int(exp(-x2),x,0,1) Isym =0.74682413281242702539946743613185 （2）format long %d=0.001;x=0:d:1;Itrapz=d*trapz(exp(-x.*x) Itrapz = 0.746824071499185 （3）fx=exp(-x.2); %Ic=quad(fx,0,1,1e-8) % Ic = 0.746824132854452 【例 4.1-6】（1）syms x ys=vpa(int(int(xy,x,0,1),y,1,2) % s =0.4054651081081643819780131
7、1546435 （2）format longs_n=dblquad(x,y)x.y,0,1,1,2) % s_n = 0.405466267243508 .1.4 函数极值的数值求解 【例4.1-7】（1）syms xy=sin(x)2*exp(-0.1*x)-0.5*sin(x)*(x+0.1);yd=diff(y,x); %xs0=solve(yd,x) % yd_xs0=vpa(subs(yd,x,xs0),6) %y_xs0=vpa(subs(y,x,xs0),6) %xs0 =0.050838341410271656880659496266968yd_xs0 =2.29589*10(
8、-41)y_xs0 =-0.00126332 （2）x1=-10;x2=10; %yx=(x)(sin(x)2*exp(-0.1*x)-0.5*sin(x)*(x+0.1); %xn0,fval,exitflag,output=fminbnd(yx,x1,x2) % % xn0 = 2.514797840754235fval = -0.499312445280039exitflag = 1output = iterations: 13 funcCount: 14 algorithm: golden section search, parabolic interpolation message:
9、 1x112 char （3）xx=-10:pi/200:10; %yxx=subs(y,x,xx);plot(xx,yxx)xlabel(x),grid on 图 4.1-5 在-10, 10区间中的函数曲线（4）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 =.023*fval = -3.568014059128578exitflag = 1output = iterations: 9 funcCount:
10、 10 algorithm: golden section search, parabolic interpolation message: 1x112 char 【例4.1-8】（1）ff=(x)(100*(x(2)-x(1).2)2+(1-x(1)2); （2）format short gx0=-5,-2,2,5;-5,-2,2,5; %提供4个搜索起点sx,sfval,sexit,soutput=fminsearch(ff,x0) %sx给出一组使优化函数值非减的局部极小点 sx = 0.99998 -0.68971 0.41507 8.0886 0.99997 -1.9168 4.96
11、43 7.8004sfval = 2.4112e-010sexit = 1soutput = iterations: 384 funcCount: 615 algorithm: Nelder-Mead simplex direct search message: 1x196 char （3）format short edisp(ff(sx(:,1),ff(sx(:,2),ff(sx(:,3),ff(sx(:,4) 2.4112e-010 5.7525e+002 2.2967e+003 3.3211e+005 .1.5 常微分方程的数值解 【例4.1-9】（1） （3）tspan=0,30; %
12、y0=1;0; %tt,yy=ode45(DyDt,tspan,y0); %plot(tt,yy(:,1)xlabel(t),title(x(t) 图 4.1-6 微分方程解（4）plot(yy(:,1),yy(:,2) %xlabel(位移),ylabel(速度) 图 4.1-7 平面相轨迹.2 矩阵和代数方程.2.1 矩阵运算和特征参数10 1 矩阵运算【例 4.2-1】（1）clearrng(default)rng(12)A=rand(2,4);B=rand(4,3);C1=zeros(size(A,1),size(B,2); % for ii=1:size(A,1) for jj=1:
13、size(B,2) for k=1:size(A,2) C1(ii,jj)=C1(ii,jj)+A(ii,k)*B(k,jj); end endend %C1 C1 = 0.7337 0.8396 0.3689 1.0624 1.1734 1.3787 （2）C2=zeros(size(A,1),size(B,2); for jj=1:size(B,2) for k=1:size(B,1) C2(:,jj)=C2(:,jj)+A(:,k)*B(k,jj); endend C2 C2 = 0.7337 0.8396 0.3689 1.0624 1.1734 1.3787 （3） C3=A*B,
14、% C3 = 0.7337 0.8396 0.3689 1.0624 1.1734 1.3787 （4 ）C3_C3=norm(C3-C1,fro) % C3_C2=norm(C3-C2,fro) % C3_C3 = 0C3_C2 = 0 【例 4.2-2】format rat %A=magic(2)+1j*pascal(2) %A = 1 + 1i 3 + 1i 4 + 1i 2 + 2i A1=A %A2=A. % A1 = 1 - 1i 4 - 1i 3 - 1i 2 - 2i A2 = 1 + 1i 4 + 1i 3 + 1i 2 + 2i B1=A*AB2=A.*AC1=A*A.C2
15、=A.*A. B1 = 12 13 - 1i 13 + 1i 25 B2 = 2 13 + 1i 13 - 1i 8 C1 = 8 + 8i 7 + 13i 7 + 13i 15 + 16i C2 = 0 + 2i 11 + 7i 11 + 7i 0 + 8i 10 2 矩阵的标量特征参数 【例4.2-3】A=reshape(1:9,3,3) %r=rank(A) %d3=det(A) %d2=det(A(1:2,1:2) %t=trace(A) % A = 1 4 7 2 5 8 3 6 9 r = 2 d3 = 0 d2 = -3 t = 15 【例4.2-4】format short %
16、rng default %A=rand(3,3); %B=rand(3,3); %C=rand(3,4);D=rand(4,3);tAB=trace(A*B) %tBA=trace(B*A) tCD=trace(C*D) %tDC=trace(D*C) tAB = 3.7479tBA = 3.7479tCD = 3.3399tDC = 3.3399 d_A_B=det(A)*det(B)dAB=det(A*B) dBA=det(B*A) d_A_B = -0.0852dAB = -0.0852dBA = -0.0852 dCD=det(C*D)dDC=det(D*C) % dCD = -0.0
17、557dDC = 1.5085e-017 .2.2 矩阵的变换和特征值分解【例4.2-5】（1）A=magic(4) %R,ci=rref(A) %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 （2）r_A=length(ci) r_A = 3 （3） aa=A(:,1:3)*R(1:3,4) %err=norm(A(:,4)-aa) % aa = 13 8 12 1err = 0 【例4.2-6】A=reshape(1:15,5,3); %X=null(A) %
18、S=A*X %n=size(A,2); %l=size(X,2); %n-l=rank(A) % X = -0.4082 0.8165 -0.4082S = 1.0e-014 * 0.2665 0.2665 0.3553 0.4441 0.6217ans = 1 【例4.2-7】 （1）A=1,-3;2,2/3V,D=eig(A) A = 1.0000 -3.0000 2.0000 0.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
19、=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) %A2=VR*DR/VRerr1=norm(A-A1,fro)err2=norm(A-A2,fro) A1 = 1.0000 -3.0000 - 0.0000i 2.0000 0.6667 A1_1 = 1.0000 -3.0000 2.0000 0.6667A2 = 1.0000 -3.0000 2.0000 0.6667err1 = 4.6613e-016err2 = 4.4409e-0
20、16 .2.3 线性方程的解10 1 线性方程解的一般结论 【例4.2-8】（1）A=reshape(1:12,4,3); %b=(13:16); % （2）ra=rank(A) %Arab=rank(A,b) % ra = 2rab = 2 （3）xs=Ab; %xg=null(A); %c=rand(1); %ba=A*(xs+c*xg) %norm(ba-b) % Warning: Rank deficient, rank = 2, tol = 1.8757e-014.ba = 13.0000 14.0000 15.0000 16.0000ans = 1.1374e-014 10 2 矩
21、阵逆【例4.2-9】（1）rng defaultA=gallery(randsvd,300,2e13,2); %x=ones(300,1); %b=A*x; %cond(A) % ans = 2.0022e+013 （2）tic %xi=inv(A)*b; % ti=toc %eri=norm(x-xi) %rei=norm(A*xi-b)/norm(b) % ti = 0.2034eri = 0.0812rei = 0.0047 （3）tic;xd=Ab; %td=tocerd=norm(x-xd)red=norm(A*xd-b)/norm(b) td = 0.0125erd = 0.027
22、4red = 7.5585e-015 .2.4 一般代数方程的解 【例 4.2-10】（1）syms tft=sin(t)2*exp(-0.1*t)-0.5*abs(t);S=solve(ft,t) % ftS=subs(ft,t,S) % S =0ftS =0 （2） （A）y_C=inline(sin(t).2.*exp(-0.1*t)-0.5*abs(t),t); % （B）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)hold off 图
23、4.2-1 函数零点分布观察图 （C） zoom on %tt,yy=ginput(5);zoom off % 图 4.2-2 局部放大和利用鼠标取值图tt % tt = -2.0039 -0.5184 -0.0042 0.6052 1.6717 （D）t4,y4=fzero(y_C,0.1) % t4 = 0.5993y4 = 1.1102e-016 .3 概率分布和统计分析10 1 二项分布(Binomial distribution) 【例 4.3-1】N=100;p=0.5; %k=0:N; %pdf=binopdf(k,N,p); %cdf=binocdf(k,N,p); %h=pl
24、otyy(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 二项分布B(100, 0.5)的概率和累计概率曲线10 2 正态分布（Normal distribu
25、tion） 【例4.3-2】mu=3;sigma=0.5; %x=mu+sigma*-3:-1,1:3;yf=normcdf(x,mu,sigma);P=yf(4)-yf(3),yf(5)-yf(2),yf(6)-yf(1); %xd=1:0.1:5;yd=normpdf(xd,mu,sigma); %clffor k=1:3 xx=x(4-k):sigma/10:x(3+k); yy=normpdf(xx,mu,sigma); subplot(3,1,k),plot(xd,yd,b); % hold on fill(x(4-k),xx,x(3+k),0,yy,0,g); % hold off
26、 if k2 text(3.8,0.6,mu-sigma,mu+sigma) else kk=int2str(k); text(3.8,0.6,mu-,kk,sigma,mu+,kk,sigma) end text(2.8,0.3,num2str(P(k);shgend xlabel(x);shg 图 4.3-2 均值两侧一、二、三倍标准差之间的概率10 3 各种概率分布的交互式观察界面 【例4.3-3】图4.3-3 概率分布交互界面.3.2 全局随机流、随机数组和统计分析10 1 全局随机流的操控表4.3-1 产生全局随机流的发生器generator算 例可取字符串含义twister默认的随机数发生器例v5u
