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

习题4

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

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． 采用数值计算方法，画出dt t

t x y x

⎰=0sin )(在]10 ,0[区间曲线，并计算)5.4(y 。（提示：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． 求函数x e x f 3sin )(=的数值积分⎰=π

0 )(dx x f s ，并请采用符号计算尝试复算。（提

示：各种数值法均可试。）

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 求取dx x e x sin 7.15⎰--π

π的数值积分，并保证积分的绝对精度为910-。（体

验：试用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． 求函数5.08.12cos 5.1)5(sin )(2

06.02++-=t t t e t t f t 在区间]5,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． 设0)0(,1)0(,1)(2)(3)(22===+-dt

dy y t y dt t dy dt t y d ，用数值法和符号法求5.0)(=t t y 。（提示：注意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． 求矩阵b Ax =的解，A 为3阶魔方阵，b 是)13(⨯的全1列向量。（提示：用rref, inv, / 体验。）

A=magic(3)

b=ones(3,1)

x=A\b

x=inv(A)*b

rref([A,b])

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

A=magic(4)

b=ones(4,1)

x=A\b

xg=null(A)

11． 求矩阵b Ax =的解，A 为4阶魔方阵，⎥⎥⎥⎥

⎦

⎤⎢⎢⎢⎢⎣⎡=4321b 。（提示：用rref, inv, / 体验。） A=magic(4)

b=(1:4)'

rref([A,b])

x=A\b

A*x

x=inv(A)*b