MATLAB--R2012a课后普习题答案全解

? ?

MATLAB R2012a 课后习题答案全解

第一章 基础准备及入门

习题1及解答

? 1.数字1.5e2，1.5e3 中的哪个与1500相同吗？

〖解答〗 1.5e3

? 2.请指出如下5个变量名中，哪些是合法的？

abcd-2 xyz_3 3chan a变量 ABCDefgh

〖解答〗

2、5是合法的。

? 3.在MATLAB环境中，比1大的最小数是多少？

〖解答〗 1+eps

? 4.设 a = -8 , 运行以下三条指令，问运行结果相同吗？为什么？

w1=a^(2/3) w2=(a^2)^(1/3) w3=(a^(1/3))^2 〖解答〗

（1）不同。具体如下

w1=a^(2/3)

%仅求出主根

1

w2=(a^2)^(1/3) w3=(a^(1/3))^2 w1 = w2 = 4.0000 w3 =

%求出(-8)^2的主根 %求出(-8)主根后再平方

-2.0000 + 3.4641i

-2.0000 + 3.4641i

（2）复数的多方根的，下面是求取全部方根的两种方法： （A）根据复数方根定义

a=-8;n=2;m=3;

ma=abs(a);aa=angle(a); for k=1:m

%m决定循环次数 %计算各根的相角 %计算各根

sa(k)=(aa+2*pi*(k-1))*n/m;

end

result=(ma^(2/3)).*exp(j*sa) result =

-2.0000 + 3.4641i 4.0000 - 0.0000i -2.0000 - 3.4641i

（B）利用多项式r?a?0求根

p=[1,0,0,-a^2]; r=roots(p) r =

-2.0000 + 3.4641i -2.0000 - 3.4641i 4.0000

32

? 5.指令clear, clf, clc各有什么用处？

〖解答〗 clear 清除工作空间中所有的变量。 clf 清除当前图形。 clc 清除命令窗口中所有显示。

? 6.以下两种说法对吗？（1）“MATLAB进行数值的表达精度与

其指令窗中的数据显示精度相同。” （2）MATLAB指令窗中显示的数值有效位数不超过7位。”

2

〖解答〗

（1）否；（2）否。

?123?456?，? 7.想要在MATLAB中产生二维数组S??下面哪些指令????789??能实现目的？

（1） S=[1,2,3;4,5,6;7,8;9]

（2） S=[1 2 3;4 5 6;7 8 9]

（3） S=[1，2，3；4，5，6；7，8，9] 〖解答〗

前两种输入方法可以，后一种方法不行。

%整个指令在中文状态下输入

? 8.试为例1.3-5编写一个解题用的M脚本文件？

〖解答〗

直接点击新文件图标，出现M文件编辑器窗口；在该M文件编辑器中，输入例1.3-5中的全部指令；并另存为p109.m，便得到所需的脚本文件。

第2章 符号运算

习题2及解答

? 1说出以下四条指令产生的结果各属于哪种数据类型，是“双精

度”对象，还是“符号”符号对象？

3/7+0.1; sym(3/7+0.1); sym('3/7+0.1'); vpa(sym(3/7+0.1))

〖目的〗

? 不能从显示形式判断数据类型，而必须依靠class指令。 〖解答〗

3

c1=3/7+0.1 c2=sym(3/7+0.1) c3=sym('3/7+0.1') c4=vpa(sym(3/7+0.1)) Cs1=class(c1) Cs2=class(c2) Cs3=class(c3) Cs4=class(c4) c1 = 0.5286 c2 = 37/70 c3 =

0.52857142857142857142857142857143 c4 =

0.52857142857142857142857142857143 Cs1 = double Cs2 = sym Cs3 = sym Cs4 = sym

? 2在不加专门指定的情况下，以下符号表达式中的哪一个变量被

认为是自由符号变量.

sym('sin(w*t)'),sym('a*exp(-X)'),sym('z*exp(j*th)')

〖目的〗

? 理解自由符号变量的确认规则。 〖解答〗

symvar(sym('sin(w*t)'),1) ans = w

symvar(sym('a*exp(-X)'),1) ans = a

4

symvar(sym('z*exp(j*th)'),1) ans = z

? 3求以下两个方程的解

3（1）试写出求三阶方程x?44.5?0正实根的程序。注意：只要正

实根，不要出现其他根。

（2）试求二阶方程x2?ax?a2?0在a?0时的根。

〖目的〗

? 体验变量限定假设的影响 〖解答〗

（1）求三阶方程x?44.5?0正实根

reset(symengine) syms x positive solve(x^3-44.5) ans =

(2^(2/3)*89^(1/3))/2

%确保下面操作不受前面指令运作的影响

3（2）求五阶方程x?ax?a?0的实根

syms a positive solve(x^2-a*x+a^2) > In solve at 83 ans =

[ empty sym ]

syms x clear syms a positive solve(x^2-a*x+a^2) ans =

a/2 + (3^(1/2)*a*i)/2 a/2 - (3^(1/2)*a*i)/2

%注意：关于x的假设没有去除

22Warning: Explicit solution could not be found.

5

f=sin(t)/t;

y=int(f,t,0,x) y5=subs(y,x,sym('4.5')) ezplot(y,[0,2*pi]) y = sinint(x) y5 =

1.6541404143792439835039224868515

sinint(x)1.81.61.41.210.80.60.40.200123x456% 将得到一个特殊经典函数

（2）数值计算复验

tt=0:0.001:4.5; tt(1)=eps;

yn=trapz(sin(tt)./tt)*0.001 yn =

1.6541

?? 12在n?0的限制下，求y(n)?13?20sinnxdx的一般积分表达式，

并计算y()的32位有效数字表达。

〖目的〗

? 一般符号解与高精度符号数值解。

11

〖解答〗

syms x

syms n positive f=sin(x)^n;

yn=int(f,x,0,pi/2)

y3s=vpa(subs(yn,n,sym('1/3'))) y3d=vpa(subs(yn,n,1/3)) yn =

beta(1/2, n/2 + 1/2)/2 y3s =

1.2935547796148952674767575125656 y3d =

1.2935547796148951782413405453553

? 13.有序列x(k)?ak，h(k)?bk，（在此k?0，a?b），求这两个序

列的卷积y(k)??h(n)x(k?n)。

n?0k〖目的〗

? 符号离散卷积直接法和变换法。 〖解答〗 （1）直接法

syms a b k n x=a^k; h=b^k;

w=symsum(subs(h,k,n)*subs(x,k,k-n),n,0,k) %据定义 y1=simple(w) w =

piecewise([a = b, b^k + b^k*k], [a <> b, (a*a^k - b*b^k)/(a - b)]) y1 =

piecewise([a = b, b^k + b^k*k], [a <> b, (a*a^k - b*b^k)/(a - b)])

（2）变换法（复验）

syms z

X=ztrans(a^k,k,z); H=ztrans(b^k,k,z); y2=iztrans(H*X,z,k) y2 =

%通过Z变换及反变换

piecewise([b <> 0, (a*a^k)/(a - b) - (b*b^k)/(a - b)])

〖说明〗

12

? 符号计算不同途径产生的结果在形式上有可能不同，而且往往无法依靠符号计算本身的

指令是它们一致。此时，必须通过手工解决。

t?0? 14.设系统的冲激响应为h(t)?e?3t，求该系统在输入u(t)?cost，

作用下的输出。

〖目的〗

? 符号连续函数卷积的直接法和变换法。 ? 符号变量限定性定义的作用。 ? laplace, ilaplace指令的应用。 〖解答〗 （1）直接法

syms t

h=exp(-3*t);u=cos(t); syms tao;

h_tao=subs(h,t,tao); u_t_tao=subs(u,t,t-tao); hu_tao=h_tao*u_t_tao;

hut=simple(int(hu_tao,tao,0,t)) %直接卷积 hut =

(3*cos(t))/10 - 3/(10*exp(3*t)) + sin(t)/10

（2）变换法（复验）

syms s;

HU=laplace(h,t,s)*laplace(u,t,s); huL=simple(ilaplace(HU,s,t)) huL =

%拉氏变换及反变换

(3*cos(t))/10 - 3/(10*exp(3*t)) + sin(t)/10

? 15.求f(t)?Ae??t,??0的Fourier变换。

〖目的〗

? 符号变量限定性定义的作用。 ? fourier指令的应用。 〖解答〗

syms A t w

a=sym('a','positive'); f=A*exp(-a*abs(t)); y=fourier(f,t,w) F=simple(y) y =

(2*A*a)/(a^2 + w^2) F =

13

(2*A*a)/(a^2 + w^2)

??t??A1?? 16.求f(t)??????0?????t??的Fourier变换，并画出A?2,??2时t??的幅频谱。

〖目的〗

? 单位阶跃符号函数heaviside的应用。 ? subs实现多变量置换。 ? ezplot的使用。 〖解答〗

syms t A w;

tao=sym('tao','positive');

f=A*((1+t/tao)*(heaviside(t+tao)-heaviside(t))+(1-t/tao)*(heaviside(t)-heaviside(t-tao))); Fw=fourier(f,t,w); Fws=simple(Fw)

Fw2=subs(Fws,[A,tao],[2,2]) ezplot(abs(Fw2)) grid Fws =

-(4*A*(cos((tao*w)/2)^2 - 1))/(tao*w^2) Fw2 =

-(8*cos(w)^2 - 8)/(2*w^2)

14

abs(8 cos(w)2 - 8)/(2 abs(w)2)43.532.521.510.50-6-4-20w246

? 17.求F(s)?〖解答〗

syms s t

s?3的Laplace反变换。

s3?3s2?6s?4F=(s+3)/(s^3+3*s^2+6*s+4); f=simple(ilaplace(F,s,t)) f =

(3^(1/2)*sin(3^(1/2)*t) - 2*cos(3^(1/2)*t) + 2)/(3*exp(t))

? 18.利用符号运算证明Laplace变换的时域求导性质：

?df(t)?L??s?L?f(t)??f(0)。 ??dt?〖目的〗

? 符号计算用于定理证明。 〖解答〗

syms t s; y=sym('f(t)'); df=diff(y,t);

Ldy=laplace(df,t,s)

15

Ldy =

s*laplace(f(t), t, s) - f(0)

?? k T? 19.求f(k)?ke的Z变换表达式。

〖目的〗

? 注意：变换中，被变换变量的约定。 〖解答〗

syms lambda k T z; f_k=k*exp(-lambda*k*T); F_z=simple(ztrans(f_k,k,z)) F_z =

(z*exp(T*lambda))/(z*exp(T*lambda) - 1)^2

? 20.求方程x2?y2?1,xy?2的解。

〖目的〗

? solve指令中，被解方程的正确书写，输出量的正确次序。 〖解答〗

eq1='x^2+y^2=1'; eq2='x*y=2';

[x,y]=solve(eq1,eq2,'x','y') x =

(1/2 + (15^(1/2)*i)/2)^(1/2)/2 - (1/2 + (15^(1/2)*i)/2)^(3/2)/2 - (1/2 + (15^(1/2)*i)/2)^(1/2)/2 + (1/2 + (15^(1/2)*i)/2)^(3/2)/2 (1/2 - (15^(1/2)*i)/2)^(1/2)/2 - (1/2 - (15^(1/2)*i)/2)^(3/2)/2 - (1/2 - (15^(1/2)*i)/2)^(1/2)/2 + (1/2 - (15^(1/2)*i)/2)^(3/2)/2 y =

(1/2 + (15^(1/2)*i)/2)^(1/2) -(1/2 + (15^(1/2)*i)/2)^(1/2) (1/2 - (15^(1/2)*i)/2)^(1/2) -(1/2 - (15^(1/2)*i)/2)^(1/2)

? 21.求图p2-1所示信号流图的系统传递函数，并对照胡寿松主编

“自动控制原理”中的例2-21结果，进行局部性验证。

16

图p2-1

〖目的〗

? 理解和掌握信号流图传递函数的“代数状态方程解法”。

? 并设法用胡寿松主编的“自动控制原理”的例2-21进行局部性验证。

〖解答〗

（1）求传递函数

syms G1 G2 G3 G4 G5 G6 G7 H1 H2 H3 H4 H5 A=[

0 0 0 0 -H3 -H4; 0 G2 0 0 -H2 G6;

G1 0 -H1 0 0 0; 0 0 G3 0 0 G7; 0 0 0 G4 0 0; 0 G5 0 0 0 -H5]; b=[ 1; 0; 0; 0; 0; 0]; c=[ 0 0 0 0 1 0]; Y2U=c*((eye(size(A))-A)\\b); %求传递函数 [NN,DD]=numden(Y2U); DD=sort(DD);

%分离出分子、分母多项式

%分母多项式排序

disp([blanks(5),'传递函数 Y2U 为']) pretty(NN/DD) 传递函数 Y2U 为

(G1 G4 (G2 G3 + G5 G7 + G3 G5 G6 + G2 G3 H5)) /

(H5 + G2 H1 + G3 G4 H2 + G1 G5 H4 + G5 G6 H1 + G2 H1 H5 + G3 G4 H2 H5 +

G1 G2 G3 G4 H3 + G1 G4 G5 G7 H3 - G4 G5 G7 H1 H2 + G1 G3 G4 G5 G6 H3

17

+

G1 G2 G3 G4 H3 H5 + G1 G3 G4 G5 H2 H4 + 1)

（2）局部性验证

syms a b c d e f g

y2u=subs(Y2U,[G1,G2,G3,G4,G5,G6,G7,H1,H2,H3,H4,H5],[a,e,f,1,b,c,0,g,0,0,0,d]);

[nn,dd]=numden(y2u); dd=sort(dd);

disp([blanks(5),'局部性验证用的传递函数y2u']) pretty(nn/dd)

局部性验证用的传递函数y2u

a (e f + b c f + d e f) --------------------------- d + e g + b c g + d e g + 1

此结果与胡寿松主编的“自动控制原理”例2-21一致。

? 22.采用代数状态方程法求图p2-2所示结构框图的传递函数

Y。 WY和U

图p2-2

〖目的〗

? 运用“代数状态方程解法”求输入和扰动同时存在的结构框图的传递函数。

〖解答〗

（1）理论演绎 对于结构框图写出状态方程

18

?x?Ax?bU?fW ??Y?cx?dU?gW此式第一个方程关于x的解可写为

（p2-1）

x?(I?A)?1bU?(I?A)?1fW

把此式代入式（p2-1）的第二个方程，加以整理后可得

（p2-2）

Y?c(I?A)?1b?dU?c(I?A)?1f?gW

据此可写出传递函数

????Y?c(I?A)?1b?d U

（p2-3） （p2-4）

Y?c(I?A)?1f?g N（2）列出“元素级”状态方程 值得提醒：在编写M码之前，最好先在草稿纸上，仔细“元素级”状态方程是避免出错的冲要措施。对此，不要掉以轻心。 本例的“元素级”状态方程如下

?x1??0?x??G?2??2?x3???0????x4??0??0?x5??? Y??01?G1??x1??G1??0??x??0??0?0?2???????0000???x3???0??U??G3??W （p2-5） ???????H1000??x4??0??0??H2000??0??????H2???x5???000??x?0?U?(?1)?W000?G2?G10（3）编写相应的M码

syms G1 G2 G3 H1 H2 A=[

0 0 0 -G1 -G1; 0 0 0 0 0;

G2 0 -G2 0 0; 0 H1 0 0 0; 0 H2 0 0 0]; b=[ G1; 0; 0; 0; 0]; f=[ 0; 0; G3; 0; -H2]; c=[ 0 1 0 0 0]; d=0; g=-1;

R=c/(eye(size(A))-A); Y2U=R*b+d; Y2W=R*f+g;

%中间变量 %计算传递函数 Y/U %计算传递函数 Y/W %分离出分子、分母多项式

%分母多项式排序

[NU,DU]=numden(Y2U); DU=sort(DU);

disp([blanks(5),'传递函数 Y2U 为'])

19

pretty(NU/DU) [NW,DW]=numden(Y2W); NW=sort(NW); DW=sort(DW);

disp([blanks(5),'传递函数 Y2W 为']) pretty(NW/DW) 传递函数 Y2U 为

G1 G2

----------------------- G1 G2 H1 + G1 G2 H2 + 1 传递函数 Y2W 为

G2 G3 + G1 G2 H1 + 1 - ----------------------- G1 G2 H1 + G1 G2 H2 + 1

? 23.求微分方程yy?5?x4?0的通解，并绘制任意常数为1时解的图

形。

〖目的〗

? 理解指令dsolve的正确使用。 ? 对dsolve输出结果的正确理解。

? ezplot指令绘图时，如何进行线色控制。 ? 如何覆盖那些不能反映图形窗内容的图名。 〖解答〗 （1）求通解

reset(symengine) clear syms y x

y=dsolve('0.2*y*Dy+0.25*x=0','x') y =

2^(1/2)*(C3 - (5*x^2)/8)^(1/2) -2^(1/2)*(C3 - (5*x^2)/8)^(1/2)

（2）根据所得通解中不定常数的符号写出“对其进行数值替代的指令”

yy=subs(y,'C3',1) yy =

2^(1/2)*(1 - (5*x^2)/8)^(1/2) -2^(1/2)*(1 - (5*x^2)/8)^(1/2)

%将通解中的C3用1代替

（3）观察通解中两个分解的平方是否相同 yy(1)^2==yy(2)^2

ans =

20

norm(A1-A2,'fro') % 求误差矩阵的F-范数，当接近eps量级时，就认为实际相等 B1 =

1.0000 1.4142 1.7321 2.0000 B2 =

0.5537 + 0.4644i 0.8070 - 0.2124i 1.2104 - 0.3186i 1.7641 + 0.1458i ans =

8.4961e-016

?0.5ty?1?ecos2t曲线。要求分别? 4.在时间区间 [0,10]中，绘制

采取“标量循环运算法”和“数组运算法”编写两段程序绘图。

〖目的〗

? 加强理解数组运算的机理和应用。

? 初步使用subplot, plot, xlabel, ylabel等指令绘图。 〖解答〗

%标量循环运算法

t=linspace(0,10,200); N=length(t); y1=zeros(size(t)); for k=1:N

y1(k)=1-exp(-0.5*t(k))*cos(2*t(k)); end

subplot(1,2,1),plot(t,y1),xlabel('t'),ylabel('y1'),grid on %数组运算法

y2=1-exp(-0.5*t).*cos(2*t);

subplot(1,2,2),plot(t,y2),xlabel('t'),ylabel('y2'),grid on

1.51.511y10.5y20.5005t10005t10

? 8.先运行clear,format long,rng('default'),A=rand(3,3)，然后根据

26

A写出两个矩阵：一个对角阵B，其相应元素由A的对角元素构成；另一个矩阵C，其对角元素全为0，而其余元素与对应的A阵元素相同。

〖目的〗

? 常用指令diag的使用场合。 〖解答〗

clear, format long rng('default') A=rand(3,3) B=diag(diag(A)) C=A-B A =

0.814723686393179 0.913375856139019 0.278498218867048 0.905791937075619 0.632359246225410 0.546881519204984 0.126986816293506 0.097540404999410 0.957506835434298 B =

0.814723686393179 0 0 0 0.632359246225410 0 0 0 0.957506835434298 C =

0 0.913375856139019 0.278498218867048 0.905791937075619 0 0.546881519204984 0.126986816293506 0.097540404999410 0

? 9.先运行指令x=-3*pi:pi/15:3*pi; y=x; [X,Y]=meshgrid(x,y);

warning off; Z=sin(X).*sin(Y)./X./Y; 产生矩阵Z。（1）请问矩阵Z中有多少个“非数”数据？（2）用指令surf(X,Y,Z); shading interp观察所绘的图形。（3）请写出绘制相应的“无裂缝”图形的全部指令。

〖目的〗

? 初步感受三维曲面的绘制方法。

? 非数NaN的产生，非数的检测，和对图形的影响。 ? sum的应用。

? eps如何克服“被零除”的尴尬。

27

〖解答〗

x=-3*pi:pi/15:3*pi; y=x;

[X,Y]=meshgrid(x,y); warning off

Z=sin(X).*sin(Y)./X./Y; NumOfNaN=sum(sum(isnan(Z))) %产生无缝图

XX=X+(X==0)*eps; YY=Y+(Y==0)*eps;

ZZ=sin(XX).*sin(YY)./XX./YY;

subplot(1,2,2),surf(XX,YY,ZZ),shading interp,title('无缝图') NumOfNaN = 181

%计算“非数”数目

subplot(1,2,1),surf(X,Y,Z),shading interp,title('有缝图')

? 10.下面有一段程序，企图用来解决如下计算任务：有矩阵

?1k?1?9k?1??2k?2?9k?2??，当k依次取10, 9, 8, 7, 6, 5, 4, 3, 2, 1时，Ak??????????k2k?10k??计算矩阵Ak“各列元素的和”，并把此求和结果存放为矩阵Sa

?14?28??25?29的第k行。例如k?3时，A阵为?此时它各列元素 的??，

??36?30??和是一个(1?10)行数组?615?87?，并把它保存为Sa的第3行。问题：该段程序的计算结果对吗？假如计算结果不正确，请指出

28

错误发生的根源，并改正之。

〖目的〗

? 正确理解sum的工作机理。 ? reshape的应用。 〖解答〗

（1）企图用以下程序完成题目要求。

for k=10:-1:1

A=reshape(1:10*k,k,10); Sa(k,:)=sum(A);

end Sa Sa =

55 55 55 55 55 55 55 55 55 55 3 7 11 15 19 23 27 31 35 39 6 15 24 33 42 51 60 69 78 87 10 26 42 58 74 90 106 122 138 154 15 40 65 90 115 140 165 190 215 240 21 57 93 129 165 201 237 273 309 345 28 77 126 175 224 273 322 371 420 469 36 100 164 228 292 356 420 484 548 612 45 126 207 288 369 450 531 612 693 774 55 155 255 355 455 555 655 755 855 955

（2）正确性分析

除k=1外，计算所得Sa所有行的结果都正确。但k=1时，A1?[1,2,?,10]，Sa的第

1行应该与A1相同。

上述程序的错误是对sum理解不正确。sum对二维数组，求和按列施行；而对一维数组，不管行数组或列数组，总是求那数组所有元素的和。 正确的程序应该写成

for k=10:-1:1

A=reshape(1:10*k,k,10); Sa(k,:)=sum(A); if k==1

Sa(k,:)=A; end

end Sa Sa =

1 2 3 4 5 6 7 8 9 10 3 7 11 15 19 23 27 31 35 39 6 15 24 33 42 51 60 69 78 87

29

10 26 42 58 74 90 106 122 138 154 15 40 65 90 115 140 165 190 215 240 21 57 93 129 165 201 237 273 309 345 28 77 126 175 224 273 322 371 420 469 36 100 164 228 292 356 420 484 548 612 45 126 207 288 369 450 531 612 693 774 55 155 255 355 455 555 655 755 855 955

第4章 数值运算

习题 4 及解答

? 1.根据题给的模拟实际测量数据的一组t和 y(t)试用数值差分diff

或数值梯度gradient指令计算y?(t)，然后把y(t)和y?(t)曲线绘制在同一张图上，观察数值求导的后果。（模拟数据从prob_data401.mat获得）

〖目的〗

? 强调：要非常慎用数值导数计算。 ? 练习mat数据文件中数据的获取。 ? 实验数据求导的后果

? 把两条曲线绘制在同一图上的一种方法。

〖解答〗

（1）从数据文件获得数据的指令 假如prob_data401.mat文件在当前目录或搜索路径上

clear

load prob_data401.mat

（2）用diff求导的指令

dt=t(2)-t(1);

30

yc=diff(y)/dt; %注意yc的长度将比y短1 plot(t,y,'b',t(2:end),yc,'r') grid on

1.510.50-0.5-1-1.5-201234567

（3）用gradent求导的指令（图形与上相似）

dt=t(2)-t(1); yc=gradient(y)/dt; plot(t,y,'b',t,yc,'r') grid on

〖说明〗

? 不到万不得已，不要进行数值求导。 ? 假若一定要计算数值导数，自变量增量dt 要取得比原有数据相对误差高1、2个量级以

上。

? 求导会使数据中原有的噪声放大。

? 2.采用数值计算方法，画出y(x)?并计算y(4.5)。

?x0sintdt在[0, 10]区间曲线，t〖提示〗

? 指定区间内的积分函数可用cumtrapz指令给出。 ?

y(4.5)在计算要求不太高的地方可用find指令算得。

〖目的〗

? 指定区间内的积分函数的数值计算法和cumtrapz指令。

31

? find指令的应用。 〖解答〗

dt=1e-4; t=0:dt:10; t=t+(t==0)*eps; f=sin(t)./t; s=cumtrapz(f)*dt; plot(t,s,'LineWidth',3) ii=find(t==4.5); s45=s(ii) s45 = 1.6541

21.81.61.41.210.80.60.40.200246810

? 3.求函数f(x)?e尝试复算。

〖提示〗

? 数值积分均可尝试。 ? 符号积分的局限性。

〖目的〗

? 符号积分的局限性。 〖解答〗

dx=pi/2000; x=0:dx:pi;

sin3x的数值积分s?? 0f(x)dx，并请采用符号计算

?s=trapz(exp(sin(x).^3))*dx

32

s =

5.1370

符号复算的尝试

syms x

f=exp(sin(x)^3); ss=int(f,x,0,pi)

Warning: Explicit integral could not be found. > In sym.int at 58 ss =

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

? 4.用quad求取??5?e?xsinxdx的数值积分，并保证积分的绝对精

度为10?9。

〖目的〗

? quadl，精度可控，计算较快。

? 近似积分指令trapz获得高精度积分的内存和时间代价较高。 〖解答〗

%精度可控的数值积分

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

sq=quadl(fx,-10*pi,1.7*pi,1e-7) sq =

1.08784993815498

1.7?

%近似积分算法

x=linspace(-10*pi,1.7*pi,1e7); dx=x(2)-x(1);

st=trapz(exp(-abs(x)).*abs(sin(x)))*dx st =

1.08784949973430

%符号积分算法

y='exp(-abs(x))*abs(sin(x))' si=vpa(int(y,-10*pi,1.7*pi),16) y =

exp(-abs(x))*abs(sin(x)) si =

1.087849499412911

33

? 5.求函数f(t)?(sin5t)2e0.06t?1.5tcos2t?1.8t?0.5在区间[?5,5]中的最

小值点。

〖目的〗

? 理解极值概念的邻域性。 ? 如何求最小值。

? 学习运用作图法求极值或最小值。 ? 感受符号法的局限性。 〖解答〗

（1）采用fminbnd找极小值点 在指令窗中多次运行以下指令，观察在不同数目子区间分割下，进行的极小值搜索。然后从一系列极小值点中，确定最小值点。 clear

ft=@(t)sin(5*t).^2.*exp(0.06*t.*t)+1.8*abs(t+0.5)-1.5*t.*cos(2*t);

disp('计算中，把[ -5,5] 分成若干搜索子区间。')

N=input(' 请输入子区间数 N，注意使N>=1 ？');%该指令只能在指令窗中运行 tt=linspace(-5,5,N+1); for k=1:N

[tmin(k),fobj(k)]=fminbnd(ft,tt(k),tt(k+1)); end

[fobj,ii]=sort(fobj); %将目标值由小到大排列 tmin=tmin(ii); %使极小值点做与目标值相应的重新排列 fobj,tmin

（2）最后确定的最小值点 在N?1,2,?,10的不同分割下，经观察，最后确定出 最小值点是 -1.28498111480531 相应目标值是 -0.18604801006545

（3）采用作图法近似确定最小值点（另一方法） （A）在指令窗中运行以下指令： clear

ft=@(t)sin(5*t).^2.*exp(0.06*t.*t)+1.8*abs(t+0.5)-1.5*t.*cos(2*t);

t=-5:0.001:5; ff=ft(t); plot(t,ff) grid on,shg

（B）经观察后，把最小值附近邻域放到足够大，然后运行以下指令，那放大图形被推向前台，与此同时光标变为“十字线”，利用它点击极值点可得到最小值数据

[tmin2,fobj2]=ginput(1)

2 34

tmin2 =

-1.28500000993975 fobj2 =

-0.18604799369136

出现具有相同数值的刻度区域表明已达最小可分辨状态

（4）符号法求最小值的尝试

syms t

fts=sin(5*t)^2*exp(0.06*t*t)-1.5*t*cos(2*t)+1.8*abs(t+0.5); dfdt=diff(fts,t); tmin=solve(dfdt,t)

%求导函数

%求导函数的零点

fobj3=subs(fts,t,tmin) %得到一个具体的极值点 tmin =

-.60100931947716486053884417850955e-2 fobj3 =

.89909908144684551670208797723124

〖说明〗

35

? 最小值是对整个区间而言的，极小值是对邻域而言的。 ? 在一个区间中寻找最小值点，对不同子区间分割进行多次搜索是必要的。这样可以避免

把极小值点误作为最小值点。最小值点是从一系列极小值点和边界点的比较中确定的。 ? 作图法求最小值点，很直观。假若绘图时，自变量步长取得足够小，那么所求得的最小

值点有相当好的精度。 ? 符号法在本例中，只求出一个极值点。其余很多极值点无法秋初，更不可能得到最小值。

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

〖目的〗

? 学习如何把高阶微分方程写成一阶微分方程组。

? ode45解算器的导数函数如何采用匿名函数形式构成。

? 如何从ode45一组数值解点，求指定自变量对应的函数值。 〖解答〗

（1）改写高阶微分方程为一阶微分方程组

令y1(t)?y(t),y2(t)?dy(t)，于是据高阶微分方程可写出 dtdy1(t)??y2(t)?dt ?dy(t)?2??2y1(t)?3y2(t)?1?dt

（2）运行以下指令求y(t)的数值解

format long ts=[0,1]; y0=[1;0];

dydt=@(t,y)[y(2);-2*y(1)+3*y(2)+1];

%<4>

%匿名函数写成的ode45所需得导数函数

[tt,yy]=ode45(dydt,ts,y0);

y_05=interp1(tt,yy(:,1),0.5,'spline'), %用一维插值求y(0.5) y_05 =

0.78958020790127

（3）符号法求解

syms t;

ys=dsolve('D2y-3*Dy+2*y=1','y(0)=1,Dy(0)=0','t') ys_05=subs(ys,t,sym('0.5')) ys =

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

36

.78958035647060552916850705213780

〖说明〗

? 第<4>条指令中的导数函数也可采用M函数文件表达，具体如下。

function S=prob_DyDt(t,y) S=[y(2);-2*y(1)+3*y(2)+1];

? 7.已知矩阵A=magic(8)，（1）求该矩阵的“值空间基阵”B ；

（2）写出“A的任何列可用基向量线性表出”的验证程序（提示：利用rref检验）。

〖目的〗

? 体验矩阵值空间的基向量组的不唯一性，但它们可以互为线性表出。 ? 利用rref检验两个矩阵能否互为表出。 〖解答〗

（1）A的值空间的三组不同“基”

A=magic(8); [R,ci]=rref(A); B1=A(:,ci)

B2=orth(A) [V,D]=eig(A); B3=V(:,1:rv) 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 -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 =

%采用8阶魔方阵作为实验矩阵

%直接从A中取基向量 %求A值空间的正交基

%非零特征值数就是矩阵的秩

%取A的非零特征值对应的特征向量作基

rv=sum(sum(abs(D))>1000*eps);

37

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

（2）验证A的任何列可用B1线性表出

B1_A=rref([B1,A]) %若B1_A矩阵的下5行全为0， B1_A =

%就表明A可以被B1的3根基向量线性表出

1 0 0 1 0 0 1 1 0 0 1 0 1 0 0 1 0 3 4 -3 -4 7 0 0 1 0 0 1 -3 -4 4 5 -7 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

B2_A=rref([B2,A]) B2_A =

Columns 1 through 7

1.0000 0 0 -91.9239 -91.9239 -91.9239 -91.9239 0 1.0000 0 51.8459 -51.8459 -51.8459 51.8459 0 0 1.0000 9.8995 -7.0711 -4.2426 1.4142 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 Columns 8 through 11

-91.9239 -91.9239 -91.9239 -91.9239 51.8459 -51.8459 -51.8459 51.8459 -1.4142 4.2426 7.0711 -9.8995 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

38

B3_A=rref([B3,A]) B3_A =

Columns 1 through 7

1.0000 0 0 91.9239 91.9239 91.9239 91.9239 0 1.0000 0 42.3447 -38.1021 -33.8594 29.6168 0 0 1.0000 12.6462 -16.8889 -21.1315 25.3741 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 Columns 8 through 11

91.9239 91.9239 91.9239 91.9239 25.3741 -21.1315 -16.8889 12.6462 29.6168 -33.8594 -38.1021 42.3447 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

〖说明〗

? magic(n)产生魔方阵。魔方阵具有很多特异的性质。就其秩而言，当n为奇数时，该矩

阵满秩；当n 为4的倍数时，矩阵的秩总是3；当 n 为偶数但不是4倍数时，则矩阵的秩等于 (n/2+2)。关于魔方阵的有关历史，请见第6.1.3节。

? 8.已知由MATLAB指令创建的矩阵A=gallery(5)，试对该矩阵进

行特征值分解，并通过验算观察发生的现象。

〖目的〗

? 展示特征值分解可能存在的数值问题。 ? condeig是比较严谨的特征值分解指令。 ? Jordan分解的作用。 〖解答〗

（1）特征值分解

A=gallery(5) [V,D]=eig(A); diag(D)' A =

%为紧凑地显示特征值而写

-9 11 -21 63 -252 70 -69 141 -421 1684 -575 575 -1149 3451 -13801 3891 -3891 7782 -23345 93365

39

1024 -1024 2048 -6144 24572 ans =

Columns 1 through 4

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

0.0144 + 0.0104i

（2）验算表明相对误差较大

AE=V*D/V

er_AE=norm(A-AE,'fro')/norm(A,'fro') AE =

1.0e+004 *

Columns 1 through 4

-0.0009 + 0.0000i 0.0011 - 0.0000i -0.0021 + 0.0000i 0.0063 - 0.0000i

0.0070 - 0.0000i -0.0069 + 0.0000i 0.0141 - 0.0000i -0.0421 + 0.0000i

-0.0575 + 0.0000i 0.0575 - 0.0000i -0.1149 + 0.0000i 0.3451 - 0.0000i

0.3891 - 0.0000i -0.3891 + 0.0000i 0.7781 - 0.0000i -2.3343 + 0.0000i

0.1024 - 0.0000i -0.1024 + 0.0000i 0.2048 - 0.0000i -0.6144 + 0.0000i Column 5

-0.0252 + 0.0000i 0.1684 - 0.0000i -1.3800 + 0.0000i 9.3359 - 0.0001i 2.4570 - 0.0000i er_AE =

6.9310e-005

%相对F范数

（3）一个更严谨的特征值分解指令

[Vc,Dc,eigc]=condeig(A) Vc =

Columns 1 through 4

-0.0000 -0.0000 + 0.0000i -0.0000 - 0.0000i 0.0000 + 0.0000i

0.0206 0.0207 + 0.0000i 0.0207 - 0.0000i 0.0207 + 0.0000i -0.1397 -0.1397 + 0.0000i -0.1397 - 0.0000i -0.1397 + 0.0000i

0.9574 0.9574 0.9574 0.9574 0.2519 0.2519 - 0.0000i 0.2519 + 0.0000i 0.2519 - 0.0000i

%eigc中的高值时，说明相应的特征值不可信。

40

Column 5

0.0000 - 0.0000i 0.0207 - 0.0000i -0.1397 - 0.0000i 0.9574 0.2519 + 0.0000i Dc =

Columns 1 through 4

-0.0181 0 0 0 0 -0.0054 + 0.0171i 0 0 0 0 -0.0054 - 0.0171i 0 0 0 0 0.0144 + 0.0104i 0 0 0 0 Column 5 0 0 0 0 0.0144 - 0.0104i eigc = 1.0e+011 * 5.2687 5.2313 5.2313 5.1725 5.1724

（4）对A采用Jordan分解并检验

[VJ,DJ]=jordan(A); DJ

AJ=VJ*DJ/VJ

er_AJ=norm(A-AJ,'fro')/norm(A,'fro') 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 AJ =

1.0e+004 *

-0.0009 0.0011 -0.0021 0.0063 -0.0252 0.0070 -0.0069 0.0141 -0.0421 0.1684 -0.0575 0.0575 -0.1149 0.3451 -1.3801 0.3891 -0.3891 0.7782 -2.3345 9.3365

%求出准确的广义特征值，使A*VJ=VJ*D成立。

41

0.1024 -0.1024 0.2048 -0.6144 2.4572 er_AJ =

2.0500e-011

〖说明〗

? 指令condeig的第3输出量eigc给出的是所谓的“矩阵特征值条件数”。当特征条件数

与1eps相当时，就意味着矩阵A可能“退化”，即矩阵可能存在“代数重数”大于

“几何重数”的特征值。此时，实施Jordan分解更适宜。 ? 顺便指出：借助condeig算得的特征值条件数与cond指令算得的矩阵条件数是两个不同

概念。前者描述特征值的问题，后者描述矩阵逆的问题。

? 本例矩阵A的特征值条件数很高，表明分解不可信。验算也表明，相对误差较大。 ? 当对矩阵A进行Jordan分解时，可以看到，A具有5重根。当对Jordan分解进行验算

时，相对误差很小。

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

〖提示〗

? 了解magic指令 ? rref 用于方程求解。

? 矩阵除法和逆阵法解方程。

〖目的〗

? 满秩方阵求解的一般过程。 ? rref 用于方程求解。

? 矩阵除法和逆阵法解方程。 〖解答〗

A=magic(3);

%产生3阶魔方阵 %（3*1）全1列向量 %矩阵除解方程 %逆阵法解方程

b=ones(3,1); x=A\\b

[R,C]=rref([A,b]) %Gauss Jordan消去法解方程，同时判断解的唯一性 xx=inv(A)*b R =

1.0000 0 0 0.0667 0 1.0000 0 0.0667 0 0 1.0000 0.0667 C =

1 2 3 x = 0.0667 0.0667 0.0667 xx = 0.0667

42

0.0667 0.0667

〖说明〗

? rref指令通过对增广矩阵进行消去法操作完成解方程。由分解得到的3根“坐标向量”

和（或）C3指示的3根基向量，可见A3满秩，因此方程解唯一。 ? 在本例情况下，矩阵除、逆阵法、rref法所得解一致。

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

〖提示〗

? 用rref 可观察A不满秩，但b在A的值空间中，这类方程用无数解。 ? 矩阵除法能正确求得这类方程的特解。 ? 逆阵法不能求得这类方程的特解。 ? 注意特解和齐次解

〖目的〗

? A不满秩，但b在A的值空间中，这类方程的求解过程。 ? rref 用于方程求解。

? 矩阵除法能正确求得这类方程的特解。 ? 逆阵法不能求得这类方程的特解。 ? 解的验证方法。 ? 齐次解的获取。 ? 全解的获得。 〖解答〗

（1）借助增广矩阵用指令rref求解

A=magic(4);

%产生3阶魔方阵 %全1列向量

b=ones(4,1);

[R,C]=rref([A,b]) R =

%求解，并判断解的唯一性

1.0000 0 0 1.0000 0.0588 0 1.0000 0 3.0000 0.1176 0 0 1.0000 -3.0000 -0.0588 0 0 0 0 0 C =

1 2 3

关于以上结果的说明： ? R阶梯阵提供的信息

? 前4列是原A阵经消元变换后的阶梯阵；而第5列是原b向量经相同变换后的结

果。

? R的前三列为“基”，说明原A阵秩为3；而第4列的前三个元素，表示原A阵

的第4列由其前三列线性组合而成时的加权系数，即方程的一个解。

43

? R的第5列表明：b可由原A阵的前三列线性表出；b给出了方程的一个解；由于

原A阵“缺秩”，所以方程的确解不唯一。

? C数组提供的信息

? 该数组中的三个元素表示变换取原A阵的第1，2，3列为基。 ? 该数组的元素总数就是“原A阵的秩”

（2）矩阵除求得的解

x=A\\b

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.306145e-017. x = 0.0588 0.1176 -0.0588 0

运行结果指示：矩阵除法给出的解与rref解相同。（实际上，MATLAB在设计“除法”程序时，针对“b在A值空间中”的情况，就是用rref求解的。）

（3）逆阵法的解

xx=inv(A)*b

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.306145e-017. xx = 0.0469 0.1875 -0.0625 -0.0156

（3）验证前面所得的解

b_rref=A(:,C)*R(C,5) %验算rref的解 b_d=A*x

%验算矩阵除的解

b_inv=A*xx b_rref = 1 1 1 1 b_d = 1 1 1 1 b_inv =

%验算逆阵法的解

44

0.7344 1.5469 1.1719 1.8594

显然，在本例中，逆阵法的解是错误的。原因是：A不满秩，A的逆阵在理论上不存在。这里所给出的逆阵是不可信的。

（4）求齐次解

xg=null(A) xg = 0.2236 0.6708 -0.6708 -0.2236

%Ax=0的齐次解

（5）方程的全解 齐次解的任何倍与特解之和就构成方程的全解。下面通过一组随机系数验证。

rng default

%为本书结果可被读者核对而设，并非必要。 %6个随机系数 %产生6个不同的特解

%所得结果的每列都应该是全1，即等于b.

f=randn(1,6) A*xx f = xx =

0.1790 0.4689 -0.4463 0.2516 0.1301 -0.2336 0.4783 1.3479 -1.3976 0.6960 0.3315 -0.7596 -0.4195 -1.2890 1.4565 -0.6372 -0.2727 0.8184 -0.1202 -0.4101 0.5051 -0.1928 -0.0713 0.2924 ans =

1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000

xx=repmat(x,1,6)+xg*f

0.5377 1.8339 -2.2588 0.8622 0.3188 -1.3077

〖解答〗

? （在用除法和逆阵法求解时出现）警告信息中RCOND = 1.306145e-017是矩阵A的估

计条件倒数。该数愈接近0，A就愈“病态”；该数接近1时，A就愈“良态”。该条件数由rcond(A)给出。注意：rcond条件倒数与cond条件数的算法不同。

45

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