ch4 - 数值计算mbook2008a

第 4 章 数值计算

与符号计算相比，数值计算在科研和工程中的应用更为广泛。MATLAB也正是凭借其卓越的数值计算能力而称雄世界。随着科研领域、工程实践的数字化进程的深入，具有数字化本质的数值计算就显得愈益重要。 较之十年、二十年前，在计算机软硬件的支持下当今人们所能拥有的计算能力已经得到了巨大的提升。这就自然地激发了人们从新的计算能力出发去学习、理解概念的欲望，鼓舞了人们用新计算能力试探解决实际问题的雄心。 鉴于当今高校本科教学偏重符号计算和便于手算简单示例的实际，也出于帮助读者克服对数值计算生疏感的考虑，本章在内容安排上仍从“微积分”开始。这一方面与第2章符号计算相呼应，另方面通过“微积分”说明数值计算离散本质的微观和宏观影响。 为便于读者学习，本章内容的展开脉络基本上沿循高校数学教程，而内容深度力求控制在高校本科水平。考虑到知识的跳跃和交叉，本书对重要概念、算式、指令进行了尽可能完整地说明。

4.1

4.1.1

数值微积分

近似数值极限及导数

1?cos2xsinx，f2(x)?，试用机器零阈值eps替代理论0计算

xsinxx极限L1(0)?limf1(x)，L2(0)?limf2(x)。

【例4.1-1】设f1(x)?x?0x?0

x=eps;

L1=(1-cos(2*x))/(x*sin(x)), L2=sin(x)/x, L1 = 0 L2 =

1

syms t

f1=(1-cos(2*t))/(t*sin(t)); f2=sin(t)/t;

Ls1=limit(f1,t,0) Ls2=limit(f2,t,0) Ls1 = 2

Ls2 = 1

d=pi/100; t=0:d:2*pi; x=sin(t);

dt=5*eps; x_eps=sin(t+dt);

【例4.1-2】已知x?sin(t)，求该函数在区间 [0, 2?]中的近似导函数。

1

dxdt_eps=(x_eps-x)/dt; plot(t,x,'LineWidth',5) hold on

plot(t,dxdt_eps) hold off

legend('x(t)','dx/dt') xlabel('t')

图 4.1-1 增量过小引起有效数字严重丢失后的毛刺曲线

x_d=sin(t+d);

dxdt_d=(x_d-x)/d; plot(t,x,'LineWidth',5) hold on

plot(t,dxdt_d) hold off

legend('x(t)','dx/dt') xlabel('t')

图 4.1-2 增量适当所得导函数比较光滑

【例4.1-3】已知x?sin(t)，采用diff和gradient计算该函数在区间 [0, 2?]中的近似导函数。。

clf

d=pi/100; t=0:d:2*pi; x=sin(t);

dxdt_diff=diff(x)/d; dxdt_grad=gradient(x)/d; subplot(1,2,1)

2

plot(t,x,'b') hold on

plot(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, 2\\pi]')

legend('x(t)','dxdt_{grad}','dxdt_{diff}','Location','North') xlabel('t'),box off hold off

subplot(1,2,2)

kk=(length(t)-10):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

图 4.1-3 diff和gradient求数值近似导数的异同比较

4.1.2 数值求和与近似数值积分

【例 4.1-4】求积分s(x)???/20y(t)dt，其中y?0.2?sin(t)。

clear

d=pi/8; t=0:d:pi/2; y=0.2+sin(t); s=sum(y); s_sa=d*s; s_ta=d*trapz(y);

disp(['sum求得积分',blanks(3),'trapz求得积分']) disp([s_sa, s_ta]) t2=[t,t(end)+d]; y2=[y,nan]; stairs(t2,y2,':k') hold on

plot(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

3

shg

sum求得积分 trapz求得积分

1.5762 1.3013

图 4.1-4 sum 和trapz求积模式示意

4.1.3 计算精度可控的数值积分

【例 4.1-5】求I??e01?x2dx 。

syms x

Isym=vpa(int(exp(-x^2),x,0,1)) Isym =

.74682413281242702539946743613185

format long d=0.001;x=0:d:1;

Itrapz=d*trapz(exp(-x.*x)) Itrapz =

0.746824071499185

fx='exp(-x.^2)'; Ic=quad(fx,0,1,1e-8) Ic =

0.746824132854452

【例 4.1-6】求s??? 1 2 1 0xydxdy。

syms x y

s=vpa(int(int(x^y,x,0,1),y,1,2)) s =

.40546510810816438197801311546432

format long

s_n=dblquad(@(x,y)x.^y,0,1,1,2) s_n =

0.405466267243508

4.1.4 函数极值的数值求解

【例4.1-7】已知

y?(x??)?esin(x??)，在??/2?x??/2区间，求函数的极小值。

4

（1）

syms x

y=(x+pi)*exp(abs(sin(x+pi))); yd=diff(y,x); xs0=solve(yd) yd_xs0=vpa(subs(yd,x,xs0),6) y_xs0=vpa(subs(y,x,xs0),6) y_m_pi=vpa(subs(y,x,-pi/2),6) y_p_pi=vpa(subs(y,x,pi/2),6)

Warning: Warning, solutions may have been lost xs0 =

-1.0676598444985783372948670854801 yd_xs0 = 0.

y_xs0 = 4.98043 y_m_pi = 4.26987 y_p_pi = 12.8096

（2）

x1=-pi/2;x2=pi/2; yx=@(x)(x+pi)*exp(abs(sin(x+pi)));

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

-1.2999e-005 fval =

3.1416 exitflag = 1 output =

iterations: 21 funcCount: 22

algorithm: 'golden section search, parabolic interpolation' message: [1x112 char]

（3）

xx=-pi/2:pi/200:pi/2; yxx=(xx+pi).*exp(abs(sin(xx+pi))); plot(xx,yxx)

xlabel('x'),grid on 131211109876543-2-1.5-1-0.50x0.511.52图 4.1-5 在[-pi/2,pi/2]区间中的函数曲线 5

[xx,yy]=ginput(1) xx= 1.5054e-008 yy= 3.1416

图 4.1-6 函数极值点附近的局部放大和交互式取值

【例4.1-8】求

f(x,y)?100(y?x2)2?(1?x)2的极小值点。它即是著名的Rosenbrock's

\测试函数，它的理论极小值是x（1）

（2）

?1,y?1。

ff=@(x)(100*(x(2)-x(1).^2)^2+(1-x(1))^2);

x0=[-5,-2,2,5;-5,-2,2,5];

[sx,sfval,sexit,soutput]=fminsearch(ff,x0) sx =

1.0000 -0.6897 0.4151 8.0886 1.0000 -1.9168 4.9643 7.8004 sfval =

2.4112e-010 sexit = 1

soutput =

iterations: 384 funcCount: 615

algorithm: [1x33 char] message: [1x196 char]

（3）检查目标函数值

format short e

disp([ff(sx(:,1)),ff(sx(:,2)),ff(sx(:,3)),ff(sx(:,4))]) Columns 1 through 3

2.4112e-010 5.7525e+002 2.2967e+003 Column 4

3.3211e+005

4.1.5

常微分方程的数值解

6

d2x2dx??(1?x)?x?0，??2，在初始条件【例 4.1-9】求微分方程2dtdtdx(0)x(0)?1,?0情况下的解，并图示。（见图4.1-7和4.1-8）

dt（1） （2） [DyDt.m]

function ydot=DyDt(t,y) mu=2;

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

（3）解算微分方程

tspan=[0,30]; y0=[1;0];

[tt,yy]=ode45(@DyDt,tspan,y0); plot(tt,yy(:,1))

xlabel('t'),title('x(t)')

图 4.1-7 微分方程解

（4）

plot(yy(:,1),yy(:,2))

xlabel('位移'),ylabel('速度')

图 4.1-8 平面相轨迹

7

4.2

4.2.1 一

矩阵和代数方程

矩阵运算和特征参数 矩阵运算

【例 4.2-1】已知矩阵A2?4，B4?3，采用三种不同的编程求这两个矩阵的乘积

C2?3?A2?4B4?3。

（1）

clear

rand('twister',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:size(B,2) for k=1:size(A,2) C1(ii,jj)=C1(ii,jj)+A(ii,k)*B(k,jj); end end end 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); end end C2 C2 =

0.7337 0.8396 0.3689 1.0624 1.1734 1.3787

（3）

C3=A*B, C3 =

0.7337 0.8396 0.3689 1.0624 1.1734 1.3787

（3）

C3_C3=norm(C3-C1,'fro') C3_C2=norm(C3-C2,'fro') C3_C3 = 0 C3_C2 = 0

【例 4.2-2】观察矩阵的转置操作和数组转置操作的差别。

format rat A=magic(2)+j*pascal(2) A =

Column 1

8

1 + 1i 4 + 1i Column 2

3 + 1i 2 + 2i %

A1=A' A2=A.' A1 =

Column 1

1 - 1i 3 - 1i Column 2

4 - 1i 2 - 2i A2 =

Column 1

1 + 1i 3 + 1i Column 2

4 + 1i 2 + 2i

B1=A*A' B2=A.*A' C1=A*A.' C2=A.*A.' B1 =

Column 1

12 13 + 1i Column 2

13 - 1i 25 B2 =

Column 1

2 13 - 1i Column 2

13 + 1i 8 C1 =

Column 1

8 + 8i 7 + 13i Column 2

7 + 13i 15 + 16i C2 =

Column 1

0 + 2i 11 + 7i Column 2

11 + 7i 0 + 8i

二 矩阵的标量特征参数

【例4.2-3】矩阵标量特征参数计算示例。A=reshape(1:9,3,3) r=rank(A)

9

d3=det(A) d2=det(A(1:2,1:2)) t=trace(A) A =

Columns 1 through 2

1 4 2 5 3 6 Column 3

7 8 9 r =

2 d3 =

0 d2 =

-3 t =

15

【例4.2-4】矩阵标量特征参数的性质。

format short

rand('twister',0) 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 =

2.6030 tBA =

2.6030 tCD =

4.1191 tDC =

4.1191

d_A_B=det(A)*det(B) dAB=det(A*B) dBA=det(B*A) d_A_B = 0.0094 dAB =

0.0094 dBA =

0.0094

dCD=det(C*D)

dDC=det(D*C) dCD =

0.0424 dDC =

-2.6800e-018

10

4.2.2 矩阵的变换和特征值分解

【例4.2-5】行阶梯阵简化指令rref计算结果的含义。 （1）

A=magic(4) [R,ci]=rref(A) A =

16 2 3 13 5 11 10 8 9 7 6 12 4 14 15 1 R =

1 0 0 1 0 1 0 3 0 0 1 -3 0 0 0 0 ci =

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 1 err = 0

【例4.2-6】矩阵零空间及其含义。

A=reshape(1:15,5,3); X=null(A) S=A*X n=size(A,2); l=size(X,2); n-l==rank(A)

X =

0.4082 -0.8165 0.4082 S =

1.0e-014 * -0.2665 -0.1776 -0.0888 -0.0888 -0.0888 ans = 1

【例4.2-7】简单实阵的特征值分解。（1）

A=[1,-3;2,2/3] [V,D]=eig(A) A =

1.0000 -3.0000

11

2.0000 0.6667 V =

0.7746 0.7746 0.0430 - 0.6310i 0.0430 + 0.6310i D =

0.8333 + 2.4438i 0 0 0.8333 - 2.4438i

（2）

[VR,DR]=cdf2rdf(V,D) VR =

0.7746 0 0.0430 -0.6310 DR =

0.8333 2.4438 -2.4438 0.8333

（3）

A1=V*D/V A1_1=real(A1) A2=VR*DR/VR

err1=norm(A-A1,'fro') err2=norm(A-A2,'fro') A1 =

1.0000 + 0.0000i -3.0000 2.0000 - 0.0000i 0.6667 A1_1 =

1.0000 -3.0000 2.0000 0.6667 A2 =

1.0000 -3.0000 2.0000 0.6667 err1 =

7.0290e-016 err2 =

4.4409e-016

4.2.3 一 二

线性方程的解 线性方程解的一般结论 除法运算解方程

59??13??14?610???x???的解。 711??15????812??16?

?1?2【例4.2-8】求方程??3??4（1）

A=reshape(1:12,4,3); b=(13:16)';

（2）

ra=rank(A) rab=rank([A,b]) ra = 2 rab = 2

12

（3）

xs=A\\b; 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.0000 ans =

1.0658e-014

三 矩阵逆

【例4.2-9】“逆阵”法和“左除”法解恰定方程的性能对比 （1）

randn('state',0);

A=gallery('randsvd',300,2e13,2); x=ones(300,1); b=A*x; cond(A) ans =

1.9997e+013

（2）

tic xi=inv(A)*b; ti=toc eri=norm(x-xi) rei=norm(A*xi-b)/norm(b) ti =

0.0728 eri =

0.0918 rei =

0.0050

（3）

tic;

xd=A\\b; td=toc

erd=norm(x-xd)

red=norm(A*xd-b)/norm(b) td =

0.0172 erd =

0.0291 red =

9.7167e-015

4.2.4 一般代数方程的解

2?0.1t【例 4.2-10】求f(t)?(sint)?e?0.5 t 的零点。

13

（1）

S=solve('sin(t)^2*exp(-0.1*t)-0.5*abs(t)','t') S = 0.

（2）

y_C=inline('sin(t).^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 on

plot(t,zeros(size(t)),'k'); xlabel('t');ylabel('y(t)') hold off 1

0-1y(t)-2-3-4-5-10-8-6-4-20t246810图 4.2-1 函数零点分布观察图

zoom on [tt,yy]=ginput(5);zoom off

图 4.2-2 局部放大和利用鼠标取值图

tt tt =

-2.0039 -0.5184 -0.0042 0.6052 1.6717

[t4,y4]=fzero(y_C,0.1) t4 =

0.5993

14

y4 =

1.1102e-016

4.3

4.3.1 一

概率分布和统计分析

概率函数、分布函数、逆分布函数和随机数的发生 二项分布(Binomial distribution)

【例 4.3-1】画出N=100, p=0.5情况下的二项分布概率特性曲线。

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 二项分布B(100, 0.5)的概率和累计概率曲线

二 正态分布（Normal distribution）

【例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); clf

for k=1:3 %------------------------------- xx=x(4-k):sigma/10:x(3+k);

15

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 if k<2

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)));shg end

xlabel('x');shg

图 4.3-2 均值两侧一、二、三倍标准差之间的概率

三 各种概率分布的交互式观察界面

【例4.3-3】概率分布交互界面使用方法介绍。

16

图4.3-3 概率分布交互界面

4.3.2 随机数发生器和 统计分析指令

【例4.3-4】随机数据的统计量。

randn('state',0) A=randn(1000,4); AMAX=max(A) AMIN=min(A) CM=mean(A) MA=mean(mean(A)) S=std(A) var(A)-S.^2 C=cov(A)

diag(C)'-var(A) p=corrcoef(A) AMAX =

2.7316 3.2025 3.4128 3.0868 AMIN =

-2.6442 -3.0737 -3.5027 -3.0461 CM =

-0.0431 0.0455 0.0177 0.0263 MA =

0.0116 S =

0.9435 1.0313 1.0248 0.9913 ans =

1.0e-015 *

0 -0.2220 0 0 C =

0.8902 -0.0528 0.0462 0.0078 -0.0528 1.0635 0.0025 0.0408 0.0462 0.0025 1.0502 -0.0150 0.0078 0.0408 -0.0150 0.9826 ans =

1.0e-014 *

-0.0111 -0.1554 -0.0888 0

17

p =

1.0000 -0.0543 0.0478 0.0083 -0.0543 1.0000 0.0024 0.0399 0.0478 0.0024 1.0000 -0.0147 0.0083 0.0399 -0.0147 1.0000

【例4.3-5】产生1000 个服从N(2,0.52)的随机数。

mu=2;s=0.5;

randn('state',22) x=randn(1000,1); y=s*x+mu; z=s*(x+mu);

subplot(3,1,1),histfit(x),axis([-5,5,0,100]),ylabel('x') subplot(3,1,2),histfit(y),axis([-5,5,0,100]),ylabel('y')

subplot(3,1,3),histfit(z),axis([-5,5,0,100]),ylabel('z')

图 4.3-4 均值为2标准差为0.5的随机数样本 z

4.4

4.4.1 一 二

多项式运算和卷积

多项式的运算函数 多项式表达方式的约定 多项式运算函数

(s2?2)(s?4)(s?1)【例4.4-1】求的“商”及“余”多项式。

s3?s?1（1）

format rat

p1=conv([1,0,2],conv([1,4],[1,1])); p2=[1 0 1 1]; [q,r]=deconv(p1,p2);

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

18

disp([cq,poly2str(q,'s')]),disp([cr,poly2str(r,'s')]) 商多项式为 s + 5

余多项式为 5 s^2 + 4 s + 3

（2）

qp2=conv(q,p2); pp1=qp2+r;

pp1==p1 ans =

1 1 1 1 1

【例 4.4-2】矩阵和特征多项式，特征值和多项式根。 （1）

A=[11 12 13;14 15 16;17 18 19]; PA=poly(A) PPA=poly2str(PA,'s') PA =

1.0000 -45.0000 -18.0000 0.0000 PPA =

s^3 - 45 s^2 - 18 s + 1.6206e-014

（2）

s=eig(A)

r=roots(PA) s =

45.3965 -0.3965 0.0000 r =

45.3965 -0.3965 0.0000

（3）

n = length(PA); AA = diag(ones(1,n-2,class(PA)),-1); AA(1,:) = -PA(2:n)/ PA(1); AA

sr = eig(AA) AA =

45.0000 18.0000 -0.0000 1.0000 0 0 0 1.0000 0 sr =

45.3965 -0.3965 0.0000

【例 4.4-3】构造指定特征根的多项式。

R=[-0.5,-0.3+0.4*i,-0.3-0.4*i]; P=poly(R) PR=real(P) PPR=poly2str(PR,'x') P =

1.0000 1.1000 0.5500 0.1250 PR =

1.0000 1.1000 0.5500 0.1250 PPR =

x^3 + 1.1 x^2 + 0.55 x + 0.125

19

【例 4.4-4】多项式求值指令polyval与polyvalm的本质差别。 （1）

clear

p=[1,2,3]; poly2str(p,'x') X=[1,2;3,4] ans =

x^2 + 2 x + 3 X =

1 2 3 4

（2）

va=X.^2+2*X+3 Va=polyval(p,X) va =

6 11 18 27 Va =

6 11 18 27

（3）

vm=X^2+2*X+3*eye(2) Vm=polyvalm(p,X) vm =

12 14 21 33 Vm =

12 14

21 33

（4）

cp=poly(X); poly2str(cp,'x')

cpXa=polyval(cp,X) cpX=polyvalm(cp,X) ans =

x^2 - 5 x - 2 cpXa =

-6 -8 -8 -6 cpX =

1.0e-015 *

0.2220 0 0 0.2220

4.4.2 一

多项式拟合和最小二乘法 多项式拟合

【例 4.4-5】 给定数据组x0 , y0 ，求拟合三阶多项式，并图示拟合情况。（见图4.4-1）

x0=0:0.1:1;

y0=[-.447,1.978,3.11,5.25,5.02,4.66,4.01,4.58,3.45,5.35,9.22];

n=3;

P=polyfit(x0,y0,n)

20

xx=0:0.01:1;

yy=polyval(P,xx);

plot(xx,yy,'-b',x0,y0,'.r','MarkerSize',20)

legend('拟合曲线','原始数据','Location','SouthEast') xlabel('x') P =

56.6915 -87.1174 40.0070 -0.9043

图 4.4-1 采用三阶多项式所得的拟合曲线

二 最小二乘问题

图 4.4-2 最小二乘的几何解释

【例 4.4-6】采用与例4.4-5相同的数据组x0 , y0 ，运用式(4.4-3)求拟合多项式的系数。

x0=(0:0.1:1)';

y0=[-.447,1.978,3.11,5.25,5.02,4.66,4.01,4.58,3.45,5.35,9.22]'; m=length(x0);

n=3;

X=zeros(m,n+1); for k=1:n

X(:,n-k+1)=(x0.^k);

21

end

X(:,n+1)=ones(m,1);

aT=(X\\y0)' aT =

56.6915 -87.1174 40.0070 -0.9043

4.4.3 两个有限长序列的卷积

?1n?3,4,?,12?1 和 B(n)??n?2,3,?,9，求这

【例4.4-7】有序列A(n)???0两个序列的卷积。 （1）

N1=3;N2=12;

A=ones(1,(N2-N1+1)); M1=2;M2=9;

B=ones(1,(M2-M1+1)); Nc1=N1+M1;Nc2=N2+M2; kcc=Nc1:Nc2; for n=Nc1:Nc2 w=0; for k=N1:N2 kk=k-N1+1; t=n-k; if t>=M1&t<=M2 tt=t-M1+1; w=w+A(kk)*B(tt); end end nn=n-Nc1+1; cc(nn)=w; end

kcc,cc kcc =

Columns 1 through 12

5 6 7 8 9 10 11 12 13 14 15 16 Columns 13 through 17

17 18 19 20 21 cc =

Columns 1 through 12

1 2 3 4 5 6 7 8 8 8 7 6 Columns 13 through 17

5 4 3 2 1

（2）

N1=3;N2=12;

a=ones(1,N2+1);a(1:N1)=0; M1=2;M2=9;

b=ones(1,M2+1);b(1:M1)=0; c=conv(a,b); kc=0:(N2+M2); kc,c kc =

Columns 1 through 12

0 1 2 3 4 5 6 7 8 9 10 11 Columns 13 through 22

12 13 14 15 16 17 18 19 20 21 c =

Columns 1 through 12

0 0 0 0 0 1 2 3 4 5 6 7

else 22

?0else Columns 13 through 22

8 8 8 7 6 5 4 3 2 1

（3）

N1=3;N2=12; M1=2;M2=9;

A=ones(1,(N2-N1+1)); B=ones(1,(M2-M1+1)); C=conv(A,B); Nc1=N1+M1;Nc2=N2+M2; KC=Nc1:Nc2; KC,C KC =

5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 C =

1 2 3 4 5 6 7 8 8 8 7 6 5 4 3 2 1

（3）

subplot(2,1,1),stem(kc,c), text(20,6,'0 起点法') CC=[zeros(1,KC(1)),C];

subplot(2,1,2),stem(kc,CC),text(18,6,'非平凡区间法') xlabel('n')

图 4.4-3 借助conv指令时两种不同序列记述法所得的卷积序列

习题4

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

〖答案〗

23

1.510.50-0.5-1-1.5-201234567

2． 采用数值计算方法，画出y(x)?

〖答案〗 s45 = 1.6541 21.81.61.41.210.80.60.40.200246810 ?x0sintdt在[0, 10]区间曲线，并计算y(4.5)。 t

3． 求函数

〖答案〗 s =

5.1370

Warning: Explicit integral could not be found.

f(x)?esin3x的数值积分s?? ? 0f(x)dx，并请采用符号计算尝试复算。

24

> In sym.int at 58 ss =

int(exp(sin(x)^3),x = 0 .. pi)

4． 用quad求取

〖答案〗 sq =

1.08784993815498

20.06t?1.5tcos2t?1.8t?0.5在区间[?5,5]中的最小值点。 5． 求函数f(t)?(sin5t)e2?1.7??5?e?xsinxdx的数值积分，并保证积分的绝对精度为10?9。

〖答案〗

最小值点是 相应目标值是

-1.28498111480531 -0.18604801006545

d2y(t)dy(t)dy(0)?3?2y(t)?1,y(0)?1,?0，用数值法和符号法求y(t)t?0.5。6． 设 2dtdtdt

〖答案〗 数值解

y_05 =

0.78958020790127

符号解

ys =

1/2-1/2*exp(2*t)+exp(t) ys_05 =

.78958035647060552916850705213780

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

〖答案〗

三组不同的基 B1 =

64 2 3 9 55 54 17 47 46 40 26 27 32 34 35 41 23 22 49 15 14 8 58 59

B2 =

-0.3536 0.5401 0.3536 -0.3536 -0.3858 -0.3536 -0.3536 -0.2315 -0.3536

25

-0.3536 0.0772 0.3536 -0.3536 -0.0772 0.3536 -0.3536 0.2315 -0.3536 -0.3536 0.3858 -0.3536 -0.3536 -0.5401 0.3536

B3 =

0.3536 0.6270 0.3913 0.3536 -0.4815 -0.2458 0.3536 -0.3361 -0.1004 0.3536 0.1906 -0.0451 0.3536 0.0451 -0.1906 0.3536 0.1004 0.3361 0.3536 0.2458 0.4815 0.3536 -0.3913 -0.6270

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

〖答案〗

用eig得到不正确的特征值 ans =

-0.0181 -0.0054 - 0.0171i -0.0054 + 0.0171i 0.0144 - 0.0104i 0.0144 + 0.0104i

应采用jordan才能得到正确的特征值

DJ =

0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0

9． 求矩阵Ax?b的解，A为3阶魔方阵，b是(3?1)的全1列向量。

〖答案〗 x =

0.0667 0.0667 0.0667

10． 求矩阵Ax?b的解，A为4阶魔方阵，b是(4?1)的全1列向量。

〖答案〗 解不唯一 x =

-0.0074 -0.0809 0.1397 0.0662

26

0.0588 0.1176 -0.0588 0

下面逆阵法的解不正确 xx =

0.0313 -0.1250 0.1250 0.0313

?1??2?11． 求矩阵Ax?b的解，A为4阶魔方阵，b???。

?3????4?〖答案〗

没有准确解，以下是伪逆阵求得的最小二乘解 x =

0.0235 0.1235 0.1235 0.0235

?0.2t12． 求?0.5?t?10esin[sint]?0的实数解。

〖答案〗

t =

2.7341 yt =

2.2204e-015 420-2-4-6-8-10-12-1012345

27

13． 求解二元函数方程组?

〖答案〗

X =

[ 1/4*pi] [ -1/4*pi] Y =

[ 1/4*pi] [ -1/4*pi]

?sin(x?y)?0的解。（用符号法解）

?cos(x?y)?0

14． 假定某窑工艺瓷器的烧制成品合格率为0.157，现该窑烧制100件瓷器，请画出合格产品数的概率分布曲线。

〖答案〗 0.120.10.080.060.040.020020406080100

?1)的正态分布随机数组 a ，15． 试产生均值为4，标准差为2的(10000 分别用hist和histfit

绘制该数组的频数直方图，观察两张图形的差异。除histfit上的拟合红线外，你能使这两个指令绘出相同的频数直方图吗？

〖答案〗

28

3000200010000-44003002001000-4-2024681012-2024681012

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

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

试问该程序所得的结果都正确吗？假如不正确，请写出正确的程序。

〖答案〗 Mx1 =

0.9997 Me1 =

0.5052 S1 =

0.2919

17． 已知有理分式R(x)?N(x)33，其中N(x)?(3x?x)(x?0.5)，

D(x)D(x)?(x2?2x?2)(5x3?2x2?1)。（1）求该分式的商多项式Q(x)和余多项式r(x)。（2）用程序验算D(x)Q(x)?r(x)?N(x)是否成立。

〖答案〗 Q =

0.6000 -1.4400 r =

-0.0000 0.0000 21.8800 -5.3400 -5.5200 4.5800 -2.8800

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

〖答案〗 P =

29

-0.0039 0.0338 -0.0227 -0.4456 0.9590 -0.0364 0.60.40.20-0.2-0.4-0.6-0.8-1-1.2-1.4-101234

19． 已知系统冲激响应为h(n)=[0.05,0.24,0.40,0.24,0.15,-0.1,0.1] ，系统输入u(n)由指令randn('state',1);u=2*(randn(1,100)>0.5)-1产生，该输入信号的起始作用时刻为0。试用直杆图（提示：用stem指令）画出分别显示该系统输入、输出信号的两张子图。

〖答案〗 Input u10.50-0.5-10204060Output y10.50-0.5-102040608010080100

30

本文来源：https://www.bwwdw.com/article/lc03.html