matlab课后习题答案（1-9章）

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

w1 = -2.0000 + 3.4641i ;w2 = 4.0000 ;w3 =-2.0000 + 3.4641i 5 指令clear, clf, clc各有什么用处？

clear 清除工作空间中所有的变量。 clf 清除当前图形。clc 第二章

清除命令窗口中所有显示。

1 说出以下四条指令产生的结果各属于哪种数据类型，是“双精度”对象，还是“符号”符号对象？

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

2 在不加专门指定的情况下，以下符号表达式中的哪一个变量被认为是自由符号变量. sym('sin(w*t)'),sym('a*exp(-X)'),sym('z*exp(j*th)') symvar(sym('sin(w*t)'),1) w a z

3 （1）试写出求三阶方程x?44.5?0正实根的程序。注意：只要正实根，不要出现其他根。 （2）试求二阶方程x2?ax?a2?0在a?0时的根。 （1）reset(symengine) syms x positive solve(x^3-44.5) ans =

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

223

（2）求五阶方程x?ax?a?0的实根 syms a positive %注意：关于x的假设没有去除 solve(x^2-a*x+a^2)

Warning: Explicit solution could not be found. > 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

4 观察一个数（在此用@记述）在以下四条不同指令作用下的异同。

a =@, b = sym( @ ), c = sym( @ ,'d ' ), d = sym( '@ ' )

在此，@ 分别代表具体数值 7/3 , pi/3 , pi*3^(1/3) ；而异同通过vpa(abs(a-d)) , vpa(abs(b-d)) , vpa(abs(c-d))等来观察。

? 理解准确符号数值的创建法。 ? 高精度误差的观察。

（1）x=7/3

x=7/3;a=x,b=sym(x),c=sym(x,'d'),d=sym('7/3'), a =

2.3333 b = 7/3 c =

1

2.3333333333333334813630699500209 d = 7/3

v1=vpa(abs(a-d)),v2=vpa(abs(b-d)),v3=vpa(abs(c-d)) v1 = 0.0 v2 = 0.0 v3 =

0.00000000000000014802973661668756666666667788716

（2）x=pi/3

x=pi/3;a=x,b=sym(x),c=sym(x,'d'),d=sym('pi/3'), a =

1.0472 b = pi/3 c =

1.047197551196597631317786181171 d = pi/3

v1=vpa(abs(a-d)),v2=vpa(abs(b-d)),v3=vpa(abs(c-d)) v1 = 0.0 v2 = 0.0 v3 =

0.00000000000000011483642827992216762806615818554

（3）x=pi*3^(1/3)

x=pi*3^(1/3);a=x,b=sym(x),c=sym(x,'d'),d=sym('pi*3^(1/3)') a =

4.5310 b =

1275352044764433/281474976710656 c =

4.5309606547207899041040946030989 d =

pi*3^(1/3)

v1=vpa(abs(a-d)),v2=vpa(abs(b-d)),v3=vpa(abs(c-d)) v1 =

0.00000000000000026601114166290944374842393221638 v2 =

0.00000000000000026601114166290944374842393221638 v3 =

0.0000000000000002660111416629094726767991785515

?a11?5 求符号矩阵A?a21???a31a12a22a32a13??a23的行列式值和逆，所得结果应采用“子表达式置换”简洁化。

?a33??

? 理解subexpr指令。

A=sym('[a11 a12 a13;a21 a22 a23;a31 a32 a33]') DA=det(A) IA=inv(A);

[IAs,d]=subexpr(IA,d) A =

[ a11, a12, a13] [ a21, a22, a23] [ a31, a32, a33] DA =

2

a11*a22*a33 - a11*a23*a32 - a12*a21*a33 + a12*a23*a31 + a13*a21*a32 - a13*a22*a31 IAs =

[ d*(a22*a33 - a23*a32), -d*(a12*a33 - a13*a32), d*(a12*a23 - a13*a22)] [ -d*(a21*a33 - a23*a31), d*(a11*a33 - a13*a31), -d*(a11*a23 - a13*a21)] [ d*(a21*a32 - a22*a31), -d*(a11*a32 - a12*a31), d*(a11*a22 - a12*a21)] d =

1/(a11*a22*a33 - a11*a23*a32 - a12*a21*a33 + a12*a23*a31 + a13*a21*a32 - a13*a22*a31)

??6 求 ? ?

?k?0xk的符号解，并进而用该符号解求

?k?0(?13?)k，

?k?0(1??)k，

?k?03k的准确值。

symsum, subs的应用。

?从实例中，感受指令所给出的关于

?k?0xk符号解的含义。

syms x k f=x^(k);

Z1=symsum(f,k,0,inf) Z1 =

piecewise([1 <= x, Inf], [abs(x) < 1, -1/(x - 1)])

subs(Z1,x,{sym('-1/3'),sym('1/pi'),sym('3')}) ans =

[ 3/4, -1/(1/pi - 1), Inf]

?7 对于x?0，求

?k?0?x?1???2k?1?x?1?22k?1。（提示：理论结果为lnx）

? 符号变量的限定性定义的作用。

syms k;

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

f_k=2/(2*k+1)*((x-1)/(x+1))^(2*k+1);

s=simple(symsum(f_k,k,0,inf)) %结果与理论值lnx相符！ s =

piecewise([abs(x - 1) < x + 1, log(x)])

〖注意〗

? 解答中，条件abs(x - 1) < x + 1意味着：

? 2>0 ? ? 约束一：x-10 此为“约束”，满足题意。 8 （1）通过符号计算求y(t)?sint的导数

dydtdy。（2）然后根据此结果，求

dy和

t?0?dtdtt??2。

? diff, limit指令的应用。 ? 如何理解运行结果。

syms t

y=abs(sin(t))

d=diff(y) %求dy/dt

d0_=limit(d,t,0,'left') %求dy/dt|t=0- dpi_2=limit(d,t,pi/2) %求dy/dt|t=pi/2 y =

abs(sin(t)) d =

3

sign(sin(t))*cos(t) d0_ = -1

dpi_2 = 0 9 求出

?1.7??10?e?xsinxdx的具有64位有效数字的积分值。

? 符号积分的解析解和符号数值解。 ? 符号计算和数值计算的相互校验。

（1）符号积分 syms x clear syms x

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

si=vpa(int(y,-10*pi,1.7*pi),64) y =

abs(sin(x))/exp(abs(x)) si =

1.087849499412904913166671875948174520895458535212845987519414166 （2）数值计算复验

xx=-10*pi:pi/100:1.7*pi;

sn=trapz(exp(-abs(xx)).*abs(sin(xx)))*pi/100 sn =

1.0877

10 计算二重积分

??12x21(x?y)dydx。

22

? 变上限二重积分的符号计算法。

syms x y f=x^2+y^2;

r=int(int(f,y,1,x^2),x,1,2) r =

1006/105

11 在[0,2?]区间，画出y(x)??xsintt0dt曲线，并计算y(4.5)。

? 在符号计算中，经常遇到计算结果是特殊经典函数的情况。 ? 如何应用subs获得超过16位有效数字的符号数值结果。 ? 初步尝试ezplot指令的简便。

（1）符号计算 syms t x; 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

4

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)??20sinnxdx的一般积分表达式，并计算y()的32位有效

13

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

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 求方程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)

5

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 =

.78958035647060552916850705213780

〖说明〗

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

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

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

〖提示〗

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

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

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

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

A=magic(3); %产生3阶魔方阵 b=ones(3,1); %（3*1）全1列向量

[R,C]=rref([A,b]) %Gauss Jordan消去法解方程，同时判断解的唯一性 x=A\\b %矩阵除解方程 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 0.0667 0.0667

〖说明〗

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

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

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

〖提示〗

16

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

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

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

（1）借助增广矩阵用指令rref求解 A=magic(4); %产生3阶魔方阵 b=ones(4,1); %全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列由其前

三列线性组合而成时的加权系数，即方程的一个解。 ? 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 %验算矩阵除的解

17

b_inv=A*xx %验算逆阵法的解 b_rref = 1 1 1 1 b_d = 1 1 1 1 b_inv = 0.7344 1.5469 1.1719 1.8594 显然，在本例中，逆阵法的解是错误的。原因是：A不满秩，A的逆阵在理论上不存在。这里所给出的逆阵是不可信的。

（4）求齐次解 xg=null(A) %Ax=0的齐次解 xg =

0.2236 0.6708 -0.6708 -0.2236

（5）方程的全解 齐次解的任何倍与特解之和就构成方程的全解。下面通过一组随机系数验证。 rng default %为本书结果可被读者核对而设，并非必要。 f=randn(1,6) %6个随机系数 xx=repmat(x,1,6)+xg*f %产生6个不同的特解 A*xx %所得结果的每列都应该是全1，即等于b. f =

0.5377 1.8339 -2.2588 0.8622 0.3188 -1.3077 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

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

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

? 求?0.5?t?10e?0.2tsin[sint]?0的实数解。

〖提示〗

? 在适当范围内，作图观察一元复杂函数的形态：观察解的存在性；解的唯一性。进而，借助图形法求

近似解。

? 匿名函数的使用方法。 ? fzero指令的用法。

? 作图法求一元复杂函数解上的作用：观察解的存在性；解的唯一性；得近似解。 ? 匿名函数的使用方法。

18

? fzero指令的用法。

（1）作图观察函数并求近似解 t=-1:0.001:5;

y=@(t)-0.5+t-10*exp(-0.2*t).*abs(sin(sin(t))); plot(t,y(t)) %利用匿名函数求y函数值 grid on,shg

[tt1,yy1]=ginput(1) %从图形获得近似解 tt1 =

2.7370 yy1 =

0.0097 420-2-4-6-8-10-12-1012345 （2）进一步利用fzero求精确解 [t,yt]=fzero(y,tt1) t =

2.7341 yt =

2.2204e-015

〖说明〗

? 假如在从图上获取数据前，先把零点附近图形放大，可以得到精度更高的近似解。

19

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