matlab课后习题答案第四章演示教学.docx
《matlab课后习题答案第四章演示教学.docx》由会员分享,可在线阅读,更多相关《matlab课后习题答案第四章演示教学.docx(34页珍藏版)》请在冰豆网上搜索。
matlab课后习题答案第四章演示教学
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值空