matlab课后习题答案第四章.docx
《matlab课后习题答案第四章.docx》由会员分享，可在线阅读，更多相关《matlab课后习题答案第四章.docx（39页珍藏版）》请在冰点文库上搜索。
matlab课后习题答案第四章
第4章数值运算
习题4及解答
11根据题给的模拟实际测量数据的一组
和
试用数值差分diff或数值梯度gradient指令计算
，然后把
和
曲线绘制在同一张图上，观察数值求导的后果。
（模拟数据从prob_data401.mat获得）
〖目的〗
●强调：
要非常慎用数值导数计算。
●练习mat数据文件中数据的获取。
●实验数据求导的后果
●把两条曲线绘制在同一图上的一种方法。
〖解答〗
（1）从数据文件获得数据的指令
假如prob_data401.mat文件在当前目录或搜索路径上
clear
loadprob_data401.mat
（2）用diff求导的指令
dt=t
（2）-t
（1）;
yc=diff（y）/dt;%注意yc的长度将比y短1
plot（t,y,'b',t（2:
end）,yc,'r'）
gridon
（3）用gradent求导的指令（图形与上相似）
dt=t
（2）-t
（1）;
yc=gradient（y）/dt;
plot（t,y,'b',t,yc,'r'）
gridon
〖说明〗
●不到万不得已，不要进行数值求导。
●假若一定要计算数值导数，自变量增量dt要取得比原有数据相对误差高1、2个量级以上。
●求导会使数据中原有的噪声放大。
12采用数值计算方法，画出
在
区间曲线，并计算
。
〖提示〗
●指定区间内的积分函数可用cumtrapz指令给出。
●
在计算要求不太高的地方可用find指令算得。
〖目的〗
●指定区间内的积分函数的数值计算法和cumtrapz指令。
●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
13求函数
的数值积分
，并请采用符号计算尝试复算。
〖提示〗
●数值积分均可尝试。
●符号积分的局限性。
〖目的〗
●符号积分的局限性。
〖解答〗
dx=pi/2000;
x=0:
dx:
pi;
s=trapz（exp（sin（x）.^3））*dx
s=
5.1370
符号复算的尝试
symsx
f=exp（sin（x）^3）;
ss=int（f,x,0,pi）
Warning:
Explicitintegralcouldnotbefound.
>Insym.intat58
ss=
int（exp（sin（x）^3）,x=0..pi）
14用quad求取
的数值积分，并保证积分的绝对精度为
。
〖目的〗
●quadl，精度可控，计算较快。
●近似积分指令trapz获得高精度积分的内存和时间代价较高。
〖解答〗
%精度可控的数值积分
fx=@（x）exp（-abs（x））.*abs（sin（x））;
formatlong
sq=quadl（fx,-10*pi,1.7*pi,1e-7）
sq=
1.08784993815498
%近似积分算法
x=linspace（-10*pi,1.7*pi,1e7）;
dx=x
（2）-x
（1）;
st=trapz（exp（-abs（x））.*abs（sin（x）））*dx
st=
.0878********
%符号积分算法
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
15求函数
在区间
中的最小值点。
〖目的〗
●理解极值概念的邻域性。
●如何求最小值。
●学习运用作图法求极值或最小值。
●感受符号法的局限性。
〖解答〗
（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）;
fork=1:
N
[tmin（k）,fobj（k）]=fminbnd（ft,tt（k）,tt（k+1））;
end
[fobj,ii]=sort（fobj）;%将目标值由小到大排列
tmin=tmin（ii）;%使极小值点做与目标值相应的重新排列
fobj,tmin
（2）最后确定的最小值点
在
的不同分割下，经观察，最后确定出
最小值点是-1.28498111480531
相应目标值是-0.186********545
（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）
gridon,shg
（B）经观察后，把最小值附近邻域放到足够大，然后运行以下指令，那放大图形被推向前台，与此同时光标变为“十字线”，利用它点击极值点可得到最小值数据
[tmin2,fobj2]=ginput
（1）
tmin2=
-1.28500000993975
fobj2=
-0.186********136
出现具有相同数值的刻度区域表明已达最小可分辨状态
（4）符号法求最小值的尝试
symst
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
〖说明〗
●最小值是对整个区间而言的，极小值是对邻域而言的。
●在一个区间中寻找最小值点，对不同子区间分割进行多次搜索是必要的。
这样可以避免把极小值点误作为最小值点。
最小值点是从一系列极小值点和边界点的比较中确定的。
●作图法求最小值点，很直观。
假若绘图时，自变量步长取得足够小，那么所求得的最小值点有相当好的精度。
●符号法在本例中，只求出一个极值点。
其余很多极值点无法秋初，更不可能得到最小值。
16设
，用数值法和符号法求
。
〖目的〗
●学习如何把高阶微分方程写成一阶微分方程组。
●ode45解算器的导数函数如何采用匿名函数形式构成。
●如何从ode45一组数值解点，求指定自变量对应的函数值。
〖解答〗
（1）改写高阶微分方程为一阶微分方程组
令
，于是据高阶微分方程可写出
（2）运行以下指令求y（t）的数值解
formatlong
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）符号法求解
symst;
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函数文件表达，具体如下。
functionS=prob_DyDt（t,y）
S=[y
（2）;-2*y
（1）+3*y
（2）+1];
17已知矩阵A=magic（8），
（1）求该矩阵的“值空间基阵”B；
（2）写出“A的任何列可用基向量线性表出”的验证程序（提示：
利用rref检验）。
〖目的〗
●体验矩阵值空间的基向量组的不唯一性，但它们可以互为线性表出。
●利用rref检验两个矩阵能否互为表出。
〖解答〗
（1）A的值空间的三组不同“基”
A=magic（8）;%采用8阶魔方阵作为实验矩阵
[R,ci]=rref（A）;
B1=A（:
ci）%直接从A中取基向量
B2=orth（A）%求A值空间的正交基
[V,D]=eig（A）;
rv=sum（sum（abs（D））>1000*eps）;%非零特征值数就是矩阵的秩
B3=V（:
1:
rv）%取A的非零特征值对应的特征向量作基
B1=
6423
95554
174746
402627
323435
412322
491514
85859
B2=
-0.35360.54010.3536
-0.3536-0.3858-0.3536
-0.3536-0.2315-0.3536
-0.35360.07720.3536
-0.3536-0.07720.3536
-0.35360.2315-0.3536
-0.35360.3858-0.3536
-0.3536-0.54010.3536
B3=
0.35360.62700.3913
0.3536-0.4815-0.2458
0.3536-0.3361-0.1004
0.35360.1906-0.0451
0.35360.0451-0.1906
0.35360.10040.3361
0.35360.24580.4815
0.3536-0.3913-0.6270
（2）验证A的任何列可用B1线性表出
B1_A=rref（[B1,A]）%若B1_A矩阵的下5行全为0，
%就表明A可以被B1的3根基向量线性表出
B1_A=
10010011001
01001034-3-47
001001-3-445-7
00000000000
00000000000
00000000000
00000000000
00000000000
B2_A=rref（[B2,A]）
B2_A=
Columns1through7
1.000000-91.9239-91.9239-91.9239-91.9239
01.0000051.8459-51.8459-51.845951.8459
001.00009.8995-7.0711-4.24261.4142
0000000
0000000
0000000
0000000
0000000
Columns8through11
-91.9239-91.9239-91.9239-91.9239
51.8459-51.8459-51.845951.8459
-1.41424.24267.0711-9.8995
0000
0000
0000
0000
0000
B3_A=rref（[B3,A]）
B3_A=
Columns1through7
1.00000091.923991.923991.923991.9239
01.0000042.3447-38.1021-33.859429.6168
001.000012.6462-16.8889-21.131525.3741
0000000
0000000
0000000
0000000
0000000
Columns8through11
91.923991.923991.923991.9239
25.3741-21.1315-16.888912.6462
29.6168-33.8594-38.102142.3447
0000
0000
0000
0000
0000
〖说明〗
●magic（n）产生魔方阵。
魔方阵具有很多特异的性质。
就其秩而言，当n为奇数时，该矩阵满秩；当n为4的倍数时，矩阵的秩总是3；当n为偶数但不是4倍数时，则矩阵的秩等于（n/2+2）。
关于魔方阵的有关历史，请见第6.1.3节。
18已知由MATLAB指令创建的矩阵A=gallery（5），试对该矩阵进行特征值分解，并通过验算观察发生的现象。
〖目的〗
●展示特征值分解可能存在的数值问题。
●condeig是比较严谨的特征值分解指令。
●Jordan分解的作用。
〖解答〗
（1）特征值分解
A=gallery（5）
[V,D]=eig（A）;
diag（D）'%为紧凑地显示特征值而写
A=
-911-2163-252
70-69141-4211684
-575575-11493451-13801
3891-38917782-2334593365
1024-10242048-614424572
ans=
Columns1through4
-0.0181-0.0054-0.0171i-0.0054+0.0171i0.0144-0.0104i
Column5
0.0144+0.0104i
（2）验算表明相对误差较大
AE=V*D/V
er_AE=norm（A-AE,'fro'）/norm（A,'fro'）%相对F范数
AE=
1.0e+004*
Columns1through4
-0.0009+0.0000i0.0011-0.0000i-0.0021+0.0000i0.0063-0.0000i
0.0070-0.0000i-0.0069+0.0000i0.0141-0.0000i-0.0421+0.0000i
-0.0575+0.0000i0.0575-0.0000i-0.1149+0.0000i0.3451-0.0000i
0.3891-0.0000i-0.3891+0.0000i0.7781-0.0000i-2.3343+0.0000i
0.1024-0.0000i-0.1024+0.0000i0.2048-0.0000i-0.6144+0.0000i
Column5
-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
（3）一个更严谨的特征值分解指令
[Vc,Dc,eigc]=condeig（A）%eigc中的高值时，说明相应的特征值不可信。
Vc=
Columns1through4
-0.0000-0.0000+0.0000i-0.0000-0.0000i0.0000+0.0000i
0.02060.0207+0.0000i0.0207-0.0000i0.0207+0.0000i
-0.1397-0.1397+0.0000i-0.1397-0.0000i-0.1397+0.0000i
0.95740.95740.95740.9574
0.25190.2519-0.0000i0.2519+0.0000i0.2519-0.0000i
Column5
0.0000-0.0000i
0.0207-0.0000i
-0.1397-0.0000i
0.9574
0.2519+0.0000i
Dc=
Columns1through4
-0.0181000
0-0.0054+0.0171i00
00-0.0054-0.0171i0
0000.0144+0.0104i
0000
Column5
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）;%求出准确的广义特征值，使A*VJ=VJ*D成立。
DJ
AJ=VJ*DJ/VJ
er_AJ=norm（A-AJ,'fro'）/norm（A,'fro'）
DJ=
01000
00100
00010
00001
00000
AJ=
1.0e+004*
-0.00090.0011-0.00210.0063-0.0252
0.0070-0.00690.0141-0.04210.1684
-0.05750.0575-0.11490.3451-1.3801
0.3891-0.38910.7782-2.33459.3365
0.1024-0.10240.2048-0.61442.4572
er_AJ=
2.0500e-011
〖说明〗
●指令condeig的第3输出量eigc给出的是所谓的“矩阵特征值条件数”。
当特征条件数与
相当时，就意味着矩阵A可能“退化”，即矩阵可能存在“代数重数”大于“几何重数”的特征值。
此时，实施Jordan分解更适宜。
●顺便指出：
借助condeig算得的特征值条件数与cond指令算得的矩阵条件数是两个不同概念。
前者描述特征值的问题，后者描述矩阵逆的问题。
●本例矩阵A的特征值条件数很高，表明分解不可信。
验算也表明，相对误差较大。
●当对矩阵A进行Jordan分解时，可以看到，A具有5重根。
当对Jordan分解进行验算时，相对误差很小。
19求矩阵
的解，A为3阶魔方阵，b是
的全1列向量。
〖提示〗
●了解magic指令
●rref用于方程求解。
●矩阵除法和逆阵法解方程。
〖目的〗
●满秩方阵求解的一般过程。
●rref用于方程求解。
●矩阵除法和逆阵法解方程。
〖解答〗
A=magic（3）;%产生3阶魔方阵
b=ones（3,1）;%（3*1）全1列向量
[R,C]=rref（[A,b]）%GaussJordan消去法解方程，同时判断解的唯一性
x=A\b%矩阵除解方程
xx=inv（A）*b%逆阵法解方程
R=
1.0000000.0667
01.000000.0667
001.00000.0667
C=
123
x=
0.0667
0.0667
0.0667
xx=
0.0667
0.0667
0.0667
〖说明〗
●rref指令通过对增广矩阵进行消去法操作完成解方程。
由分解得到的3根“坐标向量”和（或）C3指示的3根基向量，可见A3满秩，因此方程解唯一。
●在本例情况下，矩阵除、逆阵法、rref法所得解一致。
110求矩阵
的解，A为4阶魔方阵，b是
的全1列向量。
〖提示〗
●用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.0000001.00000.0588
01.000003.00000.1176
001.0000-3.0000-0.0588
00000
C=
123
关于以上结果的说明：
●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:
Matrixisclosetosingularorbadlyscaled.
Resultsmaybeinaccurate.RCOND=1.306145e-017.
x=
0.0588
0.1176
-0.0588
0
运行结果指示：
矩阵除法给出的解与rref解相同。
（实际上，MATLAB在设计“除法”程序时，针对“b在A值空间中”的情况，就是用rref求解的。
）
（3）
升级会员 