使用MATLAB进行数值计算

与符号计算相比，数值计算在科研和工程中的应用更为广泛。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

\"Banana\" 测试函数，它的理论极小值是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指令时两种不同序列记述法所得的卷积序列

23