北航Matlab教程(R2011a)习题4解答

习题4

1．根据题给的模拟实际测量数据的一组和试用数值差分diff或数值梯度gradient指令计算，然后把和曲线绘制在同一张图上，观察数值求导的后果。（模拟数据从prob_data401.mat获得）（提示：自变量采样间距太小。）

load prob_401;

N=20;

diff_y1=(diff(y(1:N:end)))./diff(t(1:N:end));

gradient_y1=(gradient(y(1:N:end)))./gradient(t(1:N:end));

t1=t(1:N:end);

length(t1)

plot(t,y,t1(1:end-1),diff_y1)

plot(t,y,t1,gradient_y1)

2．采用数值计算方法，画出在区间曲线，并计算。（提示：cumtrapz快捷，在精度要求不高处可用；quad也可试。巧用find。）

d=0.5;

tt=0:d:10;

t=tt+(tt==0)*eps;

y=sin(t)./t;

s=d*trapz(y)

ss=d*(cumtrapz(y))

plot(t,y,t,ss,'r'),hold on

y4_5=ss(find(t==4.5))

yi=interp1(t,ss,4.5),plot(4.5,yi,'r+')

3．求函数的数值积分，并请采用符号计算尝试复算。（提示：各种数值法均可试。）

d=pi/20;

x=0:d:pi;

fx=exp(sin(x).^3);

s=d*trapz(fx)

s1=quad('exp(sin(x).^3)',0,pi)

s2=quadl('exp(sin(x).^3)',0,pi)

s3=vpa(int('exp(sin(x)^3)',0,pi))

s4=vpa(int(sym('exp(sin(x)^3)'),0,pi))

4．用quad求取的数值积分，并保证积分的绝对精度为。（体验：试用trapz，如何算得同样精度的积分。）

s1=quad('exp(-abs(x)).*abs(sin(x))',-5*pi,1.7*pi,1e-10)

s2=quadl('exp(-abs(x)).*abs(sin(x))',-5*pi,1.7*pi)

syms x;

s3=vpa(int(exp(-abs(x))*abs(sin(x)),-5*pi,1.7*pi))

d=pi/1000;

x=-5*pi:d:1.7*pi;

fx=exp(-abs(x)).*abs(sin(x));

s=d*trapz(fx)

5．求函数在区间中的最小值点。（提示：作图观察。）

x1=-5;

x2=5;

yx=inline('(sin(5*t)).^2.*exp(0.06*t.^2)-

1.5.*t.*cos(2*t)+1.8.*abs(t+0.5)')

[xn0,fval]=fminbnd(yx,x1,x2)

t=x1:0.1:x2;

plot(t,yx(t)),hold on ,plot(xn0,fval,'r*')

6．设

，用数值法和符号法求。（提示：注意ode45和 dsolve的用法。）tspan=[0,0.5];

y0=[1;0];

[tt,yy]=ode45(@DyDt_6,tspan,y0);

y0_5=yy(end,1)

S = dsolve('D2y-3*Dy+2*y = 1','y(0) = 1','Dy(0) = 0')

ys0_5=subs(S,0.5)

function ydot=DyDt_6(t,y)

mu=3;

ydot=[y(2);mu*y(2)-2*y(1)+1];

7．已知矩阵A=magic(8)，（1）求该矩阵的“值空间基阵”B ；（2）写出“A的任何列可用基向量线性表出”的验证程序。（提示：方法很多；建议使用rref体验。）

A=magic(8)

B=orth(A)

rref(A)

rref(B)

8．已知由MATLAB指令创建的矩阵A=gallery(5)，试对该矩阵进行特征值分解，并通过验算观察发生的现象。（提示：condeig）

A=gallery(5)

[V,D,s]=condeig(A)

[V,D]=eig(A)

cond(A)

jordan(A)

9．求矩阵的解，A为3阶魔方阵，b是的全1列向量。（提示：用rref, inv, / 体验。）

A=magic(3)

b=ones(3,1)

x=A\b

x=inv(A)*b

rref([A,b])

10．求矩阵的解，A为4阶魔方阵，b是的全1列向量。（提示：用rref, inv, / 体验。）

A=magic(4)

b=ones(4,1)

x=A\b

xg=null(A)

11．求矩阵的解，A为4阶魔方阵，。（提示：用rref, inv, / 体验。）

A=magic(4)

b=(1:4)'

rref([A,b])

x=A\b

A*x

x=inv(A)*b

12．求的实数解。（提示：发挥作图法功用）

y_C=inline('-0.5+t-10.*exp(-0.2.*t).*abs(sin(sin(t)))','t');

t=-10:0.01:10;

Y=y_C(t);

plot(t,Y,'r'),hold on

plot(t,zeros(size(t)),'k');

xlabel('t');ylabel('y(t)')

zoom on

[tt,yy]=ginput(1),zoom off

[t1,y1]=fzero(y_C,tt)

[t2,y2]=fsolve(y_C,tt)

13．求解二元函数方程组的解。（提示：可尝试符号法解；试用contour作图求解；比较之。此题有无数解。）

S=solve('sin(x-y)=0','cos(x+y)=0','x','y')

S.x, S.y

14．假定某窑工艺瓷器的烧制成品合格率为0.157，现该窑烧制100件瓷器，请画出合格产品数的概率分布曲线。（提示：二项式分布概率指令binopdf；stem）

y = binopdf([0:100],100,0.157);

stem(1:length(y),y)

axis([0 length(y) 0 .12 ])

15．试产生均值为4，标准差为2的的正态分布随机数组 a ，分别用hist和histfit绘制该数组的频数直方图，观察两张图形的差异。除histfit上的拟合红线外，你能使这两个指令绘出相同的频数直方图吗？（提示：为保证结果的重现性，在随机数组a产生前，先运行 rng default指令；可使用指令normrnd产生正态分布随机数；理解hist(Y, m)指令格式。）

a=normrnd(4,2,10000,1);

hist(a)

histfit(a)

hist(a,sqrt(10000))

16．从数据文件prob_data416.mat得到随机数组R，下面有一段求取随机数组全部数据最大值、均值和标准差的程序。

Mx=max(max(R)),Me=mean(mean(R)),St=std(std(R)),

试问该程序所得的结果都正确吗？假如不正确，请写出正确的程序。（提示：load；R(:)。）

load prob_416;

Mx=max(max(R))

Me=mean(mean(R))

St=std(R(:))

17．已知有理分式，其中，

。（1）求该分式的商多项式和余多项式。（2）用程序验算是否成立。（提示：采用范数指令norm验算。）

format rat

NX=conv([3,0,1,0],[1,0,0,0.5]), DX=conv([1,2,-2],[5,2,0,1])

[q,r]=deconv(NX,DX)

cq='商多项式为 ';cr='余多项式为 ';

disp([cq,poly2str(q,'s')]),disp([cr,poly2str(r,'s')])

qp2=conv(q,DX), pp1=qp2+r, pp1==NX

18．现有一组实验数据x, y（数据从prob_data418.mat获得），试求这组数据的5阶拟合多项式。（提示：load, polyfit, polyval）

load prob_418, who,x

P=polyfit(x,y,5), Pt=poly2str(P,'t')

xx=-1:0.01:4, yy=polyval(P,xx), plot(xx,yy,x,y,'*r')

legend('拟合曲线','原始曲线','Location','SouthEast')

19．已知系统冲激响应为h(n)=[0.05,0.24,0.40,0.24,0.15,-0.1,0.1] ，系统输入u(n)由指令rng default;u=2*(randn(1,100)>0.5)-1产生，该输入信号的起始作用时刻为0。试画出类似图p4-1所示的系统输入、输出信号图形。（提示：注意输入信号尾部的处理；NaN的使用。）

图 p4-1

h=[0.05,0.24,0.40,0.24,0.15,-0.1,0.1]

randn('state',1);

u=2*(randn(1,100)>0.5)-1;

y=conv(u,h);

subplot(2,1,1),stem(u,'filled') axis([0 length(y) -1 1 ]) subplot(2,1,2),stem(y,'filled') axis([0 length(y) -1 1 ])